Installation of external librairies

SelectDirectSolver

Syntax :

  void SelectDirectSolver(int type);

This method sets the Cholesky solver to use, you can choose between :

CHOLMOD
PASTIX
SELDON_SOLVER : Basic sparse solver proposed by Seldon (very slow)

Example :

// you fill a symmetric sparse matrix (assumed positive definite)
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you declare a sparse Cholesky Solver
SparseCholeskySolver<double> mat_chol;
// you select which solver you want to use
mat_chol.SelectDirectSolver(mat_chol.PASTIX);
// then you can factorize with this solver
mat_chol.Factorize(A);

Location :

SetNonSymmetricIlut

Syntax :

  void SetNonSymmetricIlut();

This method tells to use non-symmetric incomplete factorization if ILUT is selected as a direct solver.

Example :

// you fill a sparse matrix
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// you select which solver you want to use
mat_lu.SelectDirectSolver(mat_lu.ILUT);
// for ILUT solver you can use non-symmetric version even though the matrix is symmetric
mat_lu.SetNonSymmetricIlut();
// then you can factorize with this solver
mat_lu.Factorize(A);

Location :

SetThresholdMatrix, GetThresholdMatrix

Syntax :

  void SetThresholdMatrix(double);
  double GetThresholdMatrix() const;

The method SetThresholdMatrix sets the threshold used to drop values during the incomplete factorization. This method is useful only if ILUT is selected as a direct solver. The method GetThresholdMatrix returns the threshold currently set.

Example :

// you fill a sparse matrix
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// you select which solver you want to use
mat_lu.SelectDirectSolver(mat_lu.ILUT);
// for ILUT solver you can use non-symmetric version even though the matrix is symmetric
mat_lu.SetNonSymmetricIlut();
// dropping threshold
mat_lu.SetThresholdMatrix(1e-2);
// which epsilon has been set ?
double eps = mat_lu.GetThresholdMatrix();
// then you can factorize with this solver
mat_lu.Factorize(A);

Location :

IsAvailableSolver

Syntax :

  bool IsAvailableSolver(int);

This method returns true if the solver is available or not (i.e. if you have compiled with this solver).

Example :

// you fill a sparse matrix
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);

// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;

// you can check if your favorite solver is present
if (mat_lu.IsAvailableSolver(mat_lu.MUMPS))
  {
    // you select which solver you want to use
    mat_lu.SelectDirectSolver(mat_lu.MUMPS);

    // then you can factorize with this solver
    mat_lu.Factorize(A);
  }
else
  cout << "recompile with Mumps" << endl;

Location :

SparseDirectSolver.cxx
SparseCholeskyFactorisation.cxx

GetDirectSolver

Syntax :

  int GetDirectSolver();

This method returns the direct solver used.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// which solver used by default ?
int type_solver = mat_lu.GetDirectSolver();
if (type_solver == mat_lu.SELDON_SOLVER)
  cout << "Warning : this solver is slow" << endl;

Location :

GetInfoFactorization

Syntax :

  void GetInfoFactorization();

This method returns the error code provided by the used direct solver. For SparseDirectSolver, the error codes are listed in an enum attribute, and can be :

FACTO_OK : successful factorization (= 0 )
STRUCTURALLY_SINGULAR_MATRIX : the matrix is structurally singular. It is probably because there is an empty row or column.
NUMERICALLY_SINGULAR_MATRIX : the matrix is numerically singular. It happens when a null pivot has been found, it may occur if the condition number of the matrix is very large.
OUT_OF_MEMORY : there is not enough memory to complete the factorization.
INVALID_ARGUMENT : the arguments provided to the direct solver are not correct.
INCORRECT_NUMBER_OF_ROWS : the number of rows specified is incorrect (usually negative or null)
MATRIX_INDICES_INCORRECT : the indices are incorrect (usually out of range, i.e. negative or greater than the dimensions of the matrix)
INVALID_PERMUTATION : the permutation array is not valid (probably, not defining a bijection).
ORDERING_FAILED : the computation of an appropriate ordering (by Metis, Scotch, etc) has failed
INTERNAL_ERROR : unknown error, you should look at the documentation of the used solver
OVERFLOW_32BIT : you should use 64-bits integers to avoid overflow of integers

For specific solvers (MatrixMumps, MatrixUmfPack, etc), the error code is described in the documentation of each solver.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// you factorize the matrix
mat_lu.Factorize(A);
// to know if there is a problem
int info = mat_lu.GetInfoFactorization();
if (info == mat_lu.OUT_OF_MEMORY)
  {
    cout << "The matrix is too large, not enough memory" << endl;
  }

// if you are using directly a given solver
MatrixMumps<double> mat_mumps;
mat_mumps.Factorize(A);
int info = mat_mumps.GetInfoFactorization();
// info is described in Mumps documentation
if (info != 0)
  {
    cout << "MUMPS Error = " << info << endl;
  }

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx
Pastix.cxx
SparseDirectSolver.cxx

UseInteger8

Syntax :

  bool UseInteger8 const();

This method returns true if the current solver is using 64-bits integers (integer*8 in Fortran).

Example :

// you declare a sparse Solver
MatrixPastix<double> mat_lu;

// you can check if the current version is using 64-bits integers or not
bool use_int8 = mat_lu.UseInteger8();

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx
Pastix.cxx
SparseSolver.cxx

Factorize

Syntax :

  void Factorize(Matrix&);
  void Factorize(Matrix&, bool);

This method factorizes the given matrix. In parallel execution, this method will consider that the linear system to solve is restricted to the current processor. For example, you can use this method to factorize independant linear systems. If the matrix is distributed overall severall processors, you should call FactorizeDistributed or use the class SparseDistributedSolver.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// you factorize the matrix
mat_lu.Factorize(A);
// to know if there is a problem
int info = mat_lu.GetInfoFactorization();
if (info == mat_lu.OUT_OF_MEMORY)
  {
    cout << "The matrix is too large, not enough memory" << endl;
  }
// once the matrix is factorized, you can solve systems
Vector<double> x(n);
mat_lu.Solve(x);

Location :

Factorize

Syntax :

  void Factorize(Matrix&);
  void Factorize(Matrix&, bool);

Example :

// you fill a sparse matrix
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you declare a sparse Solver
SparseCholeskySolver<double> mat_chol;
// you factorize the matrix
mat_chol.Factorize(A);

// once the matrix is factorized, you can solve systems
Vector<double> x(n); x.FillRand();
mat_col.Solve(SeldonNoTrans, x);
mat_col.Solve(SeldonTrans, x);

Location :

Solve

Syntax :

  void Solve(Vector&);
  void Solve(SeldonTranspose, Vector&);

  void Solve(Matrix<T, General, ColMajor>&);
  void Solve(SeldonTranspose, Matrix<T, General, ColMajor>&);

This method computes the solution of A x = b or A^T x = b, assuming that Factorize has been called before.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a sparse Solver
SparseDirectSolver<double> mat_lu;
// you factorize the matrix
mat_lu.Factorize(A);
// once the matrix is factorized, you can solve systems
Vector<double> x(n), b(n);
b.Fill();
x = b;
mat_lu.Solve(x);
// to solve A^T x = b :
x = b;
mat_lu.Solve(SeldonTrans, x);

// you can solve system with multiple right hand sides
Matrix<double, General, ColMajor> B(n, 10);
B.FillRand();
// B is overwritten by the solution
mat_lu.Solve(B);

// and transpose system can be solved with multiple right hand sides
B.FillRand();
mat_lu.Solve(SeldonTrans, B);

Location :

Solve

Syntax :

  void Solve(SeldonTranspose, Vector&);

This method computes the solution of L x = b or L^T x = b, assuming that Factorize has been called before.

Example :

// you fill a symmetric sparse matrix
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you declare a sparse Cholesky Solver
SparseCholeskySolver<double> mat_chol;
// you factorize the matrix
mat_chol.Factorize(A);
// once the matrix is factorized, you can solve systems
// for Cholesky factorisation, it will solve either L x = b or L^T x = b
// so you need to call the method twice to recover the solution (since A = L L^T)
x = b;
mat_chol.Solve(SeldonNoTrans, x);
mat_chol.Solve(SeldonTrans, x);

Location :

Mlt

Syntax :

  void Mlt(SeldonTranspose, Vector&);

This method computes the matrix-vector product L x or L^T x, once the Cholesky factorization (A = L L^T) has been computed.

Example :

// for Cholesky factorisation, it will
// compute either L x or L^T x
// the result is overwritten in x

// filling A
A.Reallocate(n, n);

// factorisation of A
SparseCholeskySolver<double> mat_chol;
mat_chol.Factorize(A);

// computation of y = A x
Vector<double> x(n), y(n);
y = x;
mat_chol.Mlt(SeldonTrans, y);
mat_chol.Mlt(SeldonNoTrans, y);

Location :

Mlt

Syntax :

  void Mlt(const SeldonTranspose, Vector&);

This method computes the matrix-vector product L x or L^T x, once the Cholesky factorization (A = L L^T) has been computed.

Example :

// for Cholesky factorisation, it will
// compute either L x or L^T x
// the result is overwritten in x

// filling A
A.Reallocate(n, n);

// factorisation of A
DistributedCholeskySolver<double> mat_chol;
mat_chol.Factorize(A);

// computation of y = A x
Vector<double> x(n), y(n);
y = x;
mat_chol.Mlt(SeldonTrans, y);
mat_chol.Mlt(SeldonNoTrans, y);

Location :

Mumps.cxx
Pastix.cxx
SparseDirectSolver.cxx

FactorizeDistributed, FactorizeDistributedMatrix

Syntax :

  void FactorizeDistributed(MPI::Comm& comm, Vector& Ptr,
                            Vector& Ind, Vector& Val, IVect& glob_num,
                            bool sym, bool keep_matrix);

  void FactorizeDistributedMatrix(MPI::Comm& comm, Vector& Ptr,
                                  Vector& Ind, Vector& Val, IVect& glob_num,
                                  bool sym);

This method factorizes a distributed matrix, which is given through arrays Ptr, Ind, Val (CSC format). glob_num is a local to global array (glob_num(i) is the global row number of local row i). Each column of the global system is assumed to be distributed to one processor and only one. Each processor is assumed to have at least one column affected to its-self. Arrays Ptr and Ind may consist of 64-bit integers (in order to be compatible with 64-bit versions of direct solvers). If sym is true, the matrix is symmetric, and we assume that Ptr, Ind, Val are representing the lower part of the matrix. If sym is false, the matrix is considered non-symmetric. The communicator given in argument regroup all the processors involved in the factorization. This method is working only if you have selected Pastix, Mumps or SuperLU (SuperLU_DIST has been compiled and used). We recommend to not call this function directly and use the class SparseDistributedSolver that will call this function.

Example :

  // initialization of MPI_Init_thread needed if Pastix has been compiled
  // with threads, otherwise you can use MPI_Init
  int required = MPI_THREAD_MULTIPLE;
  int provided = -1;
  MPI_Init_thread(&argc, &argv, required, &provided);
    
  // declaration of the sparse solve
  SparseDirectSolver<double> mat_lu;
  
  // considered global matrix
  // A = |1.5  1    0    0  -0.3 |
  //     | 1  2.0   0  -1.0   0  |
  //     | 0   0   2.0  0    0.8 |
  //     | 0  -1.0  0   3.0  1.2 |
  //     | -0.3  0.0  0.8  1.2  3.0 |
  
  // global right hand side (solution is 1, 2, 3, 4, 5)
  // B = |2 1 10 16 21.9|

  if (MPI::COMM_WORLD.Get_rank() == 0)
    {
      // on first processor, we provide columns 1, 2 and 5
      // only lower part since matrix is symmetric
      int n = 3;
      Matrix<double, General, ArrayColSparse> A(5, n); 
      A(0, 0) = 1.5; A(1,0) = 1.0; A(4,0) = -0.3;
      A(1,1) = 2.0; A(3,1) = -1.0; A(4,2) = 3.0;
      Vector<double> b(n);
      b(0) = 2; b(1) = 1; b(2) = 21.9;
      IVect num_loc(n); num_loc(0) = 0; num_loc(1) = 1; num_loc(2) = 4;
      
      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);
      
      // factorisation
      mat_lu.FactorizeDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                  num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      
      DISP(b);
    }
  else
    {
      // on second processor, we provide columns 3 and 4
      // only lower part since matrix is symmetric
      int n = 2;
      Matrix<double, General, ArrayColSparse> A(5, n);
      A(2, 0) = 2.0; A(4, 0) = 0.8; A(3,1) = 3.0; A(4,1) = 1.2;
      Vector<double> b(n);
      b(0) = 10; b(1) = 16;
      IVect num_loc(n);
      num_loc(0) = 2; num_loc(1) = 3;

      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);

      // factorisation      
      mat_lu.FactorizeDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                  num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      
      DISP(b);
    }
  
  MPI_Finalize();

Location :

PerformAnalysisDistributed, PerformFactorizationDistributed

Syntax :

  void PerformAnalysisDistributed(MPI::Comm& comm, Vector& Ptr,
                                  Vector& Ind, Vector& Val, IVect& glob_num,
                                  bool sym, bool keep_matrix);

  void PerformFactorizationDistributed(MPI::Comm& comm, Vector& Ptr,
                                       Vector& Ind, Vector& Val, IVect& glob_num,
                                       bool sym, bool keep_matrix);

Example :

  // initialization of MPI_Init_thread needed if Pastix has been compiled
  // with threads, otherwise you can use MPI_Init
  int required = MPI_THREAD_MULTIPLE;
  int provided = -1;
  MPI_Init_thread(&argc, &argv, required, &provided);
    
  // declaration of the sparse solve
  SparseDirectSolver<double> mat_lu;
  
  // considered global matrix
  // A = |1.5  1    0    0  -0.3 |
  //     | 1  2.0   0  -1.0   0  |
  //     | 0   0   2.0  0    0.8 |
  //     | 0  -1.0  0   3.0  1.2 |
  //     | -0.3  0.0  0.8  1.2  3.0 |
  
  // global right hand side (solution is 1, 2, 3, 4, 5)
  // B = |2 1 10 16 21.9|

  if (MPI::COMM_WORLD.Get_rank() == 0)
    {
      // on first processor, we provide columns 1, 2 and 5
      // only lower part since matrix is symmetric
      int n = 3;
      Matrix<double, General, ArrayColSparse> A(5, n); 
      A(0, 0) = 1.5; A(1,0) = 1.0; A(4,0) = -0.3;
      A(1,1) = 2.0; A(3,1) = -1.0; A(4,2) = 3.0;
      Vector<double> b(n);
      b(0) = 2; b(1) = 1; b(2) = 21.9;
      IVect num_loc(n); num_loc(0) = 0; num_loc(1) = 1; num_loc(2) = 4;
      
      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);

      // analysis
      mat_lu.PerformAnalysisDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                        num_loc, true);

      // factorisation
      mat_lu.PerformFactorizationDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                             num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      DISP(b);

      // if the pattern is the same as previously, but the values are different
      // only numerical factorisation need to be called
      Val.FillRand();
      mat_lu.PerformFactorizationDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                             num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);

      DISP(b);
    }
  else
    {
      // on second processor, we provide columns 3 and 4
      // only lower part since matrix is symmetric
      int n = 2;
      Matrix<double, General, ArrayColSparse> A(5, n);
      A(2, 0) = 2.0; A(4, 0) = 0.8; A(3,1) = 3.0; A(4,1) = 1.2;
      Vector<double> b(n);
      b(0) = 10; b(1) = 16;
      IVect num_loc(n);
      num_loc(0) = 2; num_loc(1) = 3;

      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);

      // analysis
      mat_lu.PerformAnalysisDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                        num_loc, true);

      // factorisation
      mat_lu.PerformFactorizationDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                             num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      DISP(b);

      // if the pattern is the same as previously, but the values are different
      // only numerical factorisation need to be called
      Val.FillRand();
      mat_lu.PerformFactorizationDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                             num_loc, true);
      
      // then solution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);

      DISP(b);
    }
  
  MPI_Finalize();

Location :

Mumps.cxx
Pastix.cxx
SparseDirectSolver.cxx

SolveDistributed

Syntax :

  void SolveDistributed(MPI::Comm& comm, Vector& x,
                        IVect& glob_num);

  void SolveDistributed(SeldonTrans, MPI::Comm& comm, Vector& x,
                        IVect& glob_num);

This method solves distributed linear system (or its transpose), once FactorizeDistributed has been called. This method is working only if you have selected Pastix or Mumps.

Example :

  // initialization of MPI_Init_thread needed if Pastix has been compiled
  // with threads, otherwise you can use MPI_Init
  int required = MPI_THREAD_MULTIPLE;
  int provided = -1;
  MPI_Init_thread(&argc, &argv, required, &provided);
    
  // declaration of the sparse solve
  SparseDirectSolver<double> mat_lu;
  
  // considered global matrix
  // A = |1.5  1    0    0  -0.3 |
  //     | 1  2.0   0  -1.0   0  |
  //     | 0   0   2.0  0    0.8 |
  //     | 0  -1.0  0   3.0  1.2 |
  //     | -0.3  0.0  0.8  1.2  3.0 |
  
  // global right hand side (solution is 1, 2, 3, 4, 5)
  // B = |2 1 10 16 21.9|

  if (MPI::COMM_WORLD.Get_rank() == 0)
    {
      // on first processor, we provide columns 1, 2 and 5
      // only lower part since matrix is symmetric
      int n = 3;
      Matrix<double, General, ArrayColSparse> A(5, n); 
      A(0, 0) = 1.5; A(1,0) = 1.0; A(4,0) = -0.3;
      A(1,1) = 2.0; A(3,1) = -1.0; A(4,2) = 3.0;
      Vector<double> b(n);
      b(0) = 2; b(1) = 1; b(2) = 21.9;
      IVect num_loc(n); num_loc(0) = 0; num_loc(1) = 1; num_loc(2) = 4;
      
      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);
      
      // factorization
      mat_lu.FactorizeDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                  num_loc, true);
      
      // then resolution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      
      DISP(b);
    }
  else
    {
      // on second processor, we provide columns 3 and 4
      // only lower part since matrix is symmetric
      int n = 2;
      Matrix<double, General, ArrayColSparse> A(5, n);
      A(2, 0) = 2.0; A(4, 0) = 0.8; A(3,1) = 3.0; A(4,1) = 1.2;
      Vector<double> b(n);
      b(0) = 10; b(1) = 16;
      IVect num_loc(n);
      num_loc(0) = 2; num_loc(1) = 3;

      // converting to csc format
      IVect Ptr, IndRow; Vector<double> Val;
      General prop;
      ConvertToCSC(A, prop, Ptr, IndRow, Val);

      // factorization
      mat_lu.FactorizeDistributed(MPI::COMM_WORLD, Ptr, IndRow, Val,
                                  num_loc, true);
      
      // then resolution
      mat_lu.SolveDistributed(MPI::COMM_WORLD, SeldonNoTrans, b, num_loc);
      
      DISP(b);
    }
  
  MPI_Finalize();

Location :

Mumps.cxx
Pastix.cxx
SparseDirectSolver.cxx

Clear

Syntax :

  void Clear();

This method releases the memory used by the factorization. It is available for every direct solver.

Example :

Matrix<double, General, ArrayRowSparse> A;
MatrixUmfPack<double> mat_lu;
// you fill A as you want
// then a first factorization is achieved
GetLU(A, mat_lu);
// then solve needed linear systems
SolveLU(mat_lu, x);
// and if you need to spare memory, you can clear factorization
mat_lu.Clear();

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx
Pastix.cxx
SparseDirectSolver.cxx
SparseCholeskyFactorisation.cxx

SetPrintLevel

Syntax :

  void SetPrintLevel(int level);

This method sets the level of verbosity. A level equal to 0 means that no messages will be displayed. A level equal to 2 will induce moderate display while a level equal to 6 will display more details.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);

// you declare a sparse Solver
SparseDistributedSolver<double> mat_lu;

// you select your print level
mat_lu.SetPrintLevel(2);

// then you can factorize
mat_lu.Factorize(A);

Location :

Factorize

Syntax :

  void Factorize(Matrix&, bool keep_matrix = false, bool scale_matrix = false);
  void Factorize(DistributedMatrix&, bool keep_matrix = false, bool scale_matrix = false);

This method factorizes a sequential or distributed matrix. The second argument is optional, if equal to true the input matrix is kept, otherwise the matrix is erased. The third argument is optional, if equal to true the matrix is scaled before factorizing it.

Example :

// you fill a distributed matrix
DistributedMatrix<double, General, ArrayRowSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);

// you declare a sparse Solver
SparseDistributedSolver<double> mat_lu;

// then you can factorize
mat_lu.Factorize(A);

// and solve as many linear systems you want
// x contains the right hand side on input, the solution on output
mat_lu.Solve(x);

Location :

DistributedCholeskySolver.cxx

Factorize

Syntax :

  void Factorize(Matrix&, bool keep_matrix = false);
  void Factorize(DistributedMatrix&, bool keep_matrix = false);

This method factorizes a sequential or distributed matrix. The second argument is optional, if equal to true the input matrix is kept, otherwise the matrix is erased.

Example :

// you fill a distributed matrix
DistributedMatrix<double, Symmetric, ArrayRowSymSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);

// you declare a sparse Solver
DistributedCholeskySolver<double> mat_lu;

// then you can factorize
mat_lu.Factorize(A);

// and solve as many linear systems you want
// x contains the right hand side on input, the solution on output
mat_lu.Solve(SeldonNoTrans, x);
mat_lu.Solve(SeldonTrans, x);

Location :

PerformAnalysis, PerformFactorization

Syntax :

  void PerformAnalysis(Matrix&);
  void PerformFactorization(Matrix&, bool scale_matrix = false);

  void PerformAnalysis(DistributedMatrix&);
  void PerformFactorization(DistributedMatrix&, bool scale_matrix = false);

This method performs either the symbolic or numerical factorization of a sequential or distributed matrix. The second argument is optional, if equal to true the matrix is scaled before factorizing it.

Example :

// you fill a distributed matrix
DistributedMatrix<double, General, ArrayRowSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);

// you declare a sparse Solver
SparseDistributedSolver<double> mat_lu;

// then you can perform a symbolic factorization
// the matrix is not erased
mat_lu.PerformAnalysis(A);

// PerformFactorization completes the numerical factorization
// and erases the input matrix
mat_lu.PerformFactorization(A);

// and solve as many linear systems as you want
// x contains the right hand side on input, the solution on output
mat_lu.Solve(x);

// then you can reconstruct A
A.Reallocate(n, n);

// fill it again
...

// if the pattern is the same, you do not need to call PerformAnalysis
mat_lu.PerformFactorization(A);

mat_lu.Solve(x);

Location :

Solve

Syntax :

  void Solve(Vector& x, const Vector& b);
  void Solve(Vector& x);
  void Solve(SeldonTranspose, Vector& x);

  void Solve(Matrix& x);
  void Solve(SeldonTranspose, Matrix& x);

This method solves a sequential or distributed linear system. For multiple right hand sides, we consider they are stored in a Matrix with storage ColMajor, each column being a different right hand side.

Example :

// you fill a distributed matrix
DistributedMatrix<double, General, ArrayRowSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);

// you declare a sparse Solver
SparseDistributedSolver<double> mat_lu;

// then you can factorize
mat_lu.Factorize(A);

// and solve as many linear systems you want
Vector<double> x(n), b(n);
b.FillRand();

// first syntax, the right hand side b and the solution x are both provided
mat_lu.Solve(x, b);

// second syntax : x contains the right hand side on input, the solution on output
x = b;
mat_lu.Solve(x);

// third syntax : you can ask to solve A^T x = b
x = b;
mat_lu.Solve(SeldonTrans, x);

// for multiple right hand sides, they are rearranged in a matrix
// each column is a different right hand side
int nb_rhs = 10;
Matrix<double, General, ColMajor> X(n, nb_rhs);
X.FillRand();

mat_lu.Solve(X);

// you can ask to solve transpose system as well
mat_lu.Solve(SeldonTrans, X);

Location :

DistributedCholeskySolver.cxx

Solve

Syntax :

  void Solve(SeldonTranspose, Vector& x);

This method solves a sequential or distributed linear system. On input, x contains the right hand side, on ouput the solution.

Example :

// you fill a distributed matrix
DistributedMatrix<double, Symmetric, ArrayRowSymSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);
// filling the matrix ...

// you declare a sparse Solver
DistributedCholeskySolver<double> mat_lu;

// then you can factorize
mat_lu.Factorize(A);

// Cholesky factorisation A = L L^T
// Solve will compute the solution of L x = y
Vector<double> y(n), x;
y.FillRand(); x = y;
mat_lu.Solve(SeldonNoTrans, x);

// or L^t x = y
x = y;
mat_lu.Solve(SeldonTrans, x);

Location :

TransSolve

Syntax :

  void TransSolve(Vector& x);
  void TransSolve(Matrix& x);

This method solves the transpose of a sequential or distributed linear system. For multiple right hand sides, we consider they are stored in a Matrix with storage ColMajor, each column being a different right hand side.

Example :

// you fill a distributed matrix
DistributedMatrix<double, General, ArrayRowSparse> A(n, n);

// informations about distributed rows/columns through function Init
A.Init(...);

// you declare a sparse Solver
SparseDistributedSolver<double> mat_lu;

// then you can factorize
mat_lu.Factorize(A);

// and solve as many linear systems you want
Vector<double> x(n), b(n);
b.FillRand();

// you can ask to solve A^T x = b
// either with Solve(SeldonTrans, x) or with TransSolve
x = b;
mat_lu.TransSolve(x);

// for multiple right hand sides, they are rearranged in a matrix
// each column is a different right hand side
int nb_rhs = 10;
Matrix<double, General, ColMajor> X(n, nb_rhs);
X.FillRand();

mat_lu.TransSolve(X);

Location :

EnableOutOfCore

Syntax :

  void EnableOutOfCore();
  void DisableOutOfCore();

This method allows the direct solver to write a part of the matrix on the disk. This option is enabled only for Mumps.

Location :

SetPivotThreshold

Syntax :

  void SetPivotThreshold(double eps);

This method is not available for each direct solver, it allows to set the threshold used for pivoting.

Example :

MatrixPastix<double> mat_lu;

mat_lu.SetPivotThreshold(1e-8);

// then you can use mat_lu to factorize and solve as usual
mat_lu.FactorizeMatrix(A);
mat_lu.Solve(x);

Location :

Pastix.cxx
Wsmp.cxx

RefineSolution, DoNotRefineSolution

Syntax :

  void RefineSolution();
  void DoNotRefineSolution();

This method is not available for each direct solver, it includes (or not) an additional refinement step in the resolution phase. The obtained solution should be more accurate, but the cost of the solving step should be twice higher at least.

Example :

MatrixPastix<double> mat_lu;

// you tell Pastix that you want a refinement
mat_lu.RefineSolution();

// then you can use mat_lu to factorize and solve as usual
mat_lu.FactorizeMatrix(A);
mat_lu.Solve(x);

Location :

Pastix.cxx
Wsmp.cxx

SetCoefficientEstimationNeededMemory, SetMaximumCoefficientEstimationNeededMemory, SetIncreaseCoefficientEstimationNeededMemory

Syntax :

  void SetCoefficientEstimationNeededMemory(double);
  void SetMaximumCoefficientEstimationNeededMemory(double);
  void SetIncreaseCoefficientEstimationNeededMemory(double)

This method is available for MUMPS only, it sets a safety coefficient to a given value (first method). then this coefficient will be multiplied by a given factor (third method) until the factorization is successful. This loop is stopped is the coefficient reaches a given maximum value (second method).

Example :

MatrixMumps<double> mat_lu;

// you set the initial safety coefficient
// to allocate needed memory for the factorisation
mat_lu.SetCoefficientEstimationNeededMemory(1.2);

// you set the increase
// here it means that if the factorisation fails (because the memory allocated is not sufficient)
// the safety coef will be multiplied by 1.5, and again until completion
mat_lu.SetIncreaseCoefficientEstimationNeededMemory(1.5);

// you set the maximum value allowed
// it means that if the safety coefficient is greater than 100.0, the solver will abort
// we consider that the matrix has a problem
mat_lu.SetMaximumCoefficientEstimationNeededMemory(100.0);

// then you can use mat_lu to factorize and solve as usual
mat_lu.FactorizeMatrix(A);
mat_lu.Solve(x);

Location :

Pastix.cxx
Wsmp.cxx

PerformAnalysis, PerformFactorization

Syntax :

  void PerformAnalysis(Matrix& A);
  void PerformFactorization(Matrix& A);

The method "PerformAnalysis" should reorder matrix, and perform a "symbolic" factorization of the matrix, whereas the method "PerformFactorization" should perform only numerical factorization. The aim here is to reduce the amount of work when we consider a family of linear systems which have the same pattern. In that configuration, the input matrix is not cleared.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a Mumps solver (or SparseDirectSolver)
MatrixMumps<double> mat_lu;
// symbolic factorization
mat_lu.PerformAnalysis(A);
// then factorization
mat_lu.PerformFactorization(A);
// you can solve a system
mat_lu.Solve(x);

// then construct another matrix A, but with the same 
// pattern. For example FillRand will modify the values
A.FillRand();
// and solve the new system
mat_lu.PerformFactorization(A);
mat_lu.Solve(x);

Location :

FindOrdering

Syntax :

  void FindOrdering(Matrix&, Vector<int>&);
  void FindOrdering(Matrix&, Vector<int>&, bool);

This method computes the new row numbers of the matrix by using available algorithms in Mumps/Pastix. It is currently not implemented for UmfPack/SuperLU.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a Mumps structure (it works also for MatrixPastix)
MatrixMumps<double> mat_lu;
IVect permutation;
// we find the best numbering of A 
// by default, the matrix A is erased
mat_lu.FindOrdering(A, permutation);
// if last argument is true, A is not modified
mat_lu.FindOrdering(A, permutation, true);

Location :

Mumps.cxx
Pastix.cxx

FindSparseOrdering

Syntax :

  void FindSparseOrdering(Matrix&, Vector<int>&, int type_ordering);

This function computes a reordering array for a given matrix. The different types of ordering are listed in the method SelectOrdering. Some orderings may be unavailable if Seldon is not interfaced with direct solvers.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
IVect permutation;
// we find the best ordering for A
FindSparseOrdering(A, permutation, SparseMatrixOrdering::SCOTCH);

Location :

Ordering.cxx

GetSchurMatrix (only for MatrixMumps)

Syntax :

  void GetSchurMatrix(Matrix&, Vector<int>&, Matrix&);
  void GetSchurMatrix(Matrix&, Vector<int>&, Matrix&, bool);

This method computes the schur complement when a matrix and row numbers of the Schur matrix are provided. It is equivalent to use the function GetSchurMatrix.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a Mumps structure
MatrixMumps<double> mat_lu;

// then you set some row numbers num that will be associated
// with a sub-matrix A22 : A = [A11 A12; A21 A22]
IVect num(3);
num(0) = 4; num(1) = 10; num(2) = 12;

// computation of Schur complement : A22 - A21 A11^-1 A12
Matrix<double> schur_cplt;
mat_lu.GetSchurMatrix(A, num, schur_cplt);

// the size of matrix schur_cplt should be the same as the size of num

Location :

GetSchurComplement

Syntax :

  void GetSchurComplement(Matrix&, Vector<int>&, Matrix&);

This method computes the schur complement when a matrix and row numbers of the Schur matrix are provided. It is equivalent to use the function GetSchurMatrix.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare a solver
SparseDistributedSolver<double> mat_lu;

// then you set some row numbers num that will be associated
// with a sub-matrix A22 : A = [A11 A12; A21 A22]
IVect num(3);
num(0) = 4; num(1) = 10; num(2) = 12;

// computation of Schur complement : A22 - A21 A11^-1 A12
Matrix<double> schur_cplt;
mat_lu.GetSchurMatrix(A, num, schur_cplt);

// the size of matrix schur_cplt should be the same as the size of num

Location :

Mumps.cxx
UmfPack.cxx
Pastix.cxx
SuperLU.cxx

GetLU

Syntax :

  void GetLU(Matrix&, SparseDirectSolver&, bool);
  void GetLU(Matrix&, MatrixMumps&, bool);
  void GetLU(Matrix&, MatrixUmfPack&, bool);
  void GetLU(Matrix&, MatrixSuperLU&, bool);
  void GetLU(Matrix&, MatrixMumps&, bool);
  void GetLU(Matrix&, MatrixUmfPack&, bool);
  void GetLU(Matrix&, MatrixSuperLU&, bool);

This method performs a LU factorization. The last argument is optional. When omitted, the initial matrix is erased, when equal to true, the initial matrix is not modified.

Example :

// sparse matrices, use of Mumps for example
MatrixMumps<double> mat_lu;
Matrix<double, General, ArrayRowSparse> Asp(n, n);
// you add all non-zero entries to matrix Asp
// then you call GetLU to factorize the matrix
GetLU(Asp, mat_lu);
// Asp is empty after GetLU
// you can solve Asp x = b 
X = B;
SolveLU(mat_lu, X);
// if you want to solve Asp^T x = b
X = B;
SolveLU(SeldonTrans, mat_lu, X);

// you can also solve Asp x = b for multiple right hand sides
// b is a matrix, each right hand side is a column of b
Matrix<double, General, ColMajor> bmat;
SolveLU(SeldonNoTrans, mat_lu, bmat);

// if you want to keep initial matrix
GetLU(Asp, mat_lu, true);

Location :

SolveLU

Syntax :

  void SolveLU(MatrixMumps&, Vector&);
  void SolveLU(const SeldonTranspose, MatrixMumps&, Vector&);
  void SolveLU(MatrixMumps&, Matrix&);
  void SolveLU(const SeldonTranspose, MatrixMumps&, Matrix&);

  void SolveLU(MatrixPastix&, Vector&);
  void SolveLU(const SeldonTranspose, MatrixPastix&, Vector&);
  void SolveLU(MatrixPastix&, Matrix&);
  void SolveLU(const SeldonTranspose, MatrixPastix&, Matrix&);

  void SolveLU(SparseDirectSolver&, Vector&);
  void SolveLU(const SeldonTranspose, SparseDirectSolver&, Vector&);
  void SolveLU(SparseDirectSolver&, Matrix&);
  void SolveLU(const SeldonTranspose, SparseDirectSolver&, Matrix&);

  void SolveLU(MatrixSuperLU&, Vector&);
  void SolveLU(const SeldonTranspose, MatrixSuperLU&, Vector&);
  void SolveLU(MatrixSuperLU&, Matrix&);
  void SolveLU(const SeldonTranspose, MatrixSuperLU&, Matrix&);

  void SolveLU(MatrixUmfPack&, Vector&);
  void SolveLU(const SeldonTranspose, MatrixUmfPack&, Vector&);
  void SolveLU(MatrixUmfPack&, Matrix&);
  void SolveLU(const SeldonTranspose, MatrixUmfPack&, Matrix&);

  void SolveLU(MatrixWsmp&, Vector&);
  void SolveLU(const SeldonTranspose, MatrixWsmp&, Vector&);
  void SolveLU(MatrixWsmp&, Matrix&);
  void SolveLU(const SeldonTranspose, MatrixWsmp&, Matrix&);

  void SolveLU(SparseSeldonSolver&, Vector&);
  void SolveLU(const SeldonTranspose, SparseSeldonSolver&, Vector&);
  void SolveLU(SparseSeldonSolver&, Matrix&);
  void SolveLU(const SeldonTranspose, SparseSeldonSolver&, Matrix&);

This method uses the LU factorization computed by GetLU in order to solve a linear system or its transpose. The right hand side is overwritten by the result. An example is given in the documentation of GetLU.

Location :

Mumps.cxx
UmfPack.cxx
Pastix.cxx
SuperLU.cxx

GetCholesky

Syntax :

  void GetCholesky(Matrix&, MatrixPastix&, bool keep = false);
  void GetCholesky(Matrix&, MatrixCholmod&, bool keep = false);
  void GetCholesky(Matrix&, int print_level = 0);

This method performs a Cholesky factorization of a symmetric positive definite matrix. The last argument is optional. When omitted, the initial matrix is erased, when equal to true, the initial matrix is not modified. The last syntax can be used to perform a Cholesky factorization of a matrix of type ArrayRowSymSparse without pivoting (very slow). The last argument is the level of verbosity (0 = no messages are displayed).

Example :

// sparse matrices, use of Cholmod for example
MatrixCholmod mat_chol;
Matrix<double, Symmetric, ArrayRowSymSparse> Asp(n, n);
// you add all non-zero entries to matrix Asp
// then you call GetCholesky to factorize the matrix
GetCholesky(Asp, mat_lu, true);

// you can solve Asp x = b 
X = B;
SolveCholesky(SeldonNoTrans, mat_chol, X);
SolveCholesky(SeldonTrans, mat_chol, X);

// or compute y = L x
Vector<double> y = X;
MltCholesky(SeldonNoTrans, mat_chol, y);

// if your matrix is already ordered correctly
// and small enough, you can call GetCholesky directly
// but usually it is slow, it is better to use
// the interface with Pastix or Cholmod
GetCholesky(Asp);

// SolveCholesy and MltCholesky are available directly
X = B;
SolveCholesky(SeldonNoTrans, Asp, X);
SolveCholesky(SeldonTrans, Asp, X);

y = X;
MltCholesky(SeldonNoTrans, Asp, y);

// for better efficiency, you can convert Asp to RowSymSparse
Matrix<double, Symmetric, RowSymSparse> C;
Copy(Asp, C);

// SolveCholesky and MltCholesky should be more efficient
X = B;
SolveCholesky(SeldonNoTrans, C, X);
SolveCholesky(SeldonTrans, C, X);

y = X;
MltCholesky(SeldonNoTrans, C, y);

Location :

Cholmod.cxx
Pastix.cxx
SparseCholeskyFactorisation.cxx

SolveCholesky

Syntax :

  void SolveCholesky(SeldonTranspose, MatrixCholmod&, Vector&);
  void SolveCholesky(SeldonTranspose, MatrixPastix&, Vector&);
  void SolveCholesky(SeldonTranspose, Matrix&, Vector&);

This method uses the Cholesky factorization (A = L L^T) computed by GetCholesky in order to solve the system L x = b or L^T x = b. The right hand side b is overwritten by the result x. An example is given in the documentation of GetCholesky.

Location :

Cholmod.cxx
Pastix.cxx
SparseCholeskyFactorisation.cxx

MltCholesky

Syntax :

  void MltCholesky(SeldonTranspose, MatrixCholmod&, Vector&);
  void MltCholesky(SeldonTranspose, MatrixPastix&, Vector&);
  void MltCholesky(SeldonTranspose, Matrix&, Vector&);

This method uses the Cholesky factorization (A = L L^T) computed by GetCholesky in order to compute y = L x or y = L^T x. The input vector x is overwritten by the result y. An example is given in the documentation of GetCholesky.

Location :

Cholmod.cxx
Pastix.cxx
SparseCholeskyFactorisation.cxx

GetSchurMatrix (only for MatrixMumps)

Syntax :

  void GetSchurMatrix(Matrix&, MatrixMumps&, Vector<int>&, Matrix&);

This method computes the so-called Schur matrix (or Schur complement) from a given matrix.

Example :

MatrixMumps<double> mat_lu;
Matrix<double, General, ArrayRowSparse> A(n, n);
// you add all non-zero entries to matrix A
// then you specify row numbers for schur matrix
IVect num(5); 
num.Fill();
// somehow, A can be written under the form A = [A11 A12; A21 A22]
// A22 is related to row numbers of the Schur matrix
// Schur matrix is dense
Matrix<double> mat_schur(5, 5);
GetSchurMatrix(A, mat_lu, num, mat_schur);

// then you should obtain A22 - A21*inv(A11)*A12 -> mat_schur

Location :

FactorizeMatrix (for MatrixCholdmod)

Syntax :

  void FactorizeMatrix(Matrix& A);

This method completes the Cholesky factorization (A = L L^T) by calling Cholmod.

Example :

MatrixCholmod<double> mat_chol;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Cholesky factorization (A = L L^T)
// the L factor is stored in mat_chol
mat_chol.FactorizeMatrix(A);

// then you can solve L y = b
Vector<double> y(n), b(n);
b.FillRand(); y = b;
mat_chol.Solve(SeldonNoTrans, y);

// and L^T x = y to obtain the solution of A x = b
Vector<double> x(y);
mat_chol.Solve(SeldonTrans, x);

Location :

Cholmod.cxx

Solve (for MatrixCholdmod)

Syntax :

  void Solve(const SeldonTranspose&, Vector& x);

This method solves L x = b or L^T x = b by calling Cholmod. The L factor is assumed to have been previously computed (by calling FactorizeMatrix). On input, the user provides the right hand side in vector x, on output x contains the solution.

Example :

MatrixCholmod<double> mat_chol;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Cholesky factorization (A = L L^T)
// the L factor is stored in mat_chol
mat_chol.FactorizeMatrix(A);

// then you can solve L y = b
Vector<double> y(n), b(n);
b.FillRand(); y = b;
mat_chol.Solve(SeldonNoTrans, y);

// and L^T x = y to obtain the solution of A x = b
Vector<double> x(y);
mat_chol.Solve(SeldonTrans, x);

Location :

Cholmod.cxx

Mlt (for MatrixCholdmod)

Syntax :

  void Mlt(const SeldonTranspose&, Vector& x);

This method computes y = L x or y = L^T x by calling Cholmod. The L factor is assumed to have been previously computed (by calling FactorizeMatrix). On input, the user provides the vector x, on output y contains the product with L or its transpose.

Example :

MatrixCholmod<double> mat_chol;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Cholesky factorization (A = L L^T)
// the L factor is stored in mat_chol
mat_chol.FactorizeMatrix(A);

// then you can compute y = L x
Vector<double> y(n), x(n);
x.FillRand(); y = x;
mat_chol.Mlt(SeldonNoTrans, y);

// or  y = L^T x
y = x;
mat_chol.Mlt(SeldonTrans, x);

Location :

Cholmod.cxx

FactorizeMatrix

Syntax :

  void FactorizeMatrix(Matrix&);
  void FactorizeMatrix(Matrix&, bool);

This method factorizes the given matrix. By default, the input matrix is cleared after the factorization. If you want to keep the input_matrix, you have to provide true in the second argument. In parallel execution, this method will consider that the linear system to solve is restricted to the current processor. For example, you can use this method to factorize independant linear systems. If the matrix is distributed overall severall processors, you should call FactorizeDistributed or use the class SparseDistributedSolver.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare the mumps solver
MatrixMumps<double> mat_lu;
// you factorize the matrix
mat_lu.FactorizeMatrix(A);
// to know if there is a problem
int info = mat_lu.GetInfoFactorization();
if (info == mat_lu.OUT_OF_MEMORY)
  {
    cout << "The matrix is too large, not enough memory" << endl;
  }
  
// once the matrix is factorized, you can solve systems
Vector<double> x(n);
mat_lu.Solve(x);

// filling another matrix
A.Reallocate(n, n); A.Get(0, 2) = 1.0; // ...

// if you want to keep the input matrix after the factorization
mat_lu.FactorizeMatrix(A, true);

Location :

Solve

Syntax :

  void Solve(Vector&);
  void Solve(SeldonTranspose, Vector&);
  void Solve(SeldonTranspose, Matrix<T, General, ColMajor>&);

This method computes the solution of A x = b or A^T x = b, assuming that GetLU (or FactorizeMatrix) has been called before. This is equivalent to use function SolveLU.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare the mumps solver
MatrixMumps<double> mat_lu;
// you factorize the matrix
mat_lu.FactorizeMatrix(A);
// once the matrix is factorized, you can solve systems
Vector<double> x(n), b(n);
b.Fill();
x = b;
mat_lu.Solve(x);
// to solve A^T x = b :
x = b;
mat_lu.Solve(SeldonTrans, x);

// you can solve system with multiple right hand sides
Matrix<double, General, ColMajor> B(n, 10);
B.FillRand();
// B is overwritten by the solution
mat_lu.Solve(SeldonNoTrans, B);

// and transpose system can be solved with multiple right hand sides
B.FillRand();
mat_lu.Solve(SeldonTrans, B);

Location :

FactorizeMatrix

Syntax :

  void FactorizeMatrix(Matrix&);
  void FactorizeMatrix(Matrix&, bool);

This method factorizes the given matrix. By default, the input matrix is cleared after the factorization. If you want to keep the input_matrix, you have to provide true in the second argument.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare the mumps solver
MatrixPardiso<double> mat_lu;
// you factorize the matrix
mat_lu.FactorizeMatrix(A);
  
// once the matrix is factorized, you can solve systems
Vector<double> x(n);
mat_lu.Solve(x);

// filling another matrix
A.Reallocate(n, n); A.Get(0, 2) = 1.0; // ...

// if you want to keep the input matrix after the factorization
mat_lu.FactorizeMatrix(A, true);

Location :

Pardiso.cxx

Solve

Syntax :

  void Solve(Vector&);
  void Solve(SeldonTranspose, Vector&);
  void Solve(SeldonTranspose, Matrix<T, General, ColMajor>&);

This method computes the solution of A x = b or A^T x = b, assuming that GetLU (or FactorizeMatrix) has been called before. This is equivalent to use function SolveLU.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A(n, n);
// you declare the mumps solver
MatrixPardiso<double> mat_lu;
// you factorize the matrix
mat_lu.FactorizeMatrix(A);
// once the matrix is factorized, you can solve systems
Vector<double> x(n), b(n);
b.Fill();
x = b;
mat_lu.Solve(x);
// to solve A^T x = b :
x = b;
mat_lu.Solve(SeldonTrans, x);

// you can solve system with multiple right hand sides
Matrix<double, General, ColMajor> B(n, 10);
B.FillRand();
// B is overwritten by the solution
mat_lu.Solve(SeldonNoTrans, B);

// and transpose system can be solved with multiple right hand sides
B.FillRand();
mat_lu.Solve(SeldonTrans, B);

Location :

Pardiso.cxx

FactorizeMatrix (for MatrixPastix)

Syntax :

  void FactorizeMatrix(Matrix& A);

This method completes the Cholesky or LU factorization (A = L L^T or A = LU) by calling Pastix.

Example :

MatrixPastix<double> mat_pastix;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Crout factorization (A = L D L^T)
// for unsymmetric matrices, it is the LU factorization (A = L U)
mat_pastix.FactorizeMatrix(A);

// then you can solve A x = b
Vector<double> x(n), b(n);
b.FillRand(); x = b;
mat_pastix.Solve(SeldonNoTrans, y);

Location :

Solve (for MatrixPastix)

Syntax :

   void Solve(const Vector& x);
  void Solve(const SeldonTranspose&, Vector& x);
  void Solve(const SeldonTranspose&, Matrix& x);

For a LU (or Crout) factorization, this method solves A x = b or A^T x = b by calling Pastix. You can provide multiple right hand sides by providing a matrix. Each column of this matrix is a right hand side. For a Cholesky factorization, this method solves L x = b or L^T x = b by calling Pastix. The L factor is assumed to have been previously computed (by calling FactorizeMatrix). On input, the user provides the right hand side in vector x, on output x contains the solution.

Example :

MatrixPastix<double> mat_pastix;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Crout factorization (A = L D L^T)
mat_pastix.FactorizeMatrix(A);

// then you can solve A x = b
Vector<double> x(n), b(n);
b.FillRand(); x = b;
mat_pastix.Solve(SeldonNoTrans, y);

// you can also solve A^T x = b
x = b;
mat_pastix.Solve(SeldonTrans, y);

// or solve A x = b (or A^T x = b) with multiple right hand sides
int nrhs = 4;
Matrix<double, General, ColMajor> bmat(n, nrhs);
bmat.FillRand();
mat_pastix.Solve(SeldonTrans, bmat);

Location :

Mlt (for MatrixPastix)

Syntax :

  void Mlt(const SeldonTranspose&, Vector& x);

This method computes y = L x or y = L^T x by calling Pastix. The L factor is assumed to have been previously computed (by calling FactorizeMatrix). On input, the user provides the vector x, on output y contains the product with L or its transpose.

Example :

MatrixPastix<double> mat_pastix;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Cholesky factorization (A = L L^T)
// the L factor is stored in mat_pastix
mat_pastix.SetCholeskyFacto(true);
mat_pastix.FactorizeMatrix(A);

// then you can compute y = L x
Vector<double> y(n), x(n);
x.FillRand(); y = x;
mat_pastix.Mlt(SeldonNoTrans, y);

// or  y = L^T x
y = x;
mat_pastix.Mlt(SeldonTrans, x);

Location :

SetCholeskyFacto

Syntax :

  void SetCholeskyFacto(bool);

By default, MatrixPastix computes a LU factorization (or LDL^T factorization for a symmetric matrix). If you want to complete a Cholesky factorization, you need to call this method.

Example :

MatrixPastix<double> mat_pastix;
Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n);
// you can fill A
A.AddInteraction(0, 0, 2.5); // etc

// and then computes the Cholesky factorization (A = L L^T)
// the L factor is stored in mat_pastix
mat_pastix.SetCholeskyFacto(true);
mat_pastix.FactorizeMatrix(A);

// then you can compute y = L x
Vector<double> y(n), x(n);
x.FillRand(); y = x;
mat_pastix.Mlt(SeldonNoTrans, y);

// or  y = L^T x
y = x;
mat_pastix.Mlt(SeldonTrans, x);

Location :