Direct solvers

Seldon is interfaced with libraries performing the direct resolution of very large sparse linear systems : MUMPS, SuperLU and UmfPack. An example file is located in test/program/direct_test.cpp. If you want to test the interface with MUMPS, assuming that MUMPS has been compiled in directory MumpsDir, you can compile this file by typing :

g++ -DSELDON_WITH_MUMPS direct_test.cpp -IMumpsDir/include -IMumpsDir/libseq \
  -IMetisDir/Lib -LMumpsDir/lib -ldmumps -lzmumps -lmumps_common -lpord \
  -LMumpsDir/libseq -lmpiseq -LScotchDir/lib -lesmumps -lfax -lsymbol \
  -ldof -lorder -lgraph -lscotch -lscotcherr -lcommon \
  -lgfortran -lm -lpthread -LMetisDir -lmetis -llapack -lblas

You can simplify the last command, if you didn't install Metis and Scotch and didn't compile MUMPS with those libraries. For UmfPack, if UmfPackDir denotes the directory where UmfPack has been installed, you have to type :

g++ -DSELDON_WITH_UMFPACK direct_test.cpp
  -IUmfPackDir/AMD/Include -IUmfPackDir/UMFPACK/Include \
  -IUmfPackDir/UMFPACK/UFconfig -LUmfPackDir/UMFPACK/Lib \
  -lumfpack -LUmfPackDir/AMD/Lib -lamd -llapack -lblas

For SuperLU, the compilation line reads (SuperLUdir is the directory where SuperLU is located) :

g++ -DSELDON_WITH_SUPERLU direct_test.cpp -ISuperLUdir/SRC -LSuperLUdir -lsuperlu

All in all, MUMPS seems more efficient, and includes more functionnalities than the other libraries.

Syntax

The syntax of all direct solvers is similar 

void GetLU(Matrix&, MatrixMumps&);
void SolveLU(MatrixMumps&, Vector&);

The interface has been done only for double precision (real or complex numbers), since single precision is not accurate enough when very large sparse linear systems are solved.

Basic use

We provide an example of direct resolution using SuperLU.

// first we construct a sparse matrix
int n = 1000; // size of linear system
// we assume that you know how to fill arrays values, ptr, ind
Matrix<double, General, RowSparse> A(n, n, nnz, values, ptr, ind);

// then we declare vectors
Vector<double> b(n), x(n);

// you fill right-hand side
b.Fill();

// you perform the factorization (real matrix)
MatrixSuperLU<double> mat_lu;
GetLU(A, mat_lu);

// then you can solve as many linear systems as you want
x = b;
SolveLU(mat_lu, x);

Methods of MatrixMumps :

Clear Releases memory used by factorization
InitSymmetricMatrix Informs the solver that the matrix is symmetric
InitUnSymmetricMatrix Informs the solver that the matrix is unsymmetric
SelectOrdering selects an ordering scheme
HideMessages requires to the direct solver to not display any message
ShowMessages requires a normal display by the direct solver
GetInfoFactorization returns the error code generated by the factorization
FindOrdering computes a new ordering of rows and columns
FactorizeMatrix performs LU factorization
GetSchurMatrix forms Schur complement
Solve uses LU factorization to solve a linear system

Methods of MatrixUmfPack :

Clear Releases memory used by factorization
HideMessages requires to the direct solver to not display any message
ShowMessages requires a normal display by the direct solver
GetInfoFactorization returns the error code generated by the factorization
FactorizeMatrix performs LU factorization
Solve uses LU factorization to solve a linear system

Methods of MatrixSuperLU :

Clear Releases memory used by factorization
HideMessages requires to the direct solver to not display any message
ShowMessages requires a normal display by the direct solver
GetInfoFactorization returns the error code generated by the factorization
FactorizeMatrix performs LU factorization
Solve uses LU factorization to solve a linear system

Functions :

GetLU performs a LU factorization
SolveLU uses LU factorization to solve a linear system
GetSchurMatrix forms Schur complement

Clear

Syntax :

  void Clear();

This method releases the memory used by the factorisation. It is available for every direct solver. A call to that method is necessary before asking for a new factorization.

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);
// you exploit this factorization
// you need to clear LU matrix before doing another factorization
mat_lu.Clear();
// you fill A with other values
// and refactorize
GetLU(A, mat_lu);

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx

InitSymmetricMatrix, InitUnSymmetricMatrix (only for MatrixMumps)

Syntax :

  void InitSymmetricMatrix();
  void InitUnSymmetricMatrix();

These methods are used before FactorizeMatrix so that Mumps knows that the input matrix is symmetric or not. UmfPack and SuperLU don't have that kind of method since they don't have specific algorithms for symmetric matrices. If you are using GetLU and SolveLU, you don't need to call those methods.

Location :

Mumps.cxx

SelectOrdering (only for MatrixMumps)

Syntax :

  void SelectOrdering(int);

You can force a specific ordering scheme to be used (Metis, Amd, Scotch, etc). In the documentation of Mumps, you can find a list of available orderings and the associated code.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a Mumps structure
MatrixMumps<double> mat_lu;
// you change the default ordering scheme
mat_lu.SelectOrdering(2);
// then you can use getlu as usual
GetLU(A, mat_lu);

Location :

Mumps.cxx

FindOrdering (only for MatrixMumps)

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.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a Mumps structure
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

ShowMessages, HideMessages

Syntax :

  void ShowMessages();
  void HideMessages();

ShowMessages allows the direct solver to display informations about the factorization and resolution phases, while HideMessages forbids any message to be displayed.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a Mumps structure
MatrixMumps<double> mat_lu;
// then before GetLU, you can require silencious mode
mat_lu.HideMessages();
GetLU(A, mat_lu);

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx

GetInfoFactorization

Syntax :

  void GetInfoFactorization();

This method returns the error code provided by the used direct solver.

Example :

// you fill a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
// you declare a Mumps structure
MatrixMumps<double> mat_lu;
// you factorize the matrix
GetLU(A, mat_lu);
// to know if there is a problem
int info = mat_lu.GetInfoFactorization();

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx

FactorizeMatrix

Syntax :

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

This method factorizes the given matrix. It is better to use GetLU, since FactorizeMatrix needs that InitSymmetricMatrix is called before (for Mumps).

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.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 better to use the function GetSchurMatrix since you don't need to call InitSymmetricMatrix before.

Location :

Mumps.cxx

Solve

Syntax :

  void Solve(Vector&);

This method computes the solution of A x = b, assuming that GetLU has been called before to factorize matrix A. This is equivalent to use function SolveLU.

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx

GetLU

Syntax :

  void GetLU(Matrix&, MatrixMumps&);
  void GetLU(Matrix&, MatrixUmfPack&);
  void GetLU(Matrix&, MatrixSuperLU&);
  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 keep initial matrix
GetLU(Asp, mat_lu, true);

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.cxx

SolveLU

Syntax :

  void SolveLU(MatrixMumps&, Vector&);
  void SolveLU(MatrixUmfPack&, Vector&);
  void SolveLU(MatrixSuperLU&, Vector&);

This method uses the LU factorization computed by GetLU in order to solve a linear system. The right hand side is overwritten by the result.

Location :

Mumps.cxx
UmfPack.cxx
SuperLU.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 :

Mumps.cxx

Email Contact