View source for Solving sparse systems

There are many methods available for solving sparse systems. We compare some of them here.
<figure id="fig:matrix1">
[[File:matrix.png|300px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Matrix of the discretized PDE. </caption>]]
</figure>
Mathematica has the following methods available (https://reference.wolfram.com/language/ref/LinearSolve.html#DetailsAndOptions)
* direct: banded, cholesky, multifrontal (direct sparse LU)
* iterative: Krylov

Matlab has the following methods:
* direct: https://www.mathworks.com/help/matlab/ref/mldivide.html#bt42omx_head
* iterative: https://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#brzoiix, including bicgstab, gmres 

Eigen has the following methods: (https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
* direct: sparse LU
* iterative: bicgstab, cg

Solving a simple sparse system $A x = b$ for steady space of heat equation in 1d with $n$ nodes, results in a matrix shown in <xr id="fig:matrix1"/>.

The following timings of solvers are given in seconds:
{| class="wikitable"
|-
! $n = 10^6$
! Matlab
! Mathematica
! Eigen
|-
! Banded
| 0.16
| 0.28
| 0.04
|-
! SparseLU
| /
| 1.73
| 0.82
|-
! BICGStab / Krylov
| 0.33
| 0.39
| 0.53
|}

Incomplete LU preconditioner was used for BICGStab.
Without the preconditioner BICGStab does not converge.

==Parallel execution ==
BICGStab can be run in parallel, as explain in the general parallelism: https://eigen.tuxfamily.org/dox/TopicMultiThreading.html, and specifically

'''"When using sparse matrices, best performance is achied for a row-major sparse matrix format.<br>Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled".'''

Eigen uses number of threads specified my OopenMP, unless <code>Eigen::setNbThreads(n);</code> was called.
Minimal working example:

<figure id="fig:eigen_par">
[[File:eigen_parallel_2.png|1200px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Memory and CPU usage. C stands for construction of the system, L stands for the calculation of ILUT preconditioner and S for BICGStab iteration.</caption>]]
</figure>

<syntaxhighlight lang="c++" line>
#include <iostream>
#include <vector>
#include "Eigen/Sparse"
#include "Eigen/IterativeLinearSolvers"

using namespace std;
using namespace Eigen;

int main(int argc, char* argv[]) {

    assert(argc == 2 && "Second argument is size of the system.");
    stringstream ss(argv[1]);
    int n;
    ss >> n;
    cout << "n = " << n << endl;

    // Fill the matrix
    VectorXd b = VectorXd::Ones(n) / n / n;
    b(0) = b(n-1) = 1;
    SparseMatrix<double, RowMajor> A(n, n);
    A.reserve(vector<int>(n, 3));  // 3 per row
    for (int i = 0; i < n-1; ++i) {
        A.insert(i, i) = -2;
        A.insert(i, i+1) = 1;
        A.insert(i+1, i) = 1;
    }
    A.coeffRef(0, 0) = 1;
    A.coeffRef(0, 1) = 0;
    A.coeffRef(n-1, n-2) = 0;
    A.coeffRef(n-1, n-1) = 1;

    // Solve the system
    BiCGSTAB<SparseMatrix<double, RowMajor>, IncompleteLUT<double>> solver;
    solver.setTolerance(1e-10);
    solver.setMaxIterations(1000);
    solver.compute(A);
    VectorXd x = solver.solve(b);
    cout << "#iterations:     " << solver.iterations() << endl;
    cout << "estimated error: " << solver.error()      << endl;
    cout << "sol: " << x.head(6).transpose() << endl;

    return 0;
}
</syntaxhighlight>

was compiled with <code>g++ -o parallel_solve -O3 -fopenmp solver_test_parallel.cpp</code>.
<xr id="fig:eigen_par"/> was produced when the program above was run as <code>./parallel_solve 10000000</code>.

==ILUT factorization==
This is a method to compute a general algebraic preconditioner which has the ability to control (to some degree) the memory and time complexity of the solution.
It is described in a paper by  

  Saad, Yousef. "ILUT: A dual threshold incomplete LU factorization." Numerical linear algebra with applications 1.4 (1994): 387-402.


It has two parameters, tolerance $\tau$ and fill factor $f$. Elements in $L$ and $U$ factors that would be below relative tolerance $\tau$ (say, to 2-norm of this row) are dropped. 
Then, only $f$ elements per row in $L$ and $U$ are allowed and the maximum $f$ (if they exist) are taken. The diagonal elements are always kept. This means that the number of non-zero
elements in the rows of the preconditioner can not exceed $2f+1$. 

Greater parameter $f$ means slover but better $ILUT$ factorization. If $\tau$ is small, having a very large $f$ means getting a complete $LU$ factorization. Typical values are $5, 10, 20, 50$.
Lower $\tau$ values mean that move small elements are kept, resulting in a more precise, but longer computation. With a large $f$, the factorization approaches the direct $LU$ factorization as $\tau \to 0$. Typical values are from 1e-2 
to 1e-6.

Whether the method converges or diverges for given parameters is very sudden and should be determined for each apllication separately. An example of 
such a convergence graph for the diffusion equation is presented in figures below:

[[File:bicgstab_err.png|600px]][[File:bicgstab_conv.png|600px]][[File:bicgstab_iter.png|600px]]

We can see that low BiCGSTAB error estimate coincides with the range of actual small error compared to the analytical solution as well as with the range of lesser number of iterations.

==Herzian contact ==
The matrix of the [[Hertzian contact]] problem is shown below:
 [[File:matrix_hertzian.png|1000px]]

The condition number iz arther large and grows superlinearly with size in the example below.
 [[File:cond.png|800px]]