CODEWITHSUNDEEP
c++sparse-matrixlinear-algebraarmadillo
user1512321
user1512321

Reputation:

Using GMRES to solve a complex sparse matrix with armadillo

I am porting a MATLAB program to C++. Matlab is able to solve a linear system, (complex values, sparse matrix) that my C++ library, armadillo, cannot.

Does armadillo have the ability to use an iterative solver such as GMRES? Unfortunately armadillo's documentation is quite spartan with details on the solver.

Otherwise, can I call LAPACK from armadillo? Or what alternative C++ library provides a strong linear solver?

(The matrix is quite ill-conditioned and I'm also looking for other solutions to this problem, including pre-conditioning matricies)

Upvotes: 1

Views: 1472

Answers (1)

francis
francis

Reputation: 9817

The sparse matrix of Armadillo features a function dedicated to solving linear system, spsolve(). As pointed out by the documentation, it uses SuperLU whenever it is installed and enabled in the configuration file config.hpp of Armadillo. Otherwise, it converts the matrix to dense format and use LAPACK. LAPACK handles various matrix formats, mostly dense and band format but not a generic sparse format like compresed sparse row.

SuperLU implements the LU factorization of a sparse matrix. It makes use of additionnal permutation matrices to obtain sparse L and U factors. But a lot of alternatives exist, implementing various strategies, such as MUMPS or UMFPACK (multifrontal LU factorization) or PaStiX or ... The page of Suitesparse tells that Matlab "backslash" makes use of UMFPACK or CHOLMOD., likely behing LU and Cholevsky factorisation in the diagram.

Furthermore, the Eigen library wraps many solvers. As a consequence, the solver can be chosen depending on the matrix properties and user choice.

Take a look at the Petsc library: it wraps a lot of preconditionners and solvers, including GMRES. See this example to understand how a matrix can be populated and solved. The user can choose its solver and preconditionner at runtime! In the linked exemple, the preconditionner is set to ICC (Incomplete Cholesky), but the solver (KSP) is set from option. This web page helps specifying a starting set of options. For gmres:

 -ksp_type gmres -ksp_monitor -ksp_monitor_true_residual -ksp_rtol <rtol> -ksp_norm_type unpreconditioned

The monitor will tell you how fast the residual decreases. rtol specifies the relative decrease in the (possibly preconditioned) residual norm. The norm type unpreconditioned specifies that the norm of the unpreconditionned residual must be computed to test the convergence. Since your matrix is ill-conditionned, monitoring the true unpreconditionned residual may be a good idea. Otherwise, if the preconditionning matrix happens to be ill-conditionned, the preconditionned problem may be solved accurately, but the residual of the initial problem may still be huge.

See this page for some options of the GMRES solver (number of Krylov directions... ) Once the program is compiled, try a run with -help, to list available options.

Upvotes: 2

Related Questions