SparseNormalSolver.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_MULTIFRONTAL_SPARSENORMALSOLVER_H
#define MYRAMATH_MULTIFRONTAL_SPARSENORMALSOLVER_H

#include <myramath/MYRAMATH_EXPORT.h>
#include <myramath/utility/Number.h>

#include <myramath/sparse/SparseMatrix.h>

#include <myramath/multifrontal/Options.h>
#include <myramath/multifrontal/rcholesky/SparseRCholeskySolver.h>
#include <myramath/multifrontal/zcholesky/SparseZCholeskySolver.h>
#include <myramath/multifrontal/rldlt/SparseRLDLTSolver.h>
#include <myramath/multifrontal/zldlh/SparseZLDLHSolver.h>

namespace myra {

// Forward declarations, A type.
template<class Number> class SparseMatrixRange;
template<class Number> class CSparseMatrixRange;

// Forward declarations, X/B types.
template <class Number> class Matrix;
template <class Number> class MatrixRange;
template <class Number> class CMatrixRange;
template <class Number> class Vector;
template <class Number> class VectorRange;
template <class Number> class CVectorRange;

// Forward declarations, serialization.
class InputStream;
class OutputStream;

// Reflects the underlying cholesky-like solver to use.
template<class Number> class ReflectNormal;
/*
template<> class ReflectNormal<float>
  { public: typedef SparseRCholeskySolver<float> type; };
template<> class ReflectNormal<double>
  { public: typedef SparseRCholeskySolver<double> type; };
template<> class ReflectNormal<std::complex<float> >
  { public: typedef SparseZCholeskySolver<float> type; };
template<> class ReflectNormal<std::complex<double> >
  { public: typedef SparseZCholeskySolver<double> type; };
*/
template<> class ReflectNormal<float>
  { public: typedef SparseRLDLTSolver<float> type; };
template<> class ReflectNormal<double>
  { public: typedef SparseRLDLTSolver<double> type; };
template<> class ReflectNormal<std::complex<float> >
  { public: typedef SparseZLDLHSolver<float> type; };
template<> class ReflectNormal<std::complex<double> >
  { public: typedef SparseZLDLHSolver<double> type; };

// Sparse direct solver for general systems, solves the normal equations.
template<class Number> class MYRAMATH_EXPORT SparseNormalSolver
  {
  public:

    // Useful typedefs for various dense/sparse ranges.
    typedef typename ReflectPrecision<Number>::type Precision;
    typedef MatrixRange<Number> DRange;        // Dense range type (e.g. right hand sides)
    typedef CSparseMatrixRange<Number> SRange; // Sparse range type (e.g. underlying A)

    // Typedef for Options pack.
    typedef ::myra::multifrontal::Options Options;

    SparseNormalSolver(const SRange& in_A, Options options = defaults());

    explicit SparseNormalSolver(InputStream& in);

    void write(OutputStream& out) const;

   ~SparseNormalSolver();

    int size() const;

    //    side  = Solve by A from the 'L'eft or from the 'R'ight?
    //    op    = Apply an operation to A? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    void solve(const DRange& B, char side = 'L', char op = 'N', Options options = defaults().set_nthreads(1)) const;

    //    side  = Solve by A from the 'L'eft or from the 'R'ight?
    //    op    = Apply an operation to A? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    //  Note, this is more expensive because it calls solve() multiple times, and spmm() too.
    //  Returns the history of the relative residual, frobenius(A*X-B)/frobenius(B)
    std::vector<Precision> refine(const DRange& B, char side = 'L', char op = 'N', Precision tolerance = default_tolerance(), int iterations = default_iterations(), Options options = defaults().set_nthreads(1)) const;

  private:

    // Returns default options.
    static Options defaults();

    // Default parameters for refine().
    static Precision default_tolerance();
    static int default_iterations();

    // Holds forward operator A.
    SparseMatrix<Number> A;
    SparseMatrix<Number> AA;
    
    // Underlying cholesky-like solver.
    typedef typename ReflectNormal<Number>::type Solver;
    Solver solver;
    
  };

// ReflectNumber for SparseNormalSolver
template<class Number> class ReflectNumber <SparseNormalSolver<Number> >
  { public: typedef Number type; };

} // namespace myra

#endif