Template Numerical Library version\ main:eacc201d
Loading...
Searching...
No Matches
Linear solvers tutorial

Introduction

Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the followiing iterative methods:

  1. Stationary methods
    1. Jacobi method (TNL::Solvers::Linear::Jacobi)
    2. Successive-overrelaxation method, SOR (TNL::Solvers::Linear::SOR)
  2. Krylov subspace methods
    1. Conjugate gradient method, CG (TNL::Solvers::Linear::CG)
    2. Biconjugate gradient stabilized method, BICGStab (TNL::Solvers::Linear::BICGStab)
    3. Biconjugate gradient stabilized method, BICGStab(l) (TNL::Solvers::Linear::BICGStabL)
    4. Transpose-free quasi-minimal residual method, TFQMR (TNL::Solvers::Linear::TFQMR)
    5. Generalized minimal residual method, GMRES (TNL::Solvers::Linear::GMRES) with various methods of orthogonalization
      1. Classical Gramm-Schmidt, CGS
      2. Classical Gramm-Schmidt with reorthogonalization, CGSR
      3. Modified Gramm-Schmidt, MGS
      4. Modified Gramm-Schmidt with reorthogonalization, MGSR
      5. Compact WY form of the Householder reflections, CWY

The iterative solvers (not the stationary solvers like TNL::Solvers::Linear::Jacobi and TNL::Solvers::Linear::SOR) can be combined with the following preconditioners:

  1. Diagonal or Jacobi - TNL::Solvers::Linear::Preconditioners::Diagonal
  2. ILU (Incomplete LU) - CPU only currently
    1. ILU(0) TNL::Solvers::Linear::Preconditioners::ILU0
    2. ILUT (ILU with thresholding) TNL::Solvers::Linear::Preconditioners::ILUT

Iterative solvers of linear systems

Basic setup

All iterative solvers for linear systems can be found in the namespace TNL::Solvers::Linear. The following example shows the use the iterative solvers:

1#include <iostream>
2#include <memory>
3#include <TNL/Algorithms/ParallelFor.h>
4#include <TNL/Matrices/SparseMatrix.h>
5#include <TNL/Devices/Host.h>
6#include <TNL/Devices/Cuda.h>
7#include <TNL/Solvers/Linear/TFQMR.h>
8
9template< typename Device >
10void iterativeLinearSolverExample()
11{
12 /***
13 * Set the following matrix (dots represent zero matrix elements):
14 *
15 * / 2.5 -1 . . . \
16 * | -1 2.5 -1 . . |
17 * | . -1 2.5 -1. . |
18 * | . . -1 2.5 -1 |
19 * \ . . . -1 2.5 /
20 */
23 const int size( 5 );
24 auto matrix_ptr = std::make_shared< MatrixType >();
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
28 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
29 const int rowIdx = row.getRowIndex();
30 if( rowIdx == 0 )
31 {
32 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
33 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
34 }
35 else if( rowIdx == size - 1 )
36 {
37 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else
41 {
42 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
43 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
44 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
45 }
46 };
47
48 /***
49 * Set the matrix elements.
50 */
51 matrix_ptr->forAllRows( f );
52 std::cout << *matrix_ptr << std::endl;
53
54 /***
55 * Set the right-hand side vector.
56 */
57 Vector x( size, 1.0 );
58 Vector b( size );
59 matrix_ptr->vectorProduct( x, b );
60 x = 0.0;
61 std::cout << "Vector b = " << b << std::endl;
62
63 /***
64 * Solve the linear system.
65 */
67 LinearSolver solver;
68 solver.setMatrix( matrix_ptr );
69 solver.setConvergenceResidue( 1.0e-6 );
70 solver.solve( b, x );
71 std::cout << "Vector x = " << x << std::endl;
72}
73
74int main( int argc, char* argv[] )
75{
76 std::cout << "Solving linear system on host: " << std::endl;
77 iterativeLinearSolverExample< TNL::Devices::Sequential >();
78
79#ifdef HAVE_CUDA
80 std::cout << "Solving linear system on CUDA device: " << std::endl;
81 iterativeLinearSolverExample< TNL::Devices::Cuda >();
82#endif
83}
#define __cuda_callable__
Definition: CudaCallable.h:22
Vector extends Array with algebraic operations.
Definition: Vector.h:40
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition: SparseMatrix.h:57
void setMatrix(const MatrixPointer &matrix)
Set the matrix of the linear system.
Definition: LinearSolver.h:127
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition: TFQMR.h:26
T endl(T... args)

In this example we solve a linear system \( A \vec x = \vec b \) where

\[ A = \left( \begin{array}{cccc} 2.5 & -1 & & & \\ -1 & 2.5 & -1 & & \\ & -1 & 2.5 & -1 & \\ & & -1 & 2.5 & -1 \\ & & & -1 & 2.5 \\ \end{array} \right) \]

The right-hand side vector \(\vec b \) is set to \(( 1.5, 0.5, 0.5, 0.5, 1.5 )^T \) so that the exact solution is \( \vec x = ( 1, 1, 1, 1, 1 )^T\). The matrix elements of \(A $\) is set on the lines 12-51 by the means of the method TNL::Matrices::SparseMatrix::forAllElements. In this example, we use the sparse matrix but any other matrix type can be used as well (see the namespace TNL::Matrices). Next we set the solution vector \( \vec x = ( 1, 1, 1, 1, 1 )^T\) (line 57) and multiply it with matrix \( A \) to get the right-hand side vector \(\vec b\) (lines 58-59). Finally, we reset the vector \(\vec x \) to zero vector.

To solve the linear system, we use TFQMR method (line 66), as an example. Other solvers can be used as well (see the namespace TNL::Solvers::Linear). The solver needs only one template parameter which is the matrix type. Next we create an instance of the solver (line 67 ) and set the matrix of the linear system (line 68). Note, that matrix is passed to the solver as a shared smart pointer (std::shared_ptr). This is why we created an instance of the smart pointer on the line 24 instead of the sparse matrix itself. On the line 69, we set the value of residue under which the convergence occur and the solver is stopped. It serves as a convergence criterion. The solver is executed on the line 70 by calling the method TNL::Solvers::Linear::LinearSolver::solve. The method accepts the right-hand side vector \( \vec b\) and the solution vector \( \vec x\).

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]

Setup with a solver monitor

Solution of large linear systems may take a lot of time. In such situations, it is useful to be able to monitor the convergence of the solver of the solver status in general. For this purpose, TNL offers solver monitors. The solver monitor prints (or somehow visualizes) the number of iterations, the residue of the current solution approximation or some other metrics. Sometimes such information is printed after each iteration or after every ten iterations. The problem of this approach is the fact that one iteration of the solver may take only few milliseconds but also several minutes. In the former case, the monitor creates overwhelming amount of output which may even slowdown the solver. In the later case, the user waits long time for update of the solver status. The monitor in TNL rather runs in separate thread and it refreshes the status of the solver in preset time periods. The use of the iterative solver monitor is demonstrated in the following example.

1#include <iostream>
2#include <memory>
3#include <chrono>
4#include <thread>
5#include <TNL/Algorithms/ParallelFor.h>
6#include <TNL/Matrices/SparseMatrix.h>
7#include <TNL/Devices/Sequential.h>
8#include <TNL/Devices/Cuda.h>
9#include <TNL/Solvers/Linear/Jacobi.h>
10
11template< typename Device >
12void iterativeLinearSolverExample()
13{
14 /***
15 * Set the following matrix (dots represent zero matrix elements):
16 *
17 * / 2.5 -1 . . . \
18 * | -1 2.5 -1 . . |
19 * | . -1 2.5 -1. . |
20 * | . . -1 2.5 -1 |
21 * \ . . . -1 2.5 /
22 */
25 const int size( 5 );
26 auto matrix_ptr = std::make_shared< MatrixType >();
27 matrix_ptr->setDimensions( size, size );
28 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
29
30 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
31 const int rowIdx = row.getRowIndex();
32 if( rowIdx == 0 )
33 {
34 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
35 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
36 }
37 else if( rowIdx == size - 1 )
38 {
39 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
40 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
41 }
42 else
43 {
44 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
45 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
46 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
47 }
48 };
49
50 /***
51 * Set the matrix elements.
52 */
53 matrix_ptr->forAllRows( f );
54 std::cout << *matrix_ptr << std::endl;
55
56 /***
57 * Set the right-hand side vector.
58 */
59 Vector x( size, 1.0 );
60 Vector b( size );
61 matrix_ptr->vectorProduct( x, b );
62 x = 0.0;
63 std::cout << "Vector b = " << b << std::endl;
64
65 /***
66 * Setup solver of the linear system.
67 */
69 LinearSolver solver;
70 solver.setMatrix( matrix_ptr );
71 solver.setOmega( 0.0005 );
72
73 /***
74 * Setup monitor of the iterative solver.
75 */
76 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
77 IterativeSolverMonitorType monitor;
78 TNL::Solvers::SolverMonitorThread monitorThread(monitor);
79 monitor.setRefreshRate(10); // refresh rate in milliseconds
80 monitor.setVerbose(1);
81 monitor.setStage( "Jacobi stage:" );
82 solver.setSolverMonitor(monitor);
83 solver.setConvergenceResidue( 1.0e-6 );
84 solver.solve( b, x );
85 monitor.stopMainLoop();
86 std::cout << "Vector x = " << x << std::endl;
87}
88
89int main( int argc, char* argv[] )
90{
91 std::cout << "Solving linear system on host: " << std::endl;
92 iterativeLinearSolverExample< TNL::Devices::Sequential >();
93
94#ifdef HAVE_CUDA
95 std::cout << "Solving linear system on CUDA device: " << std::endl;
96 iterativeLinearSolverExample< TNL::Devices::Cuda >();
97#endif
98}
Iterative solver of linear systems based on the Jacobi method.
Definition: Jacobi.h:26
A RAII wrapper for launching the SolverMonitor's main loop in a separate thread.
Definition: SolverMonitor.h:137

On the lines 1-70, we setup the same linear system as in the previous example, we create an instance of the Jacobi solver and we pass the matrix of the linear system to the solver. On the line 71, we set the relaxation parameter \( \omega \) of the Jacobi solver to 0.0005 (TNL::Solvers::Linear::Jacobi). The reason is to slowdown the convergence because we want to see some iterations in this example. Next we create an instance of the solver monitor (lines 76 and 77) and we create a special thread for the monitor (line 78, TNL::Solvers::SolverMonitorThread ). We set the refresh rate of the monitor to 10 milliseconds (line 79, TNL::Solvers::SolverMonitor::setRefreshRate). We set a verbosity of the monitor to 1 (line 80 TNL::Solvers::IterativeSolverMonitor::setVerbose). Next we set a name of the solver stage (line 81, TNL::Solvers::IterativeSolverMonitor::setStage). The monitor stages serve for distinguishing between different phases or stages of more complex solvers (for example when the linear system solver is embedded into a time dependent PDE solver). Next we connect the solver with the monitor (line 82, TNL::Solvers::IterativeSolver::setSolverMonitor) and we set the convergence criterion based on the value of the residue (line 83). Finally we start the solver (line 84, TNL::Solvers::Linear::Jacobi::solve) and when the solver finishes we have to stop the monitor (line 84, TNL::Solvers::SolverMonitor::stopMainLoop).

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Jacobi stage: ITER: 10835 RES: 0.13673
Jacobi stage: ITER: 19439 RES: 0.036484
Jacobi stage: ITER: 29327 RES: 0.0079889
Jacobi stage: ITER: 38839 RES: 0.0018533
Jacobi stage: ITER: 47531 RES: 0.00048767
Jacobi stage: ITER: 56887 RES: 0.00011588
Jacobi stage: ITER: 66807 RES: 2.525e-05
Jacobi stage: ITER: 76279 RES: 5.8938e-06
Jacobi stage: ITER: 84849 RES: 1.5797e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Jacobi stage: ITER: 2559 RES: 0.51931
Jacobi stage: ITER: 4345 RES: 0.37751
Jacobi stage: ITER: 5951 RES: 0.29143
Jacobi stage: ITER: 7855 RES: 0.21655
Jacobi stage: ITER: 8733 RES: 0.18903
Jacobi stage: ITER: 9523 RES: 0.16742
Jacobi stage: ITER: 10731 RES: 0.13902
Jacobi stage: ITER: 11323 RES: 0.12693
Jacobi stage: ITER: 13117 RES: 0.096319
Jacobi stage: ITER: 14465 RES: 0.078303
Jacobi stage: ITER: 15381 RES: 0.068025
Jacobi stage: ITER: 16023 RES: 0.061656
Jacobi stage: ITER: 16491 RES: 0.057379
Jacobi stage: ITER: 17091 RES: 0.052327
Jacobi stage: ITER: 18235 RES: 0.043895
Jacobi stage: ITER: 20177 RES: 0.032564
Jacobi stage: ITER: 21927 RES: 0.024896
Jacobi stage: ITER: 22523 RES: 0.022718
Jacobi stage: ITER: 23127 RES: 0.020705
Jacobi stage: ITER: 24729 RES: 0.016184
Jacobi stage: ITER: 26495 RES: 0.012342
Jacobi stage: ITER: 27233 RES: 0.011016
Jacobi stage: ITER: 28463 RES: 0.0091226
Jacobi stage: ITER: 29287 RES: 0.0080381
Jacobi stage: ITER: 30855 RES: 0.0063176
Jacobi stage: ITER: 32635 RES: 0.0048063
Jacobi stage: ITER: 33899 RES: 0.0039582
Jacobi stage: ITER: 35087 RES: 0.003298
Jacobi stage: ITER: 35963 RES: 0.0028828
Jacobi stage: ITER: 36675 RES: 0.0025841
Jacobi stage: ITER: 37519 RES: 0.0022699
Jacobi stage: ITER: 38371 RES: 0.0019915
Jacobi stage: ITER: 39511 RES: 0.0016716
Jacobi stage: ITER: 40843 RES: 0.0013623
Jacobi stage: ITER: 42235 RES: 0.0011001
Jacobi stage: ITER: 43447 RES: 0.00091319
Jacobi stage: ITER: 44513 RES: 0.00077503
Jacobi stage: ITER: 45697 RES: 0.00064615
Jacobi stage: ITER: 47583 RES: 0.00048379
Jacobi stage: ITER: 49153 RES: 0.00038001
Jacobi stage: ITER: 50827 RES: 0.00029394
Jacobi stage: ITER: 52675 RES: 0.0002213
Jacobi stage: ITER: 54689 RES: 0.00016237
Jacobi stage: ITER: 55971 RES: 0.00013339
Jacobi stage: ITER: 57907 RES: 9.9074e-05
Jacobi stage: ITER: 59955 RES: 7.2334e-05
Jacobi stage: ITER: 61491 RES: 5.7132e-05
Jacobi stage: ITER: 63251 RES: 4.3598e-05
Jacobi stage: ITER: 64419 RES: 3.6438e-05
Jacobi stage: ITER: 65655 RES: 3.0137e-05
Jacobi stage: ITER: 67387 RES: 2.3097e-05
Jacobi stage: ITER: 68763 RES: 1.8697e-05
Jacobi stage: ITER: 70419 RES: 1.4498e-05
Jacobi stage: ITER: 72167 RES: 1.1084e-05
Jacobi stage: ITER: 74107 RES: 8.2278e-06
Jacobi stage: ITER: 76063 RES: 6.0926e-06
Jacobi stage: ITER: 77775 RES: 4.6838e-06
Jacobi stage: ITER: 78875 RES: 3.9557e-06
Jacobi stage: ITER: 80203 RES: 3.2257e-06
Jacobi stage: ITER: 83407 RES: 1.972e-06
Jacobi stage: ITER: 86595 RES: 1.2085e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]

The monitoring of the solver can be improved by time elapsed since the beginning of the computation as demonstrated in the following example:

1#include <iostream>
2#include <memory>
3#include <chrono>
4#include <thread>
5#include <TNL/Timer.h>
6#include <TNL/Algorithms/ParallelFor.h>
7#include <TNL/Matrices/SparseMatrix.h>
8#include <TNL/Devices/Sequential.h>
9#include <TNL/Devices/Cuda.h>
10#include <TNL/Solvers/Linear/Jacobi.h>
11
12template< typename Device >
13void iterativeLinearSolverExample()
14{
15 /***
16 * Set the following matrix (dots represent zero matrix elements):
17 *
18 * / 2.5 -1 . . . \
19 * | -1 2.5 -1 . . |
20 * | . -1 2.5 -1. . |
21 * | . . -1 2.5 -1 |
22 * \ . . . -1 2.5 /
23 */
26 const int size( 5 );
27 auto matrix_ptr = std::make_shared< MatrixType >();
28 matrix_ptr->setDimensions( size, size );
29 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
30
31 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
32 const int rowIdx = row.getRowIndex();
33 if( rowIdx == 0 )
34 {
35 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
36 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
37 }
38 else if( rowIdx == size - 1 )
39 {
40 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
41 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
42 }
43 else
44 {
45 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
46 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
47 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
48 }
49 };
50
51 /***
52 * Set the matrix elements.
53 */
54 matrix_ptr->forAllRows( f );
55 std::cout << *matrix_ptr << std::endl;
56
57 /***
58 * Set the right-hand side vector.
59 */
60 Vector x( size, 1.0 );
61 Vector b( size );
62 matrix_ptr->vectorProduct( x, b );
63 x = 0.0;
64 std::cout << "Vector b = " << b << std::endl;
65
66 /***
67 * Setup solver of the linear system.
68 */
70 LinearSolver solver;
71 solver.setMatrix( matrix_ptr );
72 solver.setOmega( 0.0005 );
73
74 /***
75 * Setup monitor of the iterative solver.
76 */
77 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
78 IterativeSolverMonitorType monitor;
79 TNL::Solvers::SolverMonitorThread mmonitorThread(monitor);
80 monitor.setRefreshRate(10); // refresh rate in milliseconds
81 monitor.setVerbose(1);
82 monitor.setStage( "Jacobi stage:" );
83 TNL::Timer timer;
84 monitor.setTimer( timer );
85 timer.start();
86 solver.setSolverMonitor(monitor);
87 solver.setConvergenceResidue( 1.0e-6 );
88 solver.solve( b, x );
89 monitor.stopMainLoop();
90 std::cout << "Vector x = " << x << std::endl;
91}
92
93int main( int argc, char* argv[] )
94{
95 std::cout << "Solving linear system on host: " << std::endl;
96 iterativeLinearSolverExample< TNL::Devices::Sequential >();
97
98#ifdef HAVE_CUDA
99 std::cout << "Solving linear system on CUDA device: " << std::endl;
100 iterativeLinearSolverExample< TNL::Devices::Cuda >();
101#endif
102}
Class for real time, CPU time and CPU cycles measuring.
Definition: Timer.h:28
void start()
Starts timer.
Definition: Timer.hpp:49

The only changes happen on lines 83-85 where we create an instance of TNL timer (line 83, TNL::Timer), connect it with the monitor (line 84, TNL::Solvers::SolverMonitor::setTimer) and start the timer (line 85, TNL::Timer::start).

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
ELA:0.002284
ELA: 0.10314 Jacobi stage: ITER: 17851 RES: 0.046562
ELA: 0.20561 Jacobi stage: ITER: 36061 RES: 0.0028388
ELA: 0.30893 Jacobi stage: ITER: 54655 RES: 0.00016327
ELA: 0.41227 Jacobi stage: ITER: 69049 RES: 1.7888e-05
ELA: 0.5156 Jacobi stage: ITER: 86081 RES: 1.3073e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
ELA:0.007485
ELA: 0.1076 Jacobi stage: ITER: 1533 RES: 0.6462
ELA: 0.20771 Jacobi stage: ITER: 1979 RES: 0.58487
ELA: 0.31083 Jacobi stage: ITER: 3921 RES: 0.40545
ELA: 0.41416 Jacobi stage: ITER: 6585 RES: 0.26373
ELA: 0.51427 Jacobi stage: ITER: 9083 RES: 0.17916
ELA: 0.61749 Jacobi stage: ITER: 11311 RES: 0.12716
ELA: 0.72082 Jacobi stage: ITER: 13155 RES: 0.095788
ELA: 0.82415 Jacobi stage: ITER: 14187 RES: 0.081744
ELA: 0.92428 Jacobi stage: ITER: 14971 RES: 0.072469
ELA: 1.0244 Jacobi stage: ITER: 15857 RES: 0.063229
ELA: 1.1245 Jacobi stage: ITER: 18353 RES: 0.043093
ELA: 1.2246 Jacobi stage: ITER: 20827 RES: 0.029478
ELA: 1.3247 Jacobi stage: ITER: 23147 RES: 0.020642
ELA: 1.4249 Jacobi stage: ITER: 24515 RES: 0.01673
ELA: 1.525 Jacobi stage: ITER: 25519 RES: 0.014339
ELA: 1.6275 Jacobi stage: ITER: 26987 RES: 0.011444
ELA: 1.7309 Jacobi stage: ITER: 29185 RES: 0.0081625
ELA: 1.8311 Jacobi stage: ITER: 31229 RES: 0.0059631
ELA: 1.9312 Jacobi stage: ITER: 33051 RES: 0.0045088
ELA: 2.0313 Jacobi stage: ITER: 34519 RES: 0.0035986
ELA: 2.1316 Jacobi stage: ITER: 35245 RES: 0.0032179
ELA: 2.2317 Jacobi stage: ITER: 36419 RES: 0.0026877
ELA: 2.3318 Jacobi stage: ITER: 37283 RES: 0.0023537
ELA: 2.4342 Jacobi stage: ITER: 39021 RES: 0.0018017
ELA: 2.5343 Jacobi stage: ITER: 39923 RES: 0.0015691
ELA: 2.6344 Jacobi stage: ITER: 41963 RES: 0.001147
ELA: 2.7345 Jacobi stage: ITER: 43611 RES: 0.00089048
ELA: 2.8346 Jacobi stage: ITER: 45205 RES: 0.00069688
ELA: 2.9347 Jacobi stage: ITER: 47419 RES: 0.00049613
ELA: 3.0348 Jacobi stage: ITER: 49499 RES: 0.00036045
ELA: 3.135 Jacobi stage: ITER: 51519 RES: 0.0002643
ELA: 3.2351 Jacobi stage: ITER: 53607 RES: 0.00019178
ELA: 3.3375 Jacobi stage: ITER: 55103 RES: 0.00015241
ELA: 3.4408 Jacobi stage: ITER: 56485 RES: 0.00012322
ELA: 3.5442 Jacobi stage: ITER: 58015 RES: 9.7444e-05
ELA: 3.6443 Jacobi stage: ITER: 60087 RES: 7.0882e-05
ELA: 3.7444 Jacobi stage: ITER: 61251 RES: 5.9277e-05
ELA: 3.8445 Jacobi stage: ITER: 62371 RES: 4.9908e-05
ELA: 3.9446 Jacobi stage: ITER: 63869 RES: 3.9638e-05
ELA: 4.0447 Jacobi stage: ITER: 64729 RES: 3.4733e-05
ELA: 4.1448 Jacobi stage: ITER: 65503 RES: 3.0849e-05
ELA: 4.245 Jacobi stage: ITER: 67159 RES: 2.3921e-05
ELA: 4.3451 Jacobi stage: ITER: 68279 RES: 2.014e-05
ELA: 4.4452 Jacobi stage: ITER: 69249 RES: 1.7347e-05
ELA: 4.5453 Jacobi stage: ITER: 70095 RES: 1.5238e-05
ELA: 4.6454 Jacobi stage: ITER: 70785 RES: 1.3701e-05
ELA: 4.7456 Jacobi stage: ITER: 72269 RES: 1.0908e-05
ELA: 4.8457 Jacobi stage: ITER: 73397 RES: 9.1731e-06
ELA: 4.9459 Jacobi stage: ITER: 74659 RES: 7.559e-06
ELA: 5.046 Jacobi stage: ITER: 76095 RES: 6.0627e-06
ELA: 5.1461 Jacobi stage: ITER: 77279 RES: 5.0546e-06
ELA: 5.2462 Jacobi stage: ITER: 78967 RES: 3.9002e-06
ELA: 5.3463 Jacobi stage: ITER: 80669 RES: 3.002e-06
ELA: 5.4464 Jacobi stage: ITER: 81803 RES: 2.5229e-06
ELA: 5.5466 Jacobi stage: ITER: 83483 RES: 1.9491e-06
ELA: 5.6467 Jacobi stage: ITER: 84769 RES: 1.5992e-06
ELA: 5.7468 Jacobi stage: ITER: 86255 RES: 1.2732e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]

Setup with preconditioner

Preconditioners of iterative solvers can significantly improve the performance of the solver. In the case of the linear systems, they are used mainly with the Krylov subspace methods. Preconditioners cannot be used with the starionary methods (TNL::Solvers::Linear::Jacobi and TNL::Solvers::Linear::SOR). The following example shows how to setup an iterative solver of linear systems with preconditioning.

1#include <iostream>
2#include <memory>
3#include <TNL/Algorithms/ParallelFor.h>
4#include <TNL/Matrices/SparseMatrix.h>
5#include <TNL/Devices/Host.h>
6#include <TNL/Devices/Cuda.h>
7#include <TNL/Solvers/Linear/TFQMR.h>
8#include <TNL/Solvers/Linear/Preconditioners/Diagonal.h>
9
10template< typename Device >
11void iterativeLinearSolverExample()
12{
13 /***
14 * Set the following matrix (dots represent zero matrix elements):
15 *
16 * / 2.5 -1 . . . \
17 * | -1 2.5 -1 . . |
18 * | . -1 2.5 -1. . |
19 * | . . -1 2.5 -1 |
20 * \ . . . -1 2.5 /
21 */
24 const int size( 5 );
25 auto matrix_ptr = std::make_shared< MatrixType >();
26 matrix_ptr->setDimensions( size, size );
27 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
28
29 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
30 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 )
32 {
33 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
34 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
35 }
36 else if( rowIdx == size - 1 )
37 {
38 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
39 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
40 }
41 else
42 {
43 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
44 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
45 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
46 }
47 };
48
49 /***
50 * Set the matrix elements.
51 */
52 matrix_ptr->forAllRows( f );
53 std::cout << *matrix_ptr << std::endl;
54
55 /***
56 * Set the right-hand side vector.
57 */
58 Vector x( size, 1.0 );
59 Vector b( size );
60 matrix_ptr->vectorProduct( x, b );
61 x = 0.0;
62 std::cout << "Vector b = " << b << std::endl;
63
64 /***
65 * Solve the linear system using diagonal (Jacobi) preconditioner.
66 */
69 auto preconditioner_ptr = std::make_shared< Preconditioner >();
70 preconditioner_ptr->update( matrix_ptr );
71 LinearSolver solver;
72 solver.setMatrix( matrix_ptr );
73 solver.setPreconditioner( preconditioner_ptr );
74 solver.setConvergenceResidue( 1.0e-6 );
75 solver.solve( b, x );
76 std::cout << "Vector x = " << x << std::endl;
77}
78
79int main( int argc, char* argv[] )
80{
81 std::cout << "Solving linear system on host: " << std::endl;
82 iterativeLinearSolverExample< TNL::Devices::Sequential >();
83
84#ifdef HAVE_CUDA
85 std::cout << "Solving linear system on CUDA device: " << std::endl;
86 iterativeLinearSolverExample< TNL::Devices::Cuda >();
87#endif
88}
Diagonal (Jacobi) preconditioner for iterative solvers of linear systems.
Definition: Diagonal.h:29

In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen on the lines (68-72). Similar to the matrix of the linear system, the preconditioner is passed to the solver by the means of smart shared pointer (std::shared_ptr). The instance is created on the lines 68 and 69. Next we have to initialize the preconditioner (line 70, TNL::Solvers::Linear::Preconditioners::Preconditioner::update). The method update has to be called everytime the matrix of the linear system changes. This is important for example when solving time dependent PDEs but it does not happen in this example. Finally, we need to connect the solver with the preconditioner (line 73, TNL::Solvers::Linear::LinearSolver).

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]

Choosing the solver and preconditioner type at runtime

When developing a numerical solver, one often has to search for a combination of various methods and algorithms that fit given requirements the best. To make this easier, TNL offers choosing the type of both linear solver and preconditioner at runtime by means of functions TNL::Solvers::getLinearSolver and TNL::Solvers::getPreconditioner. The following example shows how to use these functions:

1#include <iostream>
2#include <memory>
3#include <TNL/Algorithms/ParallelFor.h>
4#include <TNL/Matrices/SparseMatrix.h>
5#include <TNL/Devices/Host.h>
6#include <TNL/Devices/Cuda.h>
7#include <TNL/Solvers/LinearSolverTypeResolver.h>
8
9template< typename Device >
10void iterativeLinearSolverExample()
11{
12 /***
13 * Set the following matrix (dots represent zero matrix elements):
14 *
15 * / 2.5 -1 . . . \
16 * | -1 2.5 -1 . . |
17 * | . -1 2.5 -1. . |
18 * | . . -1 2.5 -1 |
19 * \ . . . -1 2.5 /
20 */
23 const int size( 5 );
24 auto matrix_ptr = std::make_shared< MatrixType >();
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
28 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
29 const int rowIdx = row.getRowIndex();
30 if( rowIdx == 0 )
31 {
32 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
33 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
34 }
35 else if( rowIdx == size - 1 )
36 {
37 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else
41 {
42 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
43 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
44 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
45 }
46 };
47
48 /***
49 * Set the matrix elements.
50 */
51 matrix_ptr->forAllRows( f );
52 std::cout << *matrix_ptr << std::endl;
53
54 /***
55 * Set the right-hand side vector.
56 */
57 Vector x( size, 1.0 );
58 Vector b( size );
59 matrix_ptr->vectorProduct( x, b );
60 x = 0.0;
61 std::cout << "Vector b = " << b << std::endl;
62
63 /***
64 * Solve the linear system using diagonal (Jacobi) preconditioner.
65 */
66 auto solver_ptr = TNL::Solvers::getLinearSolver< MatrixType >( "tfqmr" );
67 auto preconditioner_ptr = TNL::Solvers::getPreconditioner< MatrixType >( "diagonal");
68 preconditioner_ptr->update( matrix_ptr );
69 solver_ptr->setMatrix( matrix_ptr );
70 solver_ptr->setPreconditioner( preconditioner_ptr );
71 solver_ptr->setConvergenceResidue( 1.0e-6 );
72 solver_ptr->solve( b, x );
73 std::cout << "Vector x = " << x << std::endl;
74}
75
76int main( int argc, char* argv[] )
77{
78 std::cout << "Solving linear system on host: " << std::endl;
79 iterativeLinearSolverExample< TNL::Devices::Sequential >();
80
81#ifdef HAVE_CUDA
82 std::cout << "Solving linear system on CUDA device: " << std::endl;
83 iterativeLinearSolverExample< TNL::Devices::Cuda >();
84#endif
85}

We still stay with the same problem and the only changes can be seen on lines 66-70. We first create an instance of shared pointer holding the solver (line 66, TNL::Solvers::getLinearSolver) and the same with the preconditioner (line 67, TNL::Solvers::getPreconditioner). The rest of the code is the same as in the previous examples with the only difference that we work with the pointer solver_ptr instead of the direct instance solver of the solver type.

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]