SparseRLDLTSolver.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_SPARSERLDLTSOLVER_H
#define MYRAMATH_MULTIFRONTAL_SPARSERLDLTSOLVER_H

#include <myramath/utility/detail/LIBPUBLIC.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/rldlt/Kernel.h>
#include <myramath/multifrontal/detail/llt/LContainer.h>

// Permutation type.
#include <vector>

namespace myra {

// Forward declarations, components.
class AssemblyTree;
template <class Precision> class RLDLTKernel;

// Forward declarations, A type.
class Permutation;
template<class Number> class SparseMatrix;
template<class Number> class SparseMatrixRange;
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;
template <class Number> class LowerMatrix;
template <class Number> class LowerMatrixRange;
template <class Number> class CLowerMatrixRange;

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

template<class Precision> class LIBPUBLIC SparseRLDLTSolver
  {
  public:

    // Useful typedefs for various dense/sparse ranges.
    typedef Precision Number;
    typedef MatrixRange<Number> DRange;        // Dense range type (e.g. right hand sides)
    typedef LowerMatrixRange<Number> LRange;   // Lower range type (e.g. symmetric schur complement)
    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.

    SparseRLDLTSolver();

    SparseRLDLTSolver(const SparseRLDLTSolver& that);

    void swap(SparseRLDLTSolver& that);

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

    SparseRLDLTSolver& operator = (SparseRLDLTSolver that);

    explicit SparseRLDLTSolver(InputStream& in);

    void write(OutputStream& out) const;

    // ----------------------------------------- Factorization and factorizing constructors.

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

    SparseRLDLTSolver(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());

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

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

    //   side  = Solve by D from the 'L'eft or from the 'R'ight?
    //   Since D is filled with +1/-1's only, there's no op modifier (they'd all do the same thing)
    void solveD(const DRange& B, char side = 'L') 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 symm() 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;
    class SparseRLDLTSolver_RefineAction;
    friend class SparseRLDLTSolver_RefineAction;

    // ----------------------------------------- Calculating schur complements.
 
    LowerMatrix<Number> schur(const SRange& B, Options options = defaults()) const;
    void schur_inplace(const SRange& B, const LRange& S, Options options = defaults()) const;
    JobGraph schur_jobgraph(const SRange& B, const LRange& S, Options options = defaults()) const;

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

    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 ij) const;

    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;

    LowerMatrix<Number> inverse(const iRange& ij, Options options = defaults()) const;
    void inverse_inplace(const iRange& ij, const LRange& Z, Options options = defaults()) const;
    JobGraph inverse_jobgraph(const iRange& ij, const LRange& 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;

    std::pair<int,int> inertia() const;

    const AssemblyTree& tree() const;

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

  private:

    // Numeric contents.    
    typedef RLDLTKernel<Precision> Kernel;
    typedef ::myra::multifrontal::detail::llt::LContainer<Kernel> LContainer;
    LContainer L;

  };

template<class Precision> class ReflectNumber <SparseRLDLTSolver<Precision> >
  { public: typedef Precision type; };

LIBPUBLIC Vector<NumberS> operator * (const SparseRLDLTSolver<NumberS>& solver, const CVectorRange<NumberS>& x);
LIBPUBLIC Vector<NumberD> operator * (const SparseRLDLTSolver<NumberD>& solver, const CVectorRange<NumberD>& x);

LIBPUBLIC Matrix<NumberS> operator * (const SparseRLDLTSolver<NumberS>& solver, const CMatrixRange<NumberS>& X);
LIBPUBLIC Matrix<NumberD> operator * (const SparseRLDLTSolver<NumberD>& solver, const CMatrixRange<NumberD>& X);

} // namespace myra

#endif