Appendix 19: GMRES Solver method

The Generalized Minimal RESidual (GMRES) method is an iterative method for the numerical solution of a system of linear equations developed by Y. Saad and Martin H. Schultz in 1986. The subroutine in Lagamine comes from the open source SPARKIT library.

The equation system is: \[[K]\underline{x}=\underline{F}\]

Where $[K]$ is the stiffness matrix, $\underline{x}$ is the displacement and $\underline{F}$ is the out of balanced forces.
One defines the residual, $\underline{r}^i$: \[\underline{r}^i =\underline{F}-[K]\underline{x}^i\] where $\underline{x}^k$ is an approximation of the solution at the iteration $k$.
At each iteration $i$, a Krylov subspace, $D_i$, is constructed by the Arnoldi’s method therefore $D_i$ is an orthogonal subspan. \[D_i = span(\underline{p}_1, [K]\underline{p}_1,\ldots,[K]^{i-1}\underline{p}_1) \text{ with } \underline{p}_1=\frac{\underline{r}^0}{\|\underline{r}^0\|_2}\] At each iteration, one searches $\underline{x}^i$ as: \[\underline{x}^i=\underline{x}^0+[V]^i\underline{y}^i\] Where $[V]^k$ is an array which contains the Krylov subspace vector and $\underline{y}^k$ is a $k$ vector which the residual norm: \[\underline{y}^i=\min{J(y)}=\min{\|\underline{F}-[K]\underline{x}^i\|}=min{\|\underline{F}-[K](\underline{x}^0+[V]^i\underline{y})\|}\] The method is considered converged if the residual norm is sufficiently diminished by the following criterion: \[\|\underline{r}^i\|<\varepsilon \|\underline{r}^0\|\] where $\varepsilon$ is a tolerance coefficient.
In case where the residual is inferior to a threshold, GMRES consider that the convergence is reached: \[\|\underline{r}^i\|<restol\] The GMRES method requires parameters:

• im: maximum dimension of the Krylov subspace
• maxit: maximum of iteration
• $\varepsilon$: tolerance coefficient
• restol: residual threshold

Preconditioning is a technique to reach quicker the convergence. The convergence depends on the distribution and the value of the eigenvalues of the stiffness matrix, $[K]$.
Suppose that $[M]$ is symmetric, positive-definite matrix that approximates $[K]$, but easier to invert. The equation system can be indirectly solved by solving: \[ [M]^{-1}[K]\underline{x}=[M]^{-1}\underline{F}\] If the set of eigenvalue of $[M]^{-1}[K]$ is lower and better clustered than the set of $[K]$, then the number of iterations to solve without iteration will be lower than the number of iterations with preconditioning.
The simplest ways of defining a preconditioner is to perform an incomplete LU factorization of the original matrix $[K]$. This entails a decomposition of the form: \[[K]=[L]_k[U]_k-[R]\] where $[L]_k$ and $[U]_k$ have nonzero structure as the lower part and upper parts of $[K]$ respectively and $[R]$ is the residual error of the factorization therefore $[M]=[L]_k[U]_k$.
The fill-in depends on the geometrical tolerance to avoid too large bandwidth for $[L]_k$ and $[U]_k$. The fill-in depends if the component, which does not belong to the diagonal is not inferior to a tolerance value. The preconditioning required two parameters:

• lfil: degree of fill in
• droptol: acceptance limit for the value of the components which are out of the diagonal of the stiffness matrix.

KNSYM must be equal to 6 at the 1st line and the 15th column in the loading file. Then the parameters for the iterative method must be entered on the second line (this line is added in the loading file) with the following format:

• Higher are IM and LFIL, lower the number of iterations before convergence is. However the memory cost is more important.
• Higher EPS is, lower the number of iterations before convergence for the iterative solver is. However if the residual is not enough low, the out-of-balanced forces diminish slower and the Newton-Raphson can converge slower.
• Higher DROTOL is, lower the memory cost of the preconditioning matrix is. However the number of iterations before convergence for the iterative solver increases.