Iterative solvers¶
Most of the solvers provided in Gmm++ come frorm ITL with slight modifications (gmres has been optimized and adapted for complex matrices). Include the file gmm/gmm_iter_solvers.h to use them.
iterations¶
The iteration object of Gmm++ is a modification of the one in ITL. This is not a template type as in ITL.
The simplest initialization is:
gmm::iteration iter(2.0E-10);
where 2.0E-10 is the (relative) residual to be obtained to have the convergence.
Some possibilities:
iter.set_noisy(n) // n = 0 : no output
// n = 1 : output of iterations on the standard output
// n = 2 : output of iterations and sub-iterations
// on the standard output
// ...
iter.get_iteration() // after a computation, gives the number of
// iterations made.
iter.converged() // true if the method converged.
iter.set_maxiter(n) // Set the maximum of iterations.
// A solver stops if the maximum of iteration is
// reached, iter.converged() is then false.
Linear solvers¶
Here is the list of available linear solvers:
gmm::row_matrix< std::vector<double> > A(10, 10); // The matrix
std::vector<double> B(10); // Right hand side
std::vector<double> X(10); // Unknown
gmm::identity_matrix PS; // Optional scalar product for cg
gmm::identity_matrix PR; // Optional preconditioner
...
gmm::iteration iter(10E-9);// Iteration object with the max residu
size_t restart = 50; // restart parameter for GMRES
gmm::cg(A, X, B, PS, PR, iter); // Conjugate gradient
gmm::bicgstab(A, X, B, PR, iter); // BICGSTAB BiConjugate Gradient Stabilized
gmm::gmres(A, X, B, PR, restart, iter) // GMRES generalized minimum residual
gmm::qmr(A, X, B, PR, iter) // Quasi-Minimal Residual method.
gmm::least_squares_cg(A, X, B, iter) // unpreconditionned least square CG.
The solver gmm::constrained_cg(A, C, X, B, PS, PR, iter); solve a system with linear constraints, C is a matrix which represents the constraints. But it is still experimental.
(Version 1.7) The solver gmm::bfgs(F, GRAD, X, restart, iter) is a BFGS quasi-Newton algorithm with a Wolfe line search for large scale problems. It minimizes the function F without constraints, be given its gradient GRAD. restart is the max number of stored update vectors.
Preconditioners¶
The following preconditioners, to be used with linear solvers, are available:
gmm::identity_matrix P; // No preconditioner
gmm::diagonal_precond<matrix_type> P(SM); // diagonal preconditioner
gmm::mr_approx_inverse_precond<matrix_type> P(SM, 10, 10E-17);
// preconditioner based on MR
// iterations
gmm::ildlt_precond<matrix_type> P(SM); // incomplete (level 0) ldlt
// preconditioner. Fast to be
// computed but less efficient than
// gmm::ildltt_precond.
// incomplete ldlt with k fill-in and threshold preconditioner.
// Efficient but could be costly.
gmm::ildltt_precond<matrix_type> P(SM, k, threshold);
gmm::ilu_precond<matrix_type> P(SM); // incomplete (level 0) ilu
// preconditioner. Very fast to be
// computed but less efficient than
// gmm::ilut_precond.
// incomplete LU with k fill-in and threshold preconditioner.
// Efficient but could be costly.
gmm::ilut_precond<matrix_type> P(SM, k, threshold);
// incomplete LU with k fill-in, threshold and column pivoting preconditioner.
// Try it when ilut encounter too small pivots.
gmm::ilutp_precond<matrix_type> P(SM, k, threshold);
Except ildltt\_precond, all these precontionners come from ITL. ilut_precond has been optimized and simplified and cholesky_precond has been corrected and transformed in an incomplete LDLT preconditioner for stability reasons (similarly, we add choleskyt_precond which is in fact an incomplete LDLT with threshold preconditioner). Of course, ildlt\_precond and ildltt_precond are designed for symmetric real or hermitian complex matrices to be use principally with cg.
Additive Schwarz method¶
The additive Schwarz method is a decomposition domain method allowing the resolution of huge linear systems (see [SCHADD] for the principle of the method).
For the moment, the method is not parallelized (this should be done …). The call is the following:
gmm::sequential_additive_schwarz(A, u, f, P, vB, iter, local_solver, global_solver)
A is the matrix of the linear system. u is the unknown vector. f is the right hand side. P is an eventual preconditioner for the local solver. vB is a vector of rectangular sparse matrices (of type const std::vector<vBMatrix>, where vBMatrix is a sparse matrix type), each of these matrices is of size \(N \times N_i\) where \(N\) is the size of A and \(N_i\) the number of variables in the \(i^{th}\) sub-domain ; each column of the matrix is a base vector of the sub-space representing the \(i^{th}\) sub-domain. iter is an iteration object. local_solver has to be chosen in the list gmm::using_gmres(), gmm::using_bicgstab(), gmm::using_cg(), gmm::using_qmr() and gmm::using_superlu() if SuperLu is installed. global_solver has to be chosen in the list gmm::using_gmres(), gmm::using_bicgstab(), gmm::using_cg(), gmm::using_qmr().
The test program schwarz_additive.C is the directory tests of GetFEM is an example of the resolution with the additive Schwarz method of an elastostatic problem with the use of coarse mesh to make a better preconditioning (i.e. one of the sub-domains represents in fact a coarser mesh).
In the case of multiple solves with the same linear system, it is possible to store the preconditioners or the LU factorizations to save computation time.
A (too) simple program in gmm/gmm_domain_decomp.h allows to build a regular domain decomposition with a certain ratio of overlap. It directly produces the vector of matrices vB for the additive Schwarz method.
Range basis function¶
The function gmm\_range\_basis(B, columns, EPS=1e-12) defined in gmm/gmm\_range\_basis.h allows to select from the columns of a sparse matrix B a basis of the range of this matrix. The result is returned in columns which should be of type std::set<size_type> and which contains the indices of the selected columns.
The algorithm is specially designed to select independent constraints from a large matrix with linearly dependent columns.
There is four step in the implemented algorithm
- Elimination of null columns.
- Selection of a set of already orthogonal columns.
- Elimination of locally dependent columns by a blockwise Gram-Schmidt algorithm.
- Computation of vectors of the remaining null space by a global restarted Lanczos algorithm and deduction of some columns to be eliminated.
The algorithm is efficient if after the local Gram-Schmidt algorithm it remains a low dimension null space. The implemented restarted Lanczos algorithm find the null space vectors one by one.
The Global restarted Lanczos algorithm may be improved or replaced by a block Lanczos method (see [ca-re-so1994] for instance), a block Wiedelann method (in order to be parallelized) or simply the computation of more than one vector of the null space at each iteration.
