Solver

*Solver, <solver name>, [<extrapolation>, <iterative solver options>]

This command is used to specify the solver to be used to solve the system of equations of the current step. As default, the Pardiso solver is used. The following solvers can be selected:

<solver name>

Use this keyword to specify the strategy for solving the system of equations. numgeo offers direct solvers as well as iterative solvers.
- Direct solvers: PARDISO, MUMPS, Simple solver, Diagonal solver
- Iterative solvers:
  - LIS library: GMRES, FGMRES, BiCGStab, ...
  - In-house: BiCGStab, BiCGStab(l)
- External solvers: Use this option to interface external solvers to numgeo.
In addition you have the option to deactive the solver.

Direct solvers

PARDISO

PARDISO is a thread-safe and memory efficient software for solving large sparse symmetric and unsymmetric linear systems of equations on shared-memory and distributed-memory multiprocessors ¹. We interface the oneMKL PARDISO implementation. PARDISO is the default solver (no need to enable it using *Solver, Pardiso)

MUMPS

MUMPS - the MUltifrontal Massively Parallel sparse direct Solver ²³. numgeo ships with a shared-memory (OpenMP enabled, MPI disabled). For an MPI-parallelised version see the External solver feature. To enable MUMPS in a step, use *Solver, Mumps.

Simple solver

A very basic (simple) implementation of a direct solver using serial forward substitution. Should only be used for very small systems, e.g. single-element tests. Compared to MUMPS or PARDISO, the overhead is little but efficiency for medium to large systems is drastically worse. To enable the simple solver in a step use *Solver, simple.

Diagonal solver

A simple implementation of a direct solver for strictly diagonal matrices such as used in explicit dynamics with mass lumping. To enable the diagonal solver in a step, use *Solver, diagnoal.

Iterative solvers:

Most numgeo calculations feature highly-nonlinear materials with only approximated Jacobians, frictional contacts with opening and closing contact surfaces, coupled two-phase or three-phase flow simulations. The arising linearised systems of equations are thus notoriously difficult to solve. So far we have not figured out good combinations of iterative solvers and preconditioners that offer robust and reliable performance. We thus recommend using direct solvers. Any help in providing reliable and efficient iterative solvers is welcome - please get in touch with us!

LIS solver library

LIS (Library of Iterative Solvers for linear systems, Pujii et al. (2005)⁴, Katekemori et al. (2005)⁵, Katekemori et al. (2008)⁶, Nishida (2010)⁷) is a parallel software library for solving discretized linear equations and eigenvalue problems that arise in the numerical solution of partial differential equations using iterative methods.

The keyword to activate the LIS solver reads:

*Solver, LIS = <iterative solver options>

The sub-keyword <iterative solver options> is used to define the iterative solver and the preconditioner. Unlike the usual numgeo convention, this command requires the definition of a string following the LIS-convention. Be careful, blanks do matter! Available solvers and preconditioners as well as their LIS-options are depicted in Figure 5.3 and Figure 5.4. The following example displays how to define a GMRES solver with the standard restart value (40). As a preconditioner, an incomplete LU factorisation with a fill in level of 2:

*Solver, LIS=$-i gmres -p ilu -ilu fill 2$

The LIS library gives access to more than 20 different iterative solvers and 10 preconditioners. A summary is given in the following tables:

Solvers:

Solver	Option	Auxiliary Options
CG	`-i {cg\|1}`
BiCG	`-i {bicg\|2}`
CGS	`-i {cgs\|3}`
BiCGSTAB	`-i {bicgstab\|4}`
BiCGSTAB(l)	`-i {bicgstabl\|5}`	`-ell [2]` — The degree l
GPBiCG	`-i {gpbicg\|6}`
TFQMR	`-i {tfqmr\|7}`
Orthomin(m)	`-i {orthomin\|8}`	`-restart [40]` — The restart value m
GMRES(m)	`-i {gmres\|9}`	`-restart [40]` — The restart value m
Jacobi	`-i {jacobi\|10}`
Gauss-Seidel	`-i {gs\|11}`
SOR	`-i {sor\|12}`	`-omega [1.9]` — The relaxation coefficient ω (0 < ω < 2)
BiCGSafe	`-i {bicgsafe\|13}`
CR	`-i {cr\|14}`
BiCR	`-i {bicr\|15}`
CRS	`-i {crs\|16}`
BiCRSTAB	`-i {bicrstab\|17}`
GPBiCR	`-i {gpbicr\|18}`
BiCRSafe	`-i {bicrsafe\|19}`
FGMRES(m)	`-i {fgmres\|20}`	`-restart [40]` — The restart value m
IDR(s)	`-i {idrs\|21}`	`-irestart [2]` — The restart value s
IDR(1)	`-i {idr1\|22}`
MINRES	`-i {minres\|23}`
COCG	`-i {cocg\|24}`
COCR	`-i {cocr\|25}`

Preconditioners

Preconditioner	Option	Auxiliary Options
None	`-p {none\|0}`
Jacobi	`-p {jacobi\|1}`
ILU(k)	`-p {ilu\|2}`	`-ilu_fill [0]` — The fill level k
SSOR	`-p {ssor\|3}`	`-ssor_omega [1.0]` — The relaxation coefficient ω (0 < ω < 2)
Hybrid	`-p {hybrid\|4}`	`-hybrid_i [sor]` — The linear solver `-hybrid_maxiter [25]` — The maximum number of iterations `-hybrid_tol [1.0e-3]` — The convergence tolerance `-hybrid_omega [1.5]` — The relaxation coefficient ω of the SOR (0 < ω < 2) `-hybrid_ell [2]` — The degree l of the BiCGSTAB(l) `-hybrid_restart [40]` — The restart values of the GMRES and Orthomin
I+S	`-p {is\|5}`	`-is_alpha [1.0]` — The parameter α of I + αS(m) `-is_m [3]` — The parameter m of I + αS(m)
SAINV	`-p {sainv\|6}`	`-sainv_drop [0.05]` — The drop criterion
SA-AMG	`-p {saamg\|7}`	`-saamg_unsym [false]` — Select the unsymmetric version (matrix must be symmetric) `-saamg_theta [0.05\\|0.12]` — The drop criterion ( a_{ij}^2 \leq \theta^2
Crout ILU	`-p {iluc\|8}`	`-iluc_drop [0.05]` — The drop criterion `-iluc_rate [5.0]` — The ratio of the maximum fill-in
ILUT	`-p {ilut\|9}`
Additive Schwarz	`-adds true`	`-adds_iter [1]` — The number of iterations

Additional settings:

Option	Description
`-maxiter [1000]`	The maximum number of iterations
`-tol [1.0e-12]`	The convergence tolerance tol
`-tol_w [1.0]`	The convergence tolerance tol_w
`-print [0]`	The output of the residual history: `-print {none\\|0}` — None `-print {mem\\|1}` — Save the residual history `-print {out\\|2}` — Output to standard output `-print {all\\|3}` — Save the residual history and output it to standard output
`-scale [0]`	The scaling (result overwrites the original matrix and vectors): `-scale {none\\|0}` — No scaling `-scale {jacobi\\|1}` — Jacobi scaling $D^{-1}Ax = D^{-1}b$, where $D$ is the diagonal of $A$ `-scale {symm_diag\\|2}` — Diagonal scaling $D^{-1/2}AD^{-1/2}x = D^{-1/2}b$, where $D^{-1/2}$ has $1/\sqrt{a_{ii}}$ on the diagonal
`-initx_zeros [1]`	The behaviour of the initial vector $x_0$: `-initx_zeros {false\\|0}` — Components given by argument `x` in `lis_solve()` `-initx_zeros {true\\|1}` — All components set to 0
`-conv_cond [0]`	The convergence condition: `-conv_cond {nrm2_r\\|0}` — $\\|b - Ax\\|_2 \le tol \cdot \\|b - Ax_0\\|_2$ `-conv_cond {nrm2_b\\|1}` — $\\|b - Ax\\|_2 \le tol \cdot \\|b\\|_2$ `-conv_cond {nrm1_b\\|2}` — $\\|b - Ax\\|_1 \le tol_w \cdot \\|b\\|_1 + tol$
`-omp_num_threads t`	Number of threads (t represents the maximum number of threads)
`-storage [0]`	The matrix storage format
`-storage_block [2]`	The block size of the BSR and BSC formats
`-f [0]`	Precision of the linear solver: `-f {double\\|0}` — Double precision `-f {quad\\|1}` — Double-double (quadruple) precision

Inhouse iterative solvers

numgeo gives access to two "in-house" implementations of BiCGStab solvers:

BiCGStab
BiCGStab(l)

BiCGStab

The BiConjugate Gradient Stabilized method, often abbreviated as BiCGStab, is an iterative method proposed by H. A. van der Vorst ⁸ for the numerical solution of nonsymmetric linear systems. A (right-)preconditioned as well as a version without preconditioning is implemented in numgeo. In cases where BiCGStab does not work well, it may be advantageous to use the BiCGStab(l) method. Example:

*Solver, BiCGStab (using default settings)
*Solver, BiCGStab, pc=none, tol=1e-4, check_rel_solution (using custom settings)

Optional sub-keywords can be used to override the default settings. All sub-keywords are separated by a comma: keyword_1, ..., keyword_n. The available keywords are:

tol = real. Define the convergence tolerance of the iterative method, this is, when the approximated solution is said to be converged. In general, given the linear problem:

\[\boldsymbol{A} \boldsymbol{x} = \boldsymbol{b}\]

with $\boldsymbol{A}$ is the global left-hand-side (lhs) of the system of equations and $\boldsymbol{b}$ the corresponding right-hand-side (rhs), $\boldsymbol{x}$ is the desired solution. For a given initial guess $\boldsymbol{x}_0$, the approximate solution $\boldsymbol{x}_k$ is refined until $\boldsymbol{r}_k = \boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k$ satisfies a given criterion. Per default numgeo checks:

$\frac{||\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k||_2}{||\boldsymbol{b}||_2} \leq \texttt{tol}$

for the BiCGStab(l) algorithm. In the BiCGStab algorithm, the residuum during the iteration process $\boldsymbol{r}_k$ is updated independent of the solution $\boldsymbol{x}_k$ (see the theory manual!), per default numgeo checks

\[\frac{||\boldsymbol{r}_k||_2}{||\boldsymbol{b}||_2} \leq \texttt{tol}\]

which corresponds to the option: check_rel_residuum. However, for the BiCGStab algorithm, there are different convergence checks available for the iterative methods: check_rel_solution, check_solution, check_rel_residuum, check_residuum.

By default, this tolerance is tol = 10$^{-3}$ for general nonlinear procedures. The default tolerance may seem loose for general nonlinear procedures, however it is important to note that the linear solver convergence tolerance is independent of the nonlinear convergence process (i.e., Newton-Raphson method) tolerances that are used to determine if analysis increments converge.
- check_residuum: solution is accepted, when $||\boldsymbol{r}_k||_2 \leq \texttt{tol}$ is satisfied. In this case, you want to choose a rather small value for the tolerance tol, e.g. tol=$10^{-9}$.
- check_rel_solution: solution is accepted, when $\frac{||\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k||_2}{||\boldsymbol{b}||_2} \leq \texttt{tol}$ is satisfied. Compared to the *_residuum controls, this requires two additional AXPY-operations and two additional vector-matrix-products.
- check_solution: solution is accepted, when $||\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k||_2 \leq \texttt{tol}$ is satisfied. In this case, you want to choose a rather small value for the tolerance tol, e.g. tol=$10^{-9}$. Compared to the *_residuum controls, this requires two additional AXPY-operations and two add. vector-matrix-products.
pc = .... Define the preconditioner to be used:
- pc = none: no preconditioner is to be used. The iterative method solves $A b = x$. This is the default preconditioner for the BiCGStab(l) algorithm.
- pc = diagonal: use the inverse diagonal $P_{i,i} = A_{i,i}$ as a preconditioner. The iterative method solves $P^{-1}A b = P^{-1}x$ with $(P_{i,i})^{-1} = 1/A_{i,i}$, but $P_{i,i} = 1$ if $|A_{i,i}| = 0$.
- pc = ls-diagonal: use the least square diagonal preconditioner $P$. $P$ is derived, such that it scales the columns of $A$ to give each unit L2-norm. The iterative method solves $P^{-1}A b = P^{-1}x$. This is the default preconditioner for the BiCGStab algorithm.
- pc = ilu0: Simple implementation of Incomplete LU (ILU) factorization for 0 levels of fill-in

BiCGStab(l)

The BiConjugate Gradient Stabilized(l) method, often abbreviated as BiCGStab(l), is an iterative method proposed by H. A. van der Vorst and D. Fokkema for the numerical solution of nonsymmetric linear systems. The version implemented in numgeo is based on the proposed implementation by D. Fokkema ⁹¹⁰¹¹, which is based on the power basis variant of BiCGstab(l). This variant is enhanced with a more stable way of determination of the iteration coefficients and with a more reliable update strategy for the residuals. Faster convergence in terms of iteration counts may be expected for increasing values of the parameter l. However, since more work is required to obtain the iterate as l increases, optimal performance in terms of computational work may be realized for quite a small value of l. Starting with the value l=2 is recommended. Example:

*Solver, BiCGStab(l) (using default settings)
*Solver, BiCGStab(l), l=4, pc=ls-diagonal, tol=1e-4 (using custom settings)

Optional sub-keywords can be used to override the default settings. All sub-keywords are separated by a comma: keyword_1, ..., keyword_n. The available keywords are:

tol = real. Define the convergence tolerance of the iterative method, this is, when the approximated solution is said to be converged. In general, given the linear problem:

\[\boldsymbol{A} \boldsymbol{x} = \boldsymbol{b}\]

with $\boldsymbol{A}$ is the global left-hand-side (lhs) of the system of equations and $\boldsymbol{b}$ the corresponding right-hand-side (rhs), $\boldsymbol{x}$ is the desired solution. For a given initial guess $\boldsymbol{x}_0$, the approximate solution $\boldsymbol{x}_k$ is refined until $\boldsymbol{r}_k = \boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k$ satisfies a given criterion. Per default numgeo checks:

$\frac{||\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_k||_2}{||\boldsymbol{b}||_2} \leq \texttt{tol}$

for the BiCGStab(l) algorithm. In the BiCGStab algorithm, the residuum during the iteration process $\boldsymbol{r}_k$ is updated independent of the solution $\boldsymbol{x}_k$ (see the theory manual!), per default numgeo checks

\[\frac{||\boldsymbol{r}_k||_2}{||\boldsymbol{b}||_2} \leq \texttt{tol}\]

which corresponds to the option: check_rel_residuum. However, for the BiCGStab algorithm, there are different convergence checks available for the iterative methods: check_rel_solution, check_solution, check_rel_residuum, check_residuum.

By default, this tolerance is tol = 10$^{-3}$ for general nonlinear procedures. The default tolerance may seem loose for general nonlinear procedures, however it is important to note that the linear solver convergence tolerance is independent of the nonlinear convergence process (i.e., Newton-Raphson method) tolerances that are used to determine if analysis increments converge.
pc = .... Define the preconditioner to be used:
- pc = none: no preconditioner is to be used. The iterative method solves $A b = x$. This is the default preconditioner for the BiCGStab(l) algorithm.
- pc = diagonal: use the inverse diagonal $P_{i,i} = A_{i,i}$ as a preconditioner. The iterative method solves $P^{-1}A b = P^{-1}x$ with $(P_{i,i})^{-1} = 1/A_{i,i}$, but $P_{i,i} = 1$ if $|A_{i,i}| = 0$.
- pc = ls-diagonal: use the least square diagonal preconditioner $P$. $P$ is derived, such that it scales the columns of $A$ to give each unit L2-norm. The iterative method solves $P^{-1}A b = P^{-1}x$. This is the default preconditioner for the BiCGStab algorithm.
- pc = ilu0: Simple implementation of Incomplete LU (ILU) factorization for 0 levels of fill-in
l = integer. Define the polynomial degree of the BiCGStab(l) algorithm. By default, the minimal applicable value l=2 is used. Currently supported are l={2,3,4,5,6,7,8}.

External solver

Use this option to interface external solvers to numgeo. For the external solver, an additional executable is required. This executable has to be placed in the same directory as the numgeo executable.

*Solver, external, exe=$executable name$ [, mpi=<n processors>, <extrapolation>]

The (sub)keywords are as follows:

executable: tells numgeo to use an external executable to solve the linear system.
exe=$executable name$: specify the name of the external executable.
[mpi=<n processors>]: specify the number of processors to be used. Notice that:
- for mpi=1, the use of mpi is disabled and the external executable is launched using ./executable name. This is the default.
- for mpi=n (with n $>$ 1), the executable is launched using mpirun -n executable name
[omp=<n threads>]: specify the number of (omp) threads to be used. Notice that:
- for [omp=1], the use of omp is disabled. This is the default.
- for [omp=n] (with n $>$ 1), the executable is launched using export OMP_NUM_THREADS=n

Existing interfaces to external solvers

Interfaces to some external solvers have already been developed and can be freely downloaded from https://github.com/j-machacek/numgeo-external-solvers:

numgeo-mumps-mpi: MPI version of the MUMPS solver
numgeo-ilupack
numgeo-pastix

Example:
*Solver, external, exe=$numgeo-mumps-mpi$, mpi=4, omp=2

Deactivate solver

Use this option to omit using a solver. However, be aware that in this case the system of equation is not solved at all, i.e. no improvement to the initial guess will be calculated (see Theory Manual). This option should best be used in combination with maxiter=0 (see *Step definition). Be aware of the consequences! To deactivate the solving stage in a step of your analysis, use *Solver, None

Arne De Coninck, Bernard De Baets, Drosos Kourounis, Fabio Verbosio, Olaf Schenk, Steven Maenhout, and Jan Fostier. Needles: toward large-scale genomic prediction with marker-by-environment interaction. Genetics, 203(1):543–555, 3 2016. URL: http://dx.doi.org/10.1534/genetics.115.179887, arXiv:http://www.genetics.org/content/203/1/543.full.pdf, doi:10.1534/genetics.115.179887. ↩
Patrick R. Amestoy, Iain S. Duff, Jean-Yves L'Excellent, and Jacko Koster. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications, 23(1):15–41, 1 2001. doi:10.1137/s0895479899358194. ↩
Patrick R. Amestoy, Abdou Guermouche, Jean-Yves L'Excellent, and Stéphane Pralet. Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 32(2):136–156, 2 2006. doi:10.1016/j.parco.2005.07.004. ↩
Akihiro Pujii, Akira Nishida, and Yoshio Oyanagi. Evaluation of Parallel Aggregate Creation Orders: Smoothed Aggregation Algebraic Multigrid Method. In Michael K. Ng, Andrei Doncescu, Laurence T. Yang, and Tau Leng, editors, High Performance Computational Science and Engineering, volume 172, pages 99–122. Springer-Verlag, 2005. doi:10.1007/0-387-24049-7_6. ↩
H. Kotakemori, H. Hasegawa, and A. Nishida. Performance evaluation of a parallel iterative method library using OpenMP. In Eighth International Conference on High-Performance Computing in Asia-Pacific Region (HPCASIA'05), 5 pp.–436. IEEE, 2005. doi:10.1109/HPCASIA.2005.74. ↩
Hisashi Kotakemori, Hidehiko Hasegawa, Tamito Kajiyama, Akira Nukada, Reiji Suda, and Akira Nishida. Performance Evaluation of Parallel Sparse Matrix–Vector Products on SGI Altix3700. In Matthias S. Mueller, Barbara M. Chapman, Bronis R. De Supinski, Allen D. Malony, and Michael Voss, editors, OpenMP Shared Memory Parallel Programming, volume 4315, pages 153–163. Springer Berlin Heidelberg, 2008. doi:10.1007/978-3-540-68555-5_13. ↩
Akira Nishida. Experience in Developing an Open Source Scalable Software Infrastructure in Japan. In David Taniar, Osvaldo Gervasi, Beniamino Murgante, Eric Pardede, and Bernady O. Apduhan, editors, Computational Science and Its Applications – ICCSA 2010, volume 6017, pages 448–462. Springer Berlin Heidelberg, 2010. doi:10.1007/978-3-642-12165-4_36. ↩
H. A. van der Vorst. Bi-CGSTAB: a fast and smoothly converging variant of bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 13(2):631–644, 1992-03. Publisher: Society for Industrial & Applied Mathematics (SIAM). doi:10.1137/0913035. ↩
Gerard LG Sleijpen and Diederik R Fokkema. BiCGstab (ell) for linear equations involving unsymmetric matrices with complex spectrum. Electronic Transactions on Numerical Analysis., 1:11–32, 1993. Publisher: Kent State University. ↩
Gerard LG Sleijpen, Henk A Van der Vorst, and Diederik R Fokkema. BiCGstab (l) and other hybrid bi-cg methods. Numerical Algorithms, 7(1):75–109, 1994. Publisher: Springer. ↩
Diederik R. Fokkema. Enhanced implementation of BiCGstab(l) for solving linear systems of equations. Citeseer, 1996. ↩