7.2 Configuring Linear Solvers¶

FEBio offers several direct and iterative linear solvers that can be used to solve the linear system of equations of the FE model. Users can configure these linear solvers in the configuration file with the default_linear_solver tag. The particular linear solver is selected via the type attribute. The following table lists the possible values, the solver type, and the supported matrix format.

solver	description	type	matrix format
pardiso	A sparse direct solver, from the MKL library. (default)	D	S/U
skyline	A sparse direct solver, but not as efficient as pardiso.	D	S
mkl_dss	The DSS solver from the MKL library.	D	S/U
fgmres	An implementation of the GMRES algorithm, from MKL.	It	U
cg	An implementation of the CG algorithm, from MKL.	It	S
boomeramg	An algebraic multigrid solver, from the HYPRE library.	It	S/U
schur	A block based iterative solver.	It	B
accelerate	Apple-specific sparse solver	D/It	S/U

D=direct solver, It=iterative solver, S = symmetric format, U = unsymmetric format, B = block format.

Iterative solvers can be configured with preconditioners, which attempt to reduce the condition number of the matrix. Preconditioning is highly recommended and often necessary for reducing the number of iterations needed by the iterative solver to converge the solution. Any of the linear solvers listed above can be used as a preconditioner. In addition, the following specialized preconditioners are also available.

preconditioner	description	matrix format
ilu0	Incomplete LU decomposition with no fill-in	unsymmetric
ilut	Incomplete LU decomposition with custom fill-in	unsymmetric
ichol	Incompletly Cholesky decomposition with no fill-in	symmetric
diagonal	Diagonal preconditioner made from the diagonal of the matrix	any

The following sections describe each of these solvers and preconditioners in more detail.

Pardiso¶

The Pardiso solver is an efficient sparse direct linear solver and is the default linear solver. FEBio uses the implementation from the MKL library. It does not require any configuration parameters. It can take symmetric, unsymmetric, and structurally symmetric sparse matrix formats.

Example:

<default_linear_solver type="pardiso"/>

Skyline¶

The Skyline solver was initially added for users who build their own copy of FEBio and did not have access to third party linear solvers. It is a sparse direct solver, but not as efficient as Pardiso and not recommended for large problems. When using this solver, it is highly recommended to set the optimize_bw flag in the model file's Control section, which attempts to minimize the matrix bandwidth. It does not require any configuration parameters and only works with symmetric matrices.

MKL_DSS¶

The MKL library also offers another direct solver, called DSS (Direct-Sparse-Solver). This solver is very promising and appears to perform better for many problems than the pardiso solver. Currently, in FEBio4, pardiso remains the default mostly for backward compatibility reasons, but it is often worthwhile to see if the mkl_dss solver performs better.

FGMRES¶

This iterative linear solver, from the MKL library, implements the GMRES algorithm, which is an efficient Krylov-based iterative solution strategy. It requires several configuration parameters, listed below, and only works with unsymmetric matrices. In the table below, N refers to the number of equations.

parameter	description	default
max_iter	The maximum number of iterations	min(N,150)
tol	relative residual tolerance	1e-6
abs_tol	absolute residual tolerance	1e-12
print_level	The amount of information printed (0,1, or 2)	0
fail_max_iters	Indicates whether reaching max iters is a failure or not (1)	1 (true)
pc_left	The left preconditioner	(none)
pc_right	The right preconditioner	(none)

Comments:

When using FGMRES as a preconditioner it is recommended to set the fail_max_iters flag to zero.

Example 1. This example shows how to set up the FGMRES solver with an ILU0 preconditioner.

<default_linear_solver type="fgmres"><br>  <max_iter>100</max_iter>
  <tol>1e-5</tol>
  <pc_left type="ilu0"/>
</default_linear_solver>

Example 2. This example sets up an FGMRES solver with a Schur preconditioner. For the Schur system we use pardiso for solving the A-block outside of the Schur complement and approximate the A-block inside the Schur complement with its diagonal. This solver requires a block-structured matrix. In order to generate the block structure you need to set the equation_scheme control parameter to 1 in the model input file.

<default_linear_solver type="fgmres"><br> <print_level>2</print_level><br> <max_iter>100</max_iter><br> <tol>1e-5</tol><br> <abs_tol>1e-9</abs_tol><br> <pc_left type="schur"><br>   <print_level>0</print_level><br>   <do_jacobi>0</do_jacobi><br>   <A_solver type="pardiso"/><br>   <schur_pc>0</schur_pc><br>   <schur_A_solver type="diagonal"/><br>   <schur_solver type="fgmres"><br>     <print_level>0</print_level><br>     <max_iter>100</max_iter><br>     <tol>0.05</tol><br>     <fail_max_iters>0</fail_max_iters><br>   </schur_solver><br>  </pc_left><br> </default_linear_solver>

CG¶

This iterative linear solver, from the MKL library, implements the Conjugate Gradient algorithm, which is an efficient Krylov-based iterative solution strategy for symmetric systems. It requires several configuration parameters, listed below, and only works with symmetric matrices. In the table below, N refers to the number of equations.

parameter	description	default
max_iter	The maximum number of iterations	min(N,150)
tol	relative residual tolerance	1e-6
print_level	The amount of information printed (0,1, or 2)	0
pc_left	The left preconditioner	(none)

BoomerAMG¶

This solver, from the HYPRE library (), is an adaptive multigrid solver. It often works better as a preconditioner to an iterative solver than as a stand-alone solver.

parameter	description	default
max_iter	The maximum number of iterations	20
tol	relative residual tolerance	1e-7
print_level	The amount of information printed [0,1, or 2]	0
max_levels	The max number of multigrid levels	25
coarsen_type	Method of grid coarsening (1)	(library default)
use_num_funcs	Use dimensionality of solution variable in coarsening	0 (false)
relax_type	relaxation type (2)	3
interp_type	interpolation type (3)	6
strong_threshold	strong threshold	0.5
p_max_elmts	P max elements	4
num_sweeps	Number of sweeps	1
agg_num_levels	Aggressive number of levels	0
nodal	Nodal option	0 (false)
do_jacobi	Do jacobi preconditioning	0 (false)
fail_max_iters	Fail on max iterations (4)	1 (true)

Comments:

_coarsen_type _selects the algorithm for the coarsening phase. The following values can be used.

value	description
0	CLJP-coarsening (a parallel coarsening algorithm using independent sets)
1	classical Ruge-Stueben coarsening, no boundary treatment (not recommended!)
3	classical Ruge-Stueben coarsening, with boundary treatment
6	Falgout coarsening
7	CLJP-coarsening (using fix random vector, for debugging purposes only)
8	PMIS-coarsening (might lead to slower convergence)
9	PMIS-coarsening (using a fixed random vector ,for debugging purposes only)
10	HMIS-coarsening
11	one-pass Ruge-Stueben coarsening (not recommended!)
21	CGC coarsening
22	CGC-E coarsening

_relax_type _selects the algorithm for the relaxation phase. The following values can be used.

value	description
0	Jacobi
1	Gauss-Seidel, sequential (very slow!)
2	Gauss-Seidel, interior points in parallel, boundary sequential (slow!)
3	hybrid Gauss-Seidel or SOR, forward solve
4	hybrid Gauss-Seidel or SOR, backward solve
5	hybrid chaotic Gauss-Seidel (works only with OpenMP)
6	hybrid symmetric Gauss-Seidel or SSOR
8	l1-scaled hybrid symmetric Gauss-Seidel
9	Gaussian elimination (only on coarsest level)
15	CG (warning - not a fixed smoother - may require RGMRES)
16	Chebyshev
17	FCF-Jacobi
18	l1-scaled Jacobi

_interp_type _selects the algorithm for the interpolation phase. The following values can be used.

value	description
0	classicial modified interpolation
1	LS interpolation (for use with GSMG)
2	classical modified interpolation for hyperbolic PDEs
3	direct interpolation (with separation of weights)
4	multipass interpolation
5	multipass interpolation (with separation of weights)
6	extended+i interpolation
7	extended+i (if no common C neighbor) interpolation
8	standard interpolation
9	standard interpolation (with separation of weights)
10	classical block interpolation (for use with nodal systems version only)
11	like 10, with diagonalized blocks
12	FF interpolation
13	FF1 interpolation
14	extended interpolation

When using this solver as a preconditioner, it is recommended to set the fail_max_iters flag to 0.

Example:

<default_linear_solver type="boomeramg">
  <max_iter>100</max_iter>
  <print_level>3</print_level>
  <tol>1e-05</tol>
  <max_levels>25</max_levels>
  <coarsen_type>-1</coarsen_type>
  <use_num_funcs>1</use_num_funcs>
  <relax_type>5</relax_type>
  <interp_type>6</interp_type>
  <strong_threshold>0.9</strong_threshold>
  <p_max_elmts>4</p_max_elmts>
  <num_sweeps>1</num_sweeps>
  <agg_interp_type>0</agg_interp_type>
  <agg_num_levels>0</agg_num_levels>
  <nodal>0</nodal>
  <do_jacobi>0</do_jacobi>
  <fail_max_iters>0</fail_max_iters>
</default_linear_solver>

Schur¶

The Schur solver is a linear solver that takes a 2x2 block structured matrix and solves the linear system by decomposing it in its Schur-complement form. Thus, instead of solving the following system directly,

\[ \left[\begin{array}{cc} A & B\\ C & D \end{array}\right]\left[\begin{array}{c} u\\ v \end{array}\right]=\left[\begin{array}{c} a\\ b \end{array}\right] \]

it first solves for v,

\[ Sv=b-CA^{-1}a \]

where \(S=D-CA^{-1}B\) is the Schur complement of A. Then the Schur solver solves for u,

\[ u=A^{-1}\left(a-Bv\right) \]

In order to use this solver, users need to select how they wish to solve the A-block, and how to solve the Schur complement. For the A-block users can choose any of the linear solvers listed above, either direct or iterative. For the Schur complement users can only choose a Krylov-based iterative solver (FGMRES or CG) because FEBio does not form the Schur complement explicitly.

This solver is often more effective as a preconditioner for FGMRES. In this case, the A-block and Schur complement only need to be solved approximately. See the examples section below for some example configurations. This solver requires a block-structured matrix. In order to generate the block structure you need to set the equation_scheme control parameter to 1 in the model input file.

parameter	description	default
print_level	The amount of information printed (0,1, or 2)	0
schur_block	Take the Schur complement of A (=0) or D (=1)	0
do_jacobi	Do Jacobi preconditioning	0
zero_D_block	Is the D-block zero?	0 (false)
A_solver	The solver for the A-block (outside the Schur complement)	user specified
schur_solver	The iterative solver for solving Schur complement	user specified
schur_A_solver	The solver for the A block in the Schur complement	(same as A_solver)

Accelerate¶

The Accelerate sparse solver is only available on Apple computers (). It uses the native Accelerate framework of Mac OS. It may be used as a direct solver or an iterative solver. Users should consult the Apple documentation to determine optimal settings for this solver, based on their application.

parameter	description	default
print_condition_number	Print the matrix condition number	0 (false)
iterative	Use iterative solver	0 (false)
factorization	Factorization method (0 to 6)	U:5=QR, S:4=DLTTPP
order_method	Ordering algorithm (0 to 2)	U:2=COLAMD, S:0=AMD
iterative_method	Iterative method (0 to 4)	4=LSMR

The factorization methods include: 0=Cholesky, 1=LDLT, 2=LDLTUnpivoted, 3=LDLTSBK, 4=LDLTTPP, 5=QR, 6=CholeskyAtA (). The ordering algorithms include: 0=AMD, 1=Metis, 2=COLAMD (). The iterative methods include: 0=ConjugateGradient, 1=GMRES, 2=DQGMRES, 3=FGMRES, 4=LSMR ().

Pardiso_64¶

The pardiso_64 solver is variant of the pardiso solver that can be used for very large problems since it uses 64 bit integers. Other than that, it is identical to the standard pardiso solver.

Example:

<default_linear_solver type="pardiso_64"/>