| Installation
|
The syntax of all iterative solvers is the same
int Gmres(const Matrix& A, Vector& x, const Vector& b,
Preconditioner& M, Iteration& iter);
The first argument A is the matrix to solve. If the flag SELDON_WITH_VIRTUAL is enabled, A must derive from the class VirtualMatrix, M from the class Preconditioner_Base. The arguments x and b are template, they can be of type Vector or DistributedVector (for parallel execution). If the flag SELDON_WITH_VIRTUAL is disabled, A and M are of generic type, you only need to define the functions Mlt and MltAdd for A, and the methods Solve and TransSolve for M. We highly recommend to activate the flag SELDON_WITH_VIRTUAL such that you can use the compiled version of Seldon for any type of matrix that derive from VirtualMatrix. If you are including SeldonLib.hxx, this flag is enabled by default. For all iterative solvers, the first argument is the matrix, the second and third argument are vectors, x contains the initial guess on input, the solution on output, b contains the right-hand-side. The fourth argument is the preconditioner. For the moment, there is an implementation of the identity preconditioner (no preconditioning), a SOR (Successive Over Relaxation) preconditioner and incomplete factorization (ILU(k) and ILUT). The last argument is a Seldon structure defining the parameters of the iteration. The input matrix is assumed to be squared for all the iterative solvers.
We provide an example of iterative solution using Seldon structures for the matrices and the vectors. An example file is provided in test/program/iterative_test.cpp, unitary tests are present in test/unit/precond_test.cc. An example file with parallel iterative solver is present in test/unit/distributed_solver.cc.
// 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 and initial guess b.Fill(); x.Zero();; // initialization of iteration parameters int nb_max_iter = 1000; // relative stopping criterion double tolerance = 1e-6; Iteration<double> iter(nb_max_iter, tolerance); // identity preconditioner -> no preconditioning Preconditioner_Base<double> precond; // then you can call an iterative solver, Cg for example Cg(A, x, b, precond, iter); // if you are using Gmres, you can set the restart parameter // by default, this parameter is equal to 10 iter.SetRestart(5); Gmres(A, x, b, precond, iter);
In this section, we assume that the flag SELDON_WITH_VIRTUAL is activated. We will show an example where we construct a new class for a matrix that we don't store. If you want to test iterative solvers with that kind of strategy, you can take a look at the file test/program/iterative_test.cpp. The main thing to do is to construct a matrix class that derive from VirtualMatrix and a preconditioning class that derive from Preconditioner_Base and overload the methods MltVector, MltAddVector (and Solve, TransSolve) that are called by iterative solvers to perform the matrix-vector product (and apply the preconditioning or its transpose).
#include "SeldonLib.hxx"
using namespace Seldon;
// Class for a new type of matrix.
template<class T>
class BlackBoxMatrix : public VirtualMatrix<T>
{
protected :
// internal data used to represent the matrix
int n;
T beta;
Seldon::Vector<T> lambda;
public :
// basic constructor , we call the constructor of VirtualMatrix
// with the number of rows and columns (to initialize m_ and n_ of VirtualMatrix)
BlackBoxMatrix(int n0, const T& beta_)
: VirtualMatrix<T>(n0, n0)
{
beta = beta_; n = n0;
lambda.Reallocate(n);
for (int i = 0; i < n; i++)
lambda(i) = i+1;
}
// this method is used by iterative solvers
// matrix vector product y = A x
void MltVector(const Vector<T>& x, Vector<T>& y)
{
// y = S^{-1} x
y = x;
for (int i = (n-2); i >= 0; i--)
y(i) -= beta*y(i+1);
// y = B y
for (int i = 0; i < n; i++)
y(i) *= lambda(i);
// y = S C
for (int i = 0; i < (n-1); i++)
y(i) += beta*y(i+1);
}
// this method is used by iterative solvers
// matrix vector product y = A^T x
void MltVector(const SeldonTranspose& trans, const Vector<T>& x, Vector<T>& y)
{
if (trans.NoTrans())
{
MltVector(x, y);
return;
}
// Transpose of S B S^{-1} is S^{-t} B S^t
y = x;
// Y = S^t Y
for (int i = (n-1); i >= 1; i--)
y(i) += beta*y(i-1);
// Y = B Y
for (int i = 0; i < n; i++)
y(i) *= lambda(i);
// Y = S^{-t} Y
for (int i = 1; i < n; i++)
y(i) -= beta*y(i-1);
}
// this method is used by iterative solvers
// matrix vector product y = gamma y + alpha A x
void MltAddVector(const T& alpha, const Vector<T>& x, const T& gamma, Vector<T>& y)
{
Vector<T> z;
MltVector(x, z);
y = gamma*y + alpha*z;
}
// this method is used by iterative solvers
// matrix vector product y = gamma y + alpha A^T x
void MltAddVector(const T& alpha, const SeldonTranspose& trans, const Vector<T>& x, const T& gamma, Vector<T>& y)
{
Vector<T> z;
MltVector(trans, x, z);
y = gamma*y + alpha*z;
}
};
// class for preconditioning
template<class T>
class MyPreconditioner : public Preconditioner_Base<T>
{
private:
// example where the preconditioning is of same type as matrix but
// with different beta
BlackBoxMatrix<T> B;
public :
MyPreconditioner(int n, double beta) : B(n, beta) {}
// apply the preconditioning, i.e. solving M r = z
void Solve(const VirtualMatrix<T>& mat, const Vector<T>& r, Vector<T>& z)
{
// mat is the matrix you want to solve (mat x = b)
B.MltVector(r, z);
}
// solving transpose(M) r = z
void TransSolve(const VirtualMatrix<T>& mat, const Vector<T>& r, Vector<T>& z)
{
B.MltVector(SeldonTrans, r, z);
}
};
int main()
{
// now it is very classical like the previous example
int n = 20; double beta = 0.5;
BlackBoxMatrix<double> A(n, beta);
Vector<double> b(n), x(n);
// you fill right-hand side and initial guess
b.Fill();
x.Fill(0);
// initialization of iteration parameters
int nb_max_iter = 1000;
double tolerance = 1e-6;
Iteration<double> iter(nb_max_iter, tolerance);
// your own preconditioner
MyPreconditioner<double> precond(n, 1.2);
// then you can call an iterative solver, Qmr for example
Qmr(A, x, b, precond, iter);
return 0;
}
By default, three preconditioning are provided : Identity (Preconditioner_Base), SOR (SorPreconditioner) and incomplete factorization (IlutPreconditioning). An interface with Hypre is provided, an example file is located in test/program/hypre_test.cc and can be compiled as follows:
mpicxx -o run -DSELDON_WITH_HYPRE test/program/hypre_test.cc -I. -IHypreDir/src/hypre/include -DSELDON_WITH_MPI -LHypreDir/src/hypre/lib -lHYPRE
where HypreDir is the directory where Hypre is located.
| Solve | Applies the preconditioner |
| TransSolve | Applies the transpose of the preconditioner |
| InitSymmetricPreconditioning | Symmetric SOR will be used |
| InitUnSymmetricPreconditioning | SOR will be used |
| SetParameterRelaxation | changes the relaxation parameter |
| SetNumberIterations | changes the number of iterations |
| Solve | Applies the preconditioner |
| TransSolve | Applies the transpose of the preconditioner |
| Clear | clears memory used by incomplete factorisation |
| ShowMessages | enables messages that will be displayed during factorization |
| HideMessages | disables messages that will be displayed during factorization |
| GetFactorisationType | returns the type of incomplete factorisation |
| GetFillLevel | returns the fill-level k (if ILU(k) is set) |
| GetAdditionalFillNumber | returns the number of additional elements per row (used for ILUT(k)) |
| GetPrintLevel | returns the verbose level |
| GetPivotBlockInteger | returns the maximum k = |i-j| allowed when pivoting |
| SetFactorisationType | sets the type of incomplete factorisation |
| SetFillLevel | sets the fill-level k (if ILU(k) is set) |
| SetAdditionalFillNumber | sets the number of additional elements per row (used for ILUT(k)) |
| SetPrintLevel | sets the verbose level |
| SetPivotBlockInteger | sets the maximum k = |i-j| allowed when pivoting |
| SetSymmetricAlgorithm | sets symmetric incomplete factorisation |
| SetUnsymmetricAlgorithm | sets unsymmetric incomplete factorisation (even for symmetric matrices) |
| GetDroppingThreshold | returns threshold to determine elements to drop |
| SetDroppingThreshold | sets threshold to determine elements to drop |
| GetDiagonalCoefficient | returns diagonal coefficient used in ILUD |
| SetDiagonalCoefficient | sets diagonal coefficient used in ILUD |
| GetPivotThreshold | returns threshold used when pivoting columns |
| SetPivotThreshold | sets threshold used when pivoting columns |
| FactorizeMatrix | performs incomplete factorisation |
| Solve | Applies the preconditioner |
| TransSolve | Applies the transpose of the preconditioner |
| SetPreconditioner | Sets preconditioning used in Hypre |
| SetSmoother | Sets smoothing used in Hypre |
| SetLevelEuclid | Sets parameter k for Euclid solver in Hypre |
| ShowMessages | Enables the display of statistics in Hypre |
| SetInputPreconditioning | Sets the preconditioning parameters in Hypre |
| ConstructPreconditioner | Constructs the preconditioning by calling Hypre |
| Solve | applies the preconditioning by calling Hypre |
| TransSolve | applies the transpose preconditioning by calling Hypre |
| Constructors | |
| GetRestart | returns restart parameter |
| SetRestart | changes restart parameter |
| GetTolerance | returns stopping criterion |
| SetTolerance | changes stopping criterion |
| SetMaxIterationNumber | changes maximum number of iterations |
| GetNumberIteration | returns iteration number |
| SetNumberIteration | changes iteration number |
| ShowMessages | displays residual each 100 iterations |
| HideMessages | displays nothing |
| ShowFullHistory | displays residual each iteration |
| SaveFullHistory | saves the history of residuals in a file |
| Init | provides right hand side |
| First | returns true for the first iteration |
| IsInitGuess_Null | returns true if the initial guess is zero |
| SetInitGuess | informs if the initial guess is zero or not |
| Finished | returns true if the stopping criteria are satisfied |
| Fail | informs that the iterative solver failed |
| ErrorCode | returns error code |
| SOR | Performs SOR iterations |
| BiCg | BIConjugate Gradient |
| BiCgcr | BIConjugate Gradient Conjugate Residual |
| BiCgStab | BIConjugate Gradient STABilized |
| BiCgStabl | BIConjugate Gradient STABilized (L) |
| Cg | Conjugate Gradient |
| Cgne | Conjugate Gradient Normal Equation |
| Cgs | Conjugate Gradient Squared |
| CoCg | Conjugate Orthogonal Conjugate Gradient |
| Gcr | Generalized Conjugate Residual |
| Gmres | Generalized Minimum RESidual |
| Lsqr | Least SQuaRes |
| MinRes | Minimum RESidual |
| QCgs | Quasi Conjugate Gradient Squared |
| Qmr | Quasi Minimum Residual |
| QmrSym | Quasi Minimum Residual SYMmetric |
| Symmlq | SYMMetric Least sQuares |
| TfQmr | Transpose Free Quasi Minimum Residual |
void Solve(const VirtualMatrix&, const Vector&, Vector&); void TransSolve(const VirtualMatrix&, const Vector&, Vector&);
These methods should be overloaded if you want to define your own preconditioner since they define the application of the preconditioner or its transpose to a vector.
// constructor of a matrix int n = 20; Matrix<double> A(n, n); A.Fill(); // declaration of a preconditioner Preconditioner_Base<double> M; // vectors Vector<double> r(n), z(n); r.Fill(); // now we apply preconditioning, i.e. solving M z = r // for Preconditioner_Base, it will set z = r (identity preconditioner) M.Solve(A, r, z); // we can also apply the transpose of preconditioner // i.e. solving transpose(M) z = r M.TransSolve(A, r, z);
Class Preconditioner_Base
Iterative.hxx
Iterative.cxx
void Solve(const Matrix&, const Vector&, Vector&); void TransSolve(const Matrix&, const Vector&, Vector&);
These methods define the application of the preconditioner and its transpose to a vector. The used preconditioner is SOR (Successive Over Relaxation). It can be used in its symmetric version (SSOR), or the unsymmetric version. In this last case, the application of the preconditioner consists of a forward sweep while the transpose consists of a backward sweep. If you are using Seldon structures of sparse matrices, the function SOR is available. If you are using other structures, it is necessary to overload the function SOR (if SELDON_WITH_VIRTUAL is not activated). If SELDON_WITH_VIRTUAL is activated, you need to overload the method ApplySor derived from the class VirtualMatrix.
// constructor of a matrix int n = 20; Matrix<double, General, RowSparse> A(n, n); A.Fill(); // declaration of a preconditioner SorPreconditioner<double> M; // by default, relaxation parameter omega = 1 // number of iterations = 1 // you can change that M.SetParameterRelaxation(1.5); M.SetNumberIterations(2); // if you prefer to use symmetric version M.InitSymmetricPreconditioning(); // vectors Vector<double> r(n), z(n); r.Fill(); // now we apply preconditioning, i.e. solving M z = r // for Preconditioner_Base, it will set z = r (identity preconditioner) M.Solve(A, r, z); // we can also apply the transpose of preconditioner // i.e. solving transpose(M) z = r M.TransSolve(A, r, z);
Class SorPreconditioner
Precond_Ssor.cxx
void InitSymmetricPreconditioning(); void InitUnsymmetricPreconditioning();
InitSymmetricPreconditioning sets SSOR as preconditioning. The symmetric SOR consists of a forward sweep followed by a backward sweep. InitUnSymmetricPreconditioning sets SOR, it consists of a forward sweep for the preconditioner, and a backward sweep for the transpose of the preconditioner.
// declaration of a preconditioner SorPreconditioner<double> M; // by default, symmetric preconditioning // use InitUnSymmetricPreconditioning to force a non-symmetric preconditioning M.InitUnSymmetricPreconditioning();
Class SorPreconditioner
Precond_Ssor.cxx
void SetParameterRelaxation(const T& omega);
This methods changes the relaxation parameter.
// declaration of a preconditioner SorPreconditioner<double> M; // relaxation parameter omega M.SetParameterRelaxation(1.6);
Class SorPreconditioner
Precond_Ssor.cxx
void SetNumberIterations(int);
This methods changes the number of SOR iterations.
// declaration of a preconditioner SorPreconditioner<double> M; // relaxation parameter omega M.SetParameterRelaxation(1.6); // number of SOR iterations each time the preconditioning is applied M.SetNumberIterations(1); // 1 is actually the default choice
Class SorPreconditioner
Precond_Ssor.cxx
void Clear()
This methods clears the factorisation stored in the structure.
Class IlutPreconditioning
IlutPreconditioning.cxx
int GetFactorisationType() const void SetFactorisationType(int type)
These methods return (or set) the type of incomplete factorisation set. The following types of incomplete factorisation are available :
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_K): // for ILU_K, you need to set the fill level k ilut.SetFillLevel(2); // we want a symmetric preconditioning ilut.SetSymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Cg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
int GetFillLevel() const void SetFillLevel(int level);
These methods return (and set) the level k associated with ILU(k) factorisation.
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_K): // for ILU_K, you need to set the fill level k ilut.SetFillLevel(2); // we want a symmetric preconditioning ilut.SetSymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Cg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
int GetAdditionalFillNumber() const void SetAdditionalFillNumber(int k)
These methods return (and set) the number of fill-in elements allowed per row during the factorisation.
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILUT): // for ILUT, you need to set the allowed additional fill-in elements // for each row. If you put n, there will be no limit. ilut.SetAdditionalFillNumber(n); ilut.SetDroppingThreshold(0.01); // we want a symmetric preconditioning ilut.SetSymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Cg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
int GetPrintLevel() const void SetPrintLevel(int print_level)
These methods return (and set) the print level.
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_K): // for ILU_K, you need to set the fill level k ilut.SetFillLevel(2); // if you want messages to be displayed during factorisation ilut.SetPrintLevel(2); // we want a symmetric preconditioning ilut.SetSymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Cg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
int GetPivotBlockInteger() const void SetPivotBlockInteger(int k)
These methods return (and set) the maximal difference between diagonal index and pivot index allowed when pivoting.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILUT): // pivot are searching within the block |i-j| <= k // if you allow any pivot, put n ilut.SetPivotBlockInteger(n); ilut.SetPivotThreshold(0.1); ilut.SetAdditionalFillNumber(n); ilut.SetDroppingThreshold(0.01); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
double GetDroppingThreshold() const void SetDroppingThreshold(double eps)
These methods return (and set) the dropping threshold.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILUT): // pivot are searching within the block |i-j| <= k // if you allow any pivot, put n ilut.SetPivotBlockInteger(n); ilut.SetPivotThreshold(0.1); ilut.SetAdditionalFillNumber(n); // threshold used to drop small values of the factorisation ilut.SetDroppingThreshold(0.01); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
double GetPivotThreshold() const void SetPivotThreshold(double eps)
These methods return (and set) the pivot threshold.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILUT): // pivot are searching within the block |i-j| <= k // if you allow any pivot, put n ilut.SetPivotBlockInteger(n); // threshold used to determine if pivoting is needed ilut.SetPivotThreshold(0.1); ilut.SetAdditionalFillNumber(n); // threshold used to drop small values of the factorisation ilut.SetDroppingThreshold(0.01); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
double GetDiagonalCoefficient() const void SetDiagonalCoefficient(double eps)
These methods return (and set) the coefficient used in the compensation of diagonal.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_D): // pivot are searching within the block |i-j| <= k // if you allow any pivot, put n ilut.SetPivotBlockInteger(n); // threshold used to determine if pivoting is needed ilut.SetPivotThreshold(0.1); ilut.SetAdditionalFillNumber(n); // threshold used to drop small values of the factorisation ilut.SetDroppingThreshold(0.01); // diagonal compensation coefficient for ILU_D ilut.SetDiagonalCoefficient(0.9); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
void FactorizeMatrix(const Vector<int>& permut, Matrix&); void FactorizeMatrix(const Vector<int>& permut, Matrix&, bool keep_matrix);
This method performs the incomplete factorisation of the given matrix. By default, the matrix is cleared during the process.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_D): // pivot are searching within the block |i-j| <= k // if you allow any pivot, put n ilut.SetPivotBlockInteger(n); // threshold used to determine if pivoting is needed ilut.SetPivotThreshold(0.1); ilut.SetAdditionalFillNumber(n); // threshold used to drop small values of the factorisation ilut.SetDroppingThreshold(0.01); // diagonal compensation coefficient for ILU_D ilut.SetDiagonalCoefficient(0.9); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); // put true if you want to keep the original matrix ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
void Solve(Vector&) void TransSolve(Vector&) void Solve(const Matrix&, const Vector&, Vector&) void TransSolve(const Matrix&, const Vector&, Vector&)
This method performs the approximate resolution of A x = b or AT x = b by using the incomplete factorisation.
// constructing a non-symmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // type of incomplete factorisation ilut.SetFactorisationType(ilut.ILU_K): ilut.SetFillLevel(2); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); // put true if you want to keep the original matrix ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system (Solve and TransSolve will be called by Qmr) Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Qmr(A, x, b, ilut, iter); // if you want to call Solve independently x = b; ilut.Solve(x);
Class IlutPreconditioning
IlutPreconditioning.cxx
void SetSymmetricAlgorithm() const
This methods informs to construct a symmetric preconditioning (if the given matrix is symmetric).
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // we want a symmetric preconditioning ilut.SetSymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); Cg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
void SetUnsymmetricAlgorithm() const
This methods informs to construct an unsymmetric preconditioning (even if the given matrix is symmetric).
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // initialization of the preconditioning IlutPreconditioning<double> ilut; // we want a non-symmetric preconditioning ilut.SetUnsymmetricAlgorithm(); // then you can compute the incomplete factorisation IVect permutation(n); permutation.Fill(); ilut.FactorizeMatrix(permutation, A, true); // and solve the linear system // with an iterative algorithm adpated to non-symmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter(1000, 1e-6); BiCg(A, x, b, ilut, iter);
Class IlutPreconditioning
IlutPreconditioning.cxx
Iteration(); Iteration(int, const T&);
The second constructor specifies the maximum number of iterations and the stopping criterion
// constructing a symmetric matrix int n = 50; Matrix<double, Symmetric, ArrayRowSymSparse> A(n, n); // identity preconditioning Preconditioner_Base<double> prec; // and solve the linear system // with an iterative algorithm adpated to symmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); // here, we impose 1000 as the maximal number of iterations // and 1e-6 as the stopping criterion // as soon as the number of iterations is larger than 1000 // or the relative residual lower than 1e-6, the iterative algorithm is stopped Iteration<double> iter(1000, 1e-6); Cg(A, x, b, prec, iter);
int GetRestart(); void SetRestart(int);
These methods give access to the restart parameter, which is used by some iterative solvers, e.g. BiCgSTAB(l), Gmres and Gcr. The default value is equal to 10.
// constructing an unsymmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // identity preconditioning Preconditioner_Base<double> prec; // and solve the linear system // with an iterative algorithm adpated to unsymmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter; // stopping criterion iter.SetTolerance(1e-8); // maximal number of iterations iter.SetMaxNumberIteration(10000); // restart parameter iter.SetRestart(20); // you can retrieve the restart parameter you gave int m = iter.GetRestart(); // the iterative solver is called Gmres(A, x, b, prec, iter);
Class Iteration
Iterative.hxx
Iterative.cxx
T GetTolerance(); void SetTolerance(T);
These methods give access to the stopping criterion.
// constructing an unsymmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // identity preconditioning Preconditioner_Base<double> prec; // and solve the linear system // with an iterative algorithm adpated to unsymmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter; // stopping criterion iter.SetTolerance(1e-8); // maximal number of iterations iter.SetMaxNumberIteration(10000); // restart parameter iter.SetRestart(20); // you can retrieve your stopping criterion double eps = iter.GetTolerance(); // the iterative solver is called Gmres(A, x, b, prec, iter);
Class Iteration
Iterative.hxx
Iterative.cxx
void SetMaxIterationNumber(int);
This method gives access to the maximum iteration number.
Class Iteration
Iterative.hxx
Iterative.cxx
int GetNumberIteration(); void SetNumberIteration(int);
This method gives access to the iteration number.
// constructing an unsymmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // identity preconditioning Preconditioner_Base<double> prec; // and solve the linear system // with an iterative algorithm adpated to unsymmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter; // stopping criterion iter.SetTolerance(1e-8); // maximal number of iterations iter.SetMaxNumberIteration(10000); // restart parameter iter.SetRestart(20); // the iterative solver is called Gmres(A, x, b, prec, iter); // if you want to know how many iterations have been needed to complete the solution int nb_iter = iter.GetNumberIteration(); // SetNumberIteration is not really useful here // since each iterative solver resets this number to 0 iter.SetNumberIteration(0);
Class Iteration
Iterative.hxx
Iterative.cxx
void ShowMessages(); void HideMessages(); void ShowFullHistory();
If ShowMessages is called, it will display residual each 100 iterations during the resolution. If HideMessages is called, there is no display at all, whereas ShowFullHistory displays residual at each iteration.
// constructing an unsymmetric matrix int n = 50; Matrix<double, General, ArrayRowSparse> A(n, n); // identity preconditioning Preconditioner_Base<double> prec; // and solve the linear system // with an iterative algorithm adpated to unsymmetric linear systems Vector<double> x(n), b(n); b.FillRand(); x.Zero(); Iteration<double> iter; // stopping criterion iter.SetTolerance(1e-8); // maximal number of iterations iter.SetMaxNumberIteration(10000); // restart parameter iter.SetRestart(20); // you can ask to display the residual at each iteration iter.ShowFullHistory(); // the iterative solver is called Gmres(A, x, b, prec, iter); // by default, no messages are displayed // to see residuals each 100 iterations : iter.ShowMessages()
Class Iteration
Iterative.hxx
Iterative.cxx
void SaveFullHistory(const string& file);
If this method is called, the history of residuals will be written in the file name given in this method.
int nb_max_iter = 100; double tol = 1e-6;
Iteration<double> iter(nb_max_iter, tol);
iter.SaveFullHistory("residue.dat");
// we consider that A, x, b and M have been constructed previously
// the residuals for Cg will be written in the file residue.dat
Cg(A, x, b, M, iter);Class Iteration
Iterative.hxx
Iterative.cxx
void Init(const Vector&);
This method is used by iterative solvers to initialize the computation of residuals. Since relative residual is computed, we need to know the norm of the first residual. This method should not be called in a regular use.
Class Iteration
Iterative.hxx
Iterative.cxx
bool First();
This method returns true for the first iteration. This method should not be called in a regular use.
Class Iteration
Iterative.hxx
Iterative.cxx
bool IsInitGuess_Null(); void SetInitGuess(bool);
SetInitGuess allows you to inform the iterative solver that your initial guess is null. This can spare one matrix-vector product, which is good if the expected number of iterations is small.
int n = 10; Matrix<double, General, ArrayRowSparse> A(n, n); Vector<double> x(n), b(n); A.Fill(); b.Fill(); Iteration<double> iter(100, 1e-6); // you inform that initial guess is null iter.SetInitGuess(true); // if the initial guess is null // and x is non-null, the solver enforces x to be 0 Preconditioner_Base<double> M; Gmres(A, x, b, M, iter);
Class Iteration
Iterative.hxx
Iterative.cxx
bool Finished(const Vector&);
This method is used by iterative solvers to know if the stopping criterion is reached. This method also displays the residual if required. This method should not be called in a regular use.
Class Iteration
Iterative.hxx
Iterative.cxx
bool Fail(int, const string&);
This method is used by iterative solvers when a breakdown occurs, often due to a division by 0. This method should not be called in a regular use.
Class Iteration
Iterative.hxx
Iterative.cxx
int ErrorCode();
This method returns the error code after the solution. If the solution process has been successful, it should return 0. This method should not be called in a regular use.
Class Iteration
Iterative.hxx
Iterative.cxx
void SOR(const Matrix& A, Vector& x, const Vector& b, const T& omega, int nb_iter); void SOR(const Matrix& A, Vector& x, const Vector& b, const T& omega, int nb_iter, int type); void SOR(SeldonTrans, const Matrix& A, Vector& x, const Vector& b, const T& omega, int nb_iter); void SOR(SeldonTrans, const Matrix& A, Vector& x, const Vector& b, const T& omega, int nb_iter, int type);
This method tries to solve A x = b with n iterations of SOR. If SeldonTrans is provided as a first argument, it will solve AT x = b. nb_iter is the number of iterations. type is optional and describes the stage of relaxation to perform. If equal to 2, a forward relaxation is performed, if equal to 3 a backward relaxation is performed, if equal to 0 a forward relaxation is performed followed by a backward relaxation (giving a symmetric algorithm called SSOR).
int n = 1000; Matrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b // 2 Sor iterations (forward sweep) with omega = 0.5 double omega = 0.5; int nb_iterations = 2; SOR(A, x, b, omega, nb_iterations); // you can ask for symmetric SOR: forward sweeps followed // by backward sweeps SOR(A, x, b, omega, nb_iterations, 0); // if you need backward sweep SOR(A, x, b, omega, nb_iterations, 3);
int BiCg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICG algorithm. This algorithm can solve complex general linear systems and calls matrix-vector products with the transpose matrix.
int BiCgcr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGCR algorithm. This algorithm can solve symmetric complex linear systems.
int BiCgStab(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int BiCgStabl(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB(l) algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int Cg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CG algorithm. This algorithm can solve real symmetric or hermitian linear systems.
int Cgne(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CGNE algorithm. This algorithm can solve complex general linear systems and calls matrix vector products with the transpose matrix.
int Cgs(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CGS algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int Cocg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm. This algorithm can solve complex symmetric linear systems.
int Gcr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int Gmres(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using restarted GMRES algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int Lsqr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using LSQR algorithm. This algorithm can solve complex general linear systems and calls matrix vector products with the transpose matrix.
int MinRes(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm. This algorithm can solve complex symmetric linear systems and doesn't call matrix vector products with the transpose matrix.
int QCgs(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using QCGS algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
int Qmr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using QMR algorithm. This algorithm can solve complex general linear systems and calls matrix vector products with the transpose matrix.
int QmrSym(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using symmetric QMR algorithm. This algorithm can solve complex symmetric linear systems.
int Symmlq(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using SYMMLQ algorithm. This algorithm can solve complex symmetric linear systems.
int TfQmr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using TFQMR algorithm. This algorithm can solve complex general linear systems and doesn't call matrix vector products with the transpose matrix.
void SetPreconditioner(int type);
This method selects which preconditioning must be used when calling Hypre. The following choices are available
You can look at the documentation of Hypre to know details about each preconditioning.</>
int n = 1000; DistributedMatrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b with an iterative solver // first the preconditioning is constructed HyprePreconditioner<double> prec; prec.SetPreconditioner(prec.BOOMER_AMG); prec.SetSmoother(prec.GS_SEQ); prec.ConstructPreconditioner(A, true); // then the solution is computed Iteration<double> iter(1000, 1e-8); BiCg(A, x, b, prec, iter);
void SetSmoother(int type);
This sets the smoother to be used in Hypre. The following smoothers are available
You can look at the documentation of Hypre to know details about each smoother.
int n = 1000; DistributedMatrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b with an iterative solver // first the preconditioning is constructed HyprePreconditioner<double> prec; prec.SetPreconditioner(prec.BOOMER_AMG); prec.SetSmoother(prec.GS_SEQ); prec.ConstructPreconditioner(A, true); // then the solution is computed Iteration<double> iter(1000, 1e-8); BiCg(A, x, b, prec, iter);
void SetLevelEuclid(int lvl);
This method sets the level k for Euclid solver implemented in Hypre.
int n = 1000; DistributedMatrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b with an iterative solver // first the preconditioning is constructed HyprePreconditioner<double> prec; prec.SetPreconditioner(prec.EUCLID); prec.SetLevelEuclid(2); prec.ConstructPreconditioner(A, true); // then the solution is computed Iteration<double> iter(1000, 1e-8); BiCg(A, x, b, prec, iter);
void ShowMessages();
This method will enable the display of messages of Hypre.
int n = 1000; DistributedMatrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b with an iterative solver // first the preconditioning is constructed HyprePreconditioner<double> prec; prec.SetPreconditioner(prec.BOOMER_AMG); prec.SetSmoother(prec.GS_SEQ); prec.ShowMessages(); prec.ConstructPreconditioner(A, true); // then the solution is computed Iteration<double> iter(1000, 1e-8); BiCg(A, x, b, prec, iter);
void SetInputPreconditioning(string, Vector<string>)
This method sets the solver to use. The first parameter is the name of the solver, the second parameter contains additional parameters. You can look at the file Hypre.cxx to see the list of parameters.
int n = 1000;
DistributedMatrix<double, General, RowSparse> A;
// you initialize A as you want (SetData for example)
// then vectors
Vector<double> x(n), b(n);
x.Zero();
b.Fill();
// we want to solve A x = b with an iterative solver
// first the preconditioning is constructed
HyprePreconditioner<double> prec;
Vector<string> param;
param.PushBack("MaximumLevel"); param.PushBack("4");
param.PushBack("NumberSweeps"); param.PushBack("2");
param.PushBack("Smoother"); param.PushBack("GaussSeidel");
prec.SetInputPreconditioning("AMG", param);
prec.ShowMessages();
prec.ConstructPreconditioner(A, true);
// then the solution is computed
Iteration<double> iter(1000, 1e-8);
BiCg(A, x, b, prec, iter);
void ConstructPreconditioner(DistributedMatrix A, bool keep_matrix);
This method constructs the preconditioning with a given matrix. If the second argument is false, the matrix is cleared during the construction.
int n = 1000;
DistributedMatrix<double, General, RowSparse> A;
// you initialize A as you want (SetData for example)
// then vectors
Vector<double> x(n), b(n);
x.Zero();
b.Fill();
// we want to solve A x = b with an iterative solver
// first the preconditioning is constructed
HyprePreconditioner<double> prec;
Vector<string> param;
param.PushBack("MaximumLevel"); param.PushBack("4");
param.PushBack("NumberSweeps"); param.PushBack("2");
param.PushBack("Smoother"); param.PushBack("GaussSeidel");
prec.SetInputPreconditioning("AMG", param);
prec.ShowMessages();
// we put true to keep the matrix A
// but if you construct a preconditioning from another matrix B
// you can clear this matrix by putting false for the second argument
prec.ConstructPreconditioner(A, true);
// then the solution is computed
Iteration<double> iter(1000, 1e-8);
BiCg(A, x, b, prec, iter);void Solve(SeldonTranspose, A, r, z) void Solve(A, r, z) void TransSolve(A, r, z)
The method Solve applies the preconditioning z = M-1 r while the method TransSolve applies the transpose preconditioning. Usually you do not have to call these methods since they will be called automatically by the iterative solver.
int n = 1000;
DistributedMatrix<double, General, RowSparse> A;
// you initialize A as you want (SetData for example)
// then vectors
Vector<double> r(n), z(n);
z.Zero();
r.Fill();
// we want to solve A x = b with an iterative solver
// first the preconditioning is constructed
HyprePreconditioner<double> prec;
Vector<string> param;
param.PushBack("MaximumLevel"); param.PushBack("4");
param.PushBack("NumberSweeps"); param.PushBack("2");
param.PushBack("Smoother"); param.PushBack("GaussSeidel");
prec.SetInputPreconditioning("AMG", param);
prec.ShowMessages();
prec.ConstructPreconditioner(A, true);
// you can apply the preconditioning for a given vector
// z = M^{-1} r
prec.Solve(A, r, z);