SparseLUSolver.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_SPARSELUSOLVER_H
#define MYRAMATH_MULTIFRONTAL_SPARSELUSOLVER_H

#include <myramath/MYRAMATH_EXPORT.h>

// For ReflectNumber<>
#include <myramath/utility/Number.h>

// Return type for all _jobgraph() methods.
#include <myramath/jobgraph/JobGraph.h>

// Options pack.
#include <myramath/multifrontal/Options.h>

// Underlying numeric storage.
#include <myramath/multifrontal/lu/Kernel.h>
#include <myramath/multifrontal/detail/lu/LUContainer.h>

// Permutation type.
#include <vector>

namespace myra {

// Forward declarations, components.
class AssemblyTree;
template<class Number> class LUKernel;

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

// Forward declarations, X/B types.
class intCRange;
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;

template<class Number> class MYRAMATH_EXPORT SparseLUSolver
  {
  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 intCRange iRange;                  // Indices range type (e.g. stencils for partial solve etc)

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

    // ----------------------------------------- Construction, serialization, value semantics.

    SparseLUSolver();

    SparseLUSolver(const SparseLUSolver& that);

    void swap(SparseLUSolver& that);

#ifdef MYRAMATH_ENABLE_CPP11
    SparseLUSolver(SparseLUSolver&& that);
#endif

    SparseLUSolver& operator = (SparseLUSolver that);

    explicit SparseLUSolver(InputStream& in);

    void write(OutputStream& out) const;

   ~SparseLUSolver();
   
    // ----------------------------------------- Factorization and factorizing constructors.

    SparseLUSolver(const SRange& in_A, Options options = defaults());
    void factor(const SRange& A, Options options = defaults());
    JobGraph factor_jobgraph(const SRange& A, Options options = defaults());

    SparseLUSolver(const SRange& in_A, const Permutation& P, Options options = defaults());
    void factor(const SRange& in_A, const Permutation& P, Options options = defaults());
    JobGraph factor_jobgraph(const SRange& in_A, const Permutation& P, Options options = defaults());

    SparseLUSolver(const SRange& in_A, const AssemblyTree& tree, Options options = defaults());
    void factor(const SRange& in_A, const AssemblyTree& tree, Options options = defaults());
    JobGraph factor_jobgraph(const SRange& in_A, const AssemblyTree& tree, Options options = defaults());

    // ----------------------------------------- Linear solution.

    //    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()) const;
    JobGraph solve_jobgraph(const DRange& B, char side = 'L', char op = 'N', Options options = defaults()) const;

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

    //   side  = Solve by U from the 'L'eft or from the 'R'ight?
    //   op    = Apply an operation to U? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    void solveU(const DRange& B, char side, char op, Options options) const;
    JobGraph solveU_jobgraph(const DRange& B, char side, char op, Options options) const;

    // ----------------------------------------- Solution with refinement.

    //    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 sparse gemm() too.
    //  Returns the residual history, frobenius(A*X-B), evaluated at each refinement iteration.
    std::vector<Precision> refine(const DRange& B, char side = 'L', char op = 'N', Precision tolerance = default_tolerance(), int iterations = default_iterations(), Options options = defaults()) const;
    class SparseLUSolver_RefineAction;
    friend class SparseLUSolver_RefineAction;

    // ----------------------------------------- Calculating schur complements.

    // Note that in the complex case, ' denotes (non-conjugated) 'T'ranspose, not 'H'ermitian.
    Matrix<Number> schur(const SRange& B, const SRange& C, Options options = defaults()) const;
    void schur_inplace(const SRange& B, const SRange& C, const DRange& S, Options options = defaults()) const;
    JobGraph schur_jobgraph(const SRange& B, const SRange& C, const DRange& S, Options options = defaults()) const;    

    Matrix<Number> schur(const iRange& Bi, const DRange& Bv, const iRange& Ci, const DRange& Cv, Options options = defaults()) const;    
    void schur_inplace(const iRange& Bi, const DRange& Bv, const iRange& Ci, const DRange& Cv, const DRange& S, Options options = defaults()) const;
    JobGraph schur_jobgraph(const iRange& Bi, const DRange& Bv, const iRange& Ci, const DRange& Cv, const DRange& S, Options options = defaults()) const;

    // ----------------------------------------- Sampling inv(A).

    Number inverse(int i, int j) const;

    Matrix<Number> inverse(const iRange& i, const iRange& j, Options options = defaults()) const;
    void inverse_inplace(const iRange& i, const iRange& j, const DRange& Z, Options options = defaults()) const;
    JobGraph inverse_jobgraph(const iRange& i, const iRange& j, const DRange& Z, Options options = defaults()) const;

    // ----------------------------------------- Partial linear solution.

    // Note X = partialsolve(linspace(0,N),linspace(0,N),B,side,op) yields the same result X = solve(B,side,op)
    // Note X = partialsolve(i,j,B,'L'eft, op) yields the same result as X =   op(this->inverse(i,j))*B
    // Note X = partialsolve(i,j,B,'R'ight,op) yields the same result as X = B*op(this->inverse(i,j))
    Matrix<Number> partialsolve(const iRange& i, const iRange& j, const CMatrixRange<Number>& B, char side = 'L', char op = 'N', Options options = defaults().set_nthreads(1)) const;
    void partialsolve_inplace(const iRange& i, const iRange& j, const CMatrixRange<Number>& B, const MatrixRange<Number>& X, char side = 'L', char op = 'N', Options options = defaults().set_nthreads(1)) const;
    JobGraph partialsolve_jobgraph(const iRange& i, const iRange& j, const CMatrixRange<Number>& B, const MatrixRange<Number>& X, char side = 'L', char op = 'N', Options options = defaults().set_nthreads(1)) const;    

    // ----------------------------------------- Miscellany.

    int size() const;

    const AssemblyTree& tree() const;

    static Options defaults();
    static Precision default_tolerance();
    static int default_iterations();

  private:

    // Numeric contents.
    typedef LUKernel<Number> Kernel;
    typedef ::myra::multifrontal::detail::lu::LUContainer<Kernel> LUContainer;
    LUContainer LU;

  };

template<class Number> class ReflectNumber <SparseLUSolver<Number> >
  { public: typedef Number type; };

MYRAMATH_EXPORT Vector<NumberS> operator * (const SparseLUSolver<NumberS>& solver, const CVectorRange<NumberS>& x);
MYRAMATH_EXPORT Vector<NumberD> operator * (const SparseLUSolver<NumberD>& solver, const CVectorRange<NumberD>& x);
MYRAMATH_EXPORT Vector<NumberC> operator * (const SparseLUSolver<NumberC>& solver, const CVectorRange<NumberC>& x);
MYRAMATH_EXPORT Vector<NumberZ> operator * (const SparseLUSolver<NumberZ>& solver, const CVectorRange<NumberZ>& x);

MYRAMATH_EXPORT Matrix<NumberS> operator * (const SparseLUSolver<NumberS>& solver, const CMatrixRange<NumberS>& X);
MYRAMATH_EXPORT Matrix<NumberD> operator * (const SparseLUSolver<NumberD>& solver, const CMatrixRange<NumberD>& X);
MYRAMATH_EXPORT Matrix<NumberC> operator * (const SparseLUSolver<NumberC>& solver, const CMatrixRange<NumberC>& X);
MYRAMATH_EXPORT Matrix<NumberZ> operator * (const SparseLUSolver<NumberZ>& solver, const CMatrixRange<NumberZ>& X);

} // namespace myra

#endif