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

Table of Contents

Introduction

Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the following 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/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Host.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/Linear/TFQMR.h>
7
8template< typename Device >
9void
10iterativeLinearSolverExample()
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 */
21 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
22 using Vector = TNL::Containers::Vector< double, Device >;
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 {
30 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
37 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
38 }
39 else {
40 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
41 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
42 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
43 }
44 };
45
46 /***
47 * Set the matrix elements.
48 */
49 matrix_ptr->forAllRows( f );
50 std::cout << *matrix_ptr << '\n';
51
52 /***
53 * Set the right-hand side vector.
54 */
55 Vector x( size, 1.0 );
56 Vector b( size );
57 matrix_ptr->vectorProduct( x, b );
58 x = 0.0;
59 std::cout << "Vector b = " << b << '\n';
60
61 /***
62 * Solve the linear system.
63 */
64 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
65 LinearSolver solver;
66 solver.setMatrix( matrix_ptr );
67 solver.setConvergenceResidue( 1.0e-6 );
68 solver.solve( b, x );
69 std::cout << "Vector x = " << x << '\n';
70}
71
72int
73main( int argc, char* argv[] )
74{
75 std::cout << "Solving linear system on host:\n";
76 iterativeLinearSolverExample< TNL::Devices::Sequential >();
77
78#ifdef __CUDACC__
79 std::cout << "Solving linear system on CUDA device:\n";
80 iterativeLinearSolverExample< TNL::Devices::Cuda >();
81#endif
82}
__cuda_callable__
#define __cuda_callable__
This macro serves for annotating functions which are supposed to be called even from the GPU device.
Definition Macros.h:50
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition Vector.h:37
TNL::Matrices::SparseMatrix
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:59
TNL::Solvers::Linear::LinearSolver::setMatrix
virtual void setMatrix(const MatrixPointer &matrix)
Set the matrix of the linear system.
Definition LinearSolver.h:120
TNL::Solvers::Linear::TFQMR
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition TFQMR.h:21
std::make_shared
T make_shared(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 elements of the matrix \( A \) are set using the method TNL::Matrices::SparseMatrix::forAllRows. 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 \) and multiply it with matrix \( A \) to get the right-hand side vector \( \vec b \). Finally, we reset the vector \( \vec x \) to zero.

To solve the linear system in the example, we use the TFQMR solver. 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 and set the matrix of the linear system. Note that the matrix is passed to the solver as a std::shared_ptr. Then we set the stopping criterion for the iterative method in terms of the relative residue, i.e. \( \lVert \vec b - A \vec x \rVert / \lVert b \rVert \). The solver is executed by calling the TNL::Solvers::Linear::LinearSolver::solve method which 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 or the solver status in general. For this purpose, TNL provides a solver monitor which can show various metrics in real time, such as current number of iterations, current residue of the approximate solution, etc. The solver monitor in TNL runs in a separate thread and it refreshes the status of the solver with a configurable refresh rate (once per 500 ms by default). The use of the solver monitor is demonstrated in the following example.

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

First, we set up 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. Then, we set the relaxation parameter \( \omega \) of the Jacobi solver to 0.0005. The reason is to artificially slow down the convergence, because we want to see some iterations in this example. Next, we create an instance of the solver monitor and a special thread for the monitor (an instance of the TNL::Solvers::SolverMonitorThread class). We use the following methods to configure the solver monitor:

Next, we call TNL::Solvers::IterativeSolver::setSolverMonitor to connect the solver with the monitor and we set the convergence criterion based on the relative residue. Finally, we start the solver using the TNL::Solvers::Linear::Jacobi::solve method and when the solver finishes, we stop the monitor using 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: 25679 RES: 0.013991
Jacobi stage: ITER: 51597 RES: 0.00026107
Jacobi stage: ITER: 77657 RES: 4.768e-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: 4757 RES: 0.35274
Jacobi stage: ITER: 8663 RES: 0.19114
Jacobi stage: ITER: 9131 RES: 0.17784
Jacobi stage: ITER: 9559 RES: 0.16649
Jacobi stage: ITER: 9867 RES: 0.15878
Jacobi stage: ITER: 10175 RES: 0.15143
Jacobi stage: ITER: 10481 RES: 0.14443
Jacobi stage: ITER: 10789 RES: 0.13775
Jacobi stage: ITER: 11095 RES: 0.13146
Jacobi stage: ITER: 11407 RES: 0.1253
Jacobi stage: ITER: 11715 RES: 0.11951
Jacobi stage: ITER: 12023 RES: 0.11398
Jacobi stage: ITER: 12327 RES: 0.10878
Jacobi stage: ITER: 12653 RES: 0.10344
Jacobi stage: ITER: 12987 RES: 0.098292
Jacobi stage: ITER: 13357 RES: 0.092832
Jacobi stage: ITER: 13853 RES: 0.086021
Jacobi stage: ITER: 14099 RES: 0.082857
Jacobi stage: ITER: 14309 RES: 0.080202
Jacobi stage: ITER: 14479 RES: 0.078159
Jacobi stage: ITER: 14701 RES: 0.075515
Jacobi stage: ITER: 14963 RES: 0.072558
Jacobi stage: ITER: 15199 RES: 0.069975
Jacobi stage: ITER: 15323 RES: 0.068655
Jacobi stage: ITER: 15535 RES: 0.066455
Jacobi stage: ITER: 15769 RES: 0.064089
Jacobi stage: ITER: 15965 RES: 0.062189
Jacobi stage: ITER: 16187 RES: 0.060122
Jacobi stage: ITER: 16307 RES: 0.059024
Jacobi stage: ITER: 16559 RES: 0.056783
Jacobi stage: ITER: 16783 RES: 0.054862
Jacobi stage: ITER: 17045 RES: 0.052682
Jacobi stage: ITER: 17215 RES: 0.05134
Jacobi stage: ITER: 17463 RES: 0.049421
Jacobi stage: ITER: 17775 RES: 0.047109
Jacobi stage: ITER: 18087 RES: 0.044904
Jacobi stage: ITER: 18399 RES: 0.042803
Jacobi stage: ITER: 18711 RES: 0.0408
Jacobi stage: ITER: 19023 RES: 0.038891
Jacobi stage: ITER: 19323 RES: 0.037139
Jacobi stage: ITER: 19641 RES: 0.035358
Jacobi stage: ITER: 19953 RES: 0.033704
Jacobi stage: ITER: 20263 RES: 0.032146
Jacobi stage: ITER: 20579 RES: 0.030623
Jacobi stage: ITER: 20873 RES: 0.029262
Jacobi stage: ITER: 21185 RES: 0.027893
Jacobi stage: ITER: 21505 RES: 0.026555
Jacobi stage: ITER: 21811 RES: 0.025343
Jacobi stage: ITER: 22127 RES: 0.024143
Jacobi stage: ITER: 22443 RES: 0.022999
Jacobi stage: ITER: 22745 RES: 0.02195
Jacobi stage: ITER: 23041 RES: 0.020974
Jacobi stage: ITER: 23351 RES: 0.020005
Jacobi stage: ITER: 23663 RES: 0.019069
Jacobi stage: ITER: 23973 RES: 0.018176
Jacobi stage: ITER: 24281 RES: 0.017336
Jacobi stage: ITER: 24599 RES: 0.016515
Jacobi stage: ITER: 24905 RES: 0.015752
Jacobi stage: ITER: 25221 RES: 0.015006
Jacobi stage: ITER: 25523 RES: 0.01433
Jacobi stage: ITER: 25823 RES: 0.013685
Jacobi stage: ITER: 26135 RES: 0.013044
Jacobi stage: ITER: 26449 RES: 0.012426
Jacobi stage: ITER: 26759 RES: 0.011852
Jacobi stage: ITER: 27075 RES: 0.01129
Jacobi stage: ITER: 27375 RES: 0.010782
Jacobi stage: ITER: 27681 RES: 0.010284
Jacobi stage: ITER: 27985 RES: 0.0098146
Jacobi stage: ITER: 28299 RES: 0.0093554
Jacobi stage: ITER: 28607 RES: 0.0089231
Jacobi stage: ITER: 28917 RES: 0.0085055
Jacobi stage: ITER: 29225 RES: 0.0081125
Jacobi stage: ITER: 29543 RES: 0.0077281
Jacobi stage: ITER: 29851 RES: 0.007371
Jacobi stage: ITER: 30147 RES: 0.0070434
Jacobi stage: ITER: 30447 RES: 0.0067262
Jacobi stage: ITER: 30763 RES: 0.0064075
Jacobi stage: ITER: 31075 RES: 0.0061077
Jacobi stage: ITER: 31389 RES: 0.0058183
Jacobi stage: ITER: 31701 RES: 0.0055461
Jacobi stage: ITER: 32003 RES: 0.0052963
Jacobi stage: ITER: 32313 RES: 0.0050485
Jacobi stage: ITER: 32615 RES: 0.0048211
Jacobi stage: ITER: 32927 RES: 0.0045955
Jacobi stage: ITER: 33239 RES: 0.0043805
Jacobi stage: ITER: 33557 RES: 0.0041704
Jacobi stage: ITER: 33875 RES: 0.0039728
Jacobi stage: ITER: 34187 RES: 0.0037869
Jacobi stage: ITER: 34465 RES: 0.0036275
Jacobi stage: ITER: 34755 RES: 0.0034705
Jacobi stage: ITER: 35077 RES: 0.003302
Jacobi stage: ITER: 35383 RES: 0.0031514
Jacobi stage: ITER: 35691 RES: 0.0030057
Jacobi stage: ITER: 35875 RES: 0.002922
Jacobi stage: ITER: 36219 RES: 0.0027716
Jacobi stage: ITER: 36715 RES: 0.0025683
Jacobi stage: ITER: 37141 RES: 0.0024049
Jacobi stage: ITER: 37513 RES: 0.0022713
Jacobi stage: ITER: 37915 RES: 0.002136
Jacobi stage: ITER: 38203 RES: 0.0020435
Jacobi stage: ITER: 38987 RES: 0.0018117
Jacobi stage: ITER: 39295 RES: 0.001728
Jacobi stage: ITER: 40487 RES: 0.0014389
Jacobi stage: ITER: 43135 RES: 0.00095802
Jacobi stage: ITER: 43965 RES: 0.00084309
Jacobi stage: ITER: 44567 RES: 0.00076886
Jacobi stage: ITER: 45029 RES: 0.00071597
Jacobi stage: ITER: 45461 RES: 0.00067001
Jacobi stage: ITER: 45973 RES: 0.00061933
Jacobi stage: ITER: 46599 RES: 0.00056273
Jacobi stage: ITER: 47051 RES: 0.00052498
Jacobi stage: ITER: 47663 RES: 0.00047788
Jacobi stage: ITER: 47893 RES: 0.00046115
Jacobi stage: ITER: 48061 RES: 0.0004494
Jacobi stage: ITER: 48191 RES: 0.00044065
Jacobi stage: ITER: 48331 RES: 0.00043128
Jacobi stage: ITER: 48491 RES: 0.00042081
Jacobi stage: ITER: 48795 RES: 0.00040161
Jacobi stage: ITER: 49213 RES: 0.00037652
Jacobi stage: ITER: 49343 RES: 0.00036919
Jacobi stage: ITER: 49517 RES: 0.00035934
Jacobi stage: ITER: 49737 RES: 0.0003474
Jacobi stage: ITER: 49915 RES: 0.00033814
Jacobi stage: ITER: 50041 RES: 0.00033155
Jacobi stage: ITER: 50303 RES: 0.00031857
Jacobi stage: ITER: 50579 RES: 0.00030535
Jacobi stage: ITER: 51027 RES: 0.00028504
Jacobi stage: ITER: 51661 RES: 0.00025852
Jacobi stage: ITER: 52287 RES: 0.00023489
Jacobi stage: ITER: 52651 RES: 0.00022212
Jacobi stage: ITER: 52945 RES: 0.00021224
Jacobi stage: ITER: 53267 RES: 0.00020206
Jacobi stage: ITER: 53747 RES: 0.0001877
Jacobi stage: ITER: 54159 RES: 0.00017619
Jacobi stage: ITER: 54343 RES: 0.00017128
Jacobi stage: ITER: 54565 RES: 0.00016549
Jacobi stage: ITER: 54945 RES: 0.0001561
Jacobi stage: ITER: 55311 RES: 0.00014762
Jacobi stage: ITER: 55731 RES: 0.00013839
Jacobi stage: ITER: 56085 RES: 0.00013103
Jacobi stage: ITER: 56491 RES: 0.00012315
Jacobi stage: ITER: 56731 RES: 0.00011869
Jacobi stage: ITER: 57215 RES: 0.00011018
Jacobi stage: ITER: 57823 RES: 0.00010036
Jacobi stage: ITER: 58153 RES: 9.5371e-05
Jacobi stage: ITER: 58683 RES: 8.7942e-05
Jacobi stage: ITER: 58879 RES: 8.5333e-05
Jacobi stage: ITER: 59339 RES: 7.9512e-05
Jacobi stage: ITER: 59815 RES: 7.3906e-05
Jacobi stage: ITER: 60947 RES: 6.2111e-05
Jacobi stage: ITER: 61263 RES: 5.9168e-05
Jacobi stage: ITER: 61587 RES: 5.6295e-05
Jacobi stage: ITER: 61947 RES: 5.3267e-05
Jacobi stage: ITER: 62251 RES: 5.0837e-05
Jacobi stage: ITER: 62459 RES: 4.9238e-05
Jacobi stage: ITER: 62691 RES: 4.7515e-05
Jacobi stage: ITER: 62961 RES: 4.557e-05
Jacobi stage: ITER: 63187 RES: 4.4029e-05
Jacobi stage: ITER: 63511 RES: 4.1892e-05
Jacobi stage: ITER: 63707 RES: 4.0649e-05
Jacobi stage: ITER: 63933 RES: 3.925e-05
Jacobi stage: ITER: 64179 RES: 3.7806e-05
Jacobi stage: ITER: 64427 RES: 3.6393e-05
Jacobi stage: ITER: 64819 RES: 3.4267e-05
Jacobi stage: ITER: 65175 RES: 3.2443e-05
Jacobi stage: ITER: 65843 RES: 2.9279e-05
Jacobi stage: ITER: 67155 RES: 2.3935e-05
Jacobi stage: ITER: 68917 RES: 1.8254e-05
Jacobi stage: ITER: 70147 RES: 1.5116e-05
Jacobi stage: ITER: 72079 RES: 1.1235e-05
Jacobi stage: ITER: 73607 RES: 8.8846e-06
Jacobi stage: ITER: 74105 RES: 8.2278e-06
Jacobi stage: ITER: 75139 RES: 7.0217e-06
Jacobi stage: ITER: 76195 RES: 5.9703e-06
Jacobi stage: ITER: 77023 RES: 5.2573e-06
Jacobi stage: ITER: 77395 RES: 4.9653e-06
Jacobi stage: ITER: 77639 RES: 4.7827e-06
Jacobi stage: ITER: 77863 RES: 4.6209e-06
Jacobi stage: ITER: 78083 RES: 4.4674e-06
Jacobi stage: ITER: 78373 RES: 4.2714e-06
Jacobi stage: ITER: 79427 RES: 3.6341e-06
Jacobi stage: ITER: 80691 RES: 2.9928e-06
Jacobi stage: ITER: 81069 RES: 2.8231e-06
Jacobi stage: ITER: 81463 RES: 2.6581e-06
Jacobi stage: ITER: 81831 RES: 2.5121e-06
Jacobi stage: ITER: 82031 RES: 2.4361e-06
Jacobi stage: ITER: 82307 RES: 2.3349e-06
Jacobi stage: ITER: 82673 RES: 2.2066e-06
Jacobi stage: ITER: 83039 RES: 2.0866e-06
Jacobi stage: ITER: 83381 RES: 1.9792e-06
Jacobi stage: ITER: 83751 RES: 1.8705e-06
Jacobi stage: ITER: 84127 RES: 1.7655e-06
Jacobi stage: ITER: 84499 RES: 1.6674e-06
Jacobi stage: ITER: 84835 RES: 1.5836e-06
Jacobi stage: ITER: 85165 RES: 1.5048e-06
Jacobi stage: ITER: 85515 RES: 1.4265e-06
Jacobi stage: ITER: 85739 RES: 1.3783e-06
Jacobi stage: ITER: 85957 RES: 1.3325e-06
Jacobi stage: ITER: 86181 RES: 1.2874e-06
Jacobi stage: ITER: 86387 RES: 1.2477e-06
Jacobi stage: ITER: 86575 RES: 1.2122e-06
Jacobi stage: ITER: 86887 RES: 1.1555e-06
Jacobi stage: ITER: 87211 RES: 1.0994e-06
Jacobi stage: ITER: 87519 RES: 1.0486e-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 <TNL/Timer.h>
4#include <TNL/Matrices/SparseMatrix.h>
5#include <TNL/Devices/Sequential.h>
6#include <TNL/Devices/Cuda.h>
7#include <TNL/Solvers/Linear/Jacobi.h>
8
9template< typename Device >
10void
11iterativeLinearSolverExample()
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 */
22 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
23 using Vector = TNL::Containers::Vector< double, Device >;
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 {
31 const int rowIdx = row.getRowIndex();
32 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else {
41 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
42 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
43 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
44 }
45 };
46
47 /***
48 * Set the matrix elements.
49 */
50 matrix_ptr->forAllRows( f );
51 std::cout << *matrix_ptr << '\n';
52
53 /***
54 * Set the right-hand side vector.
55 */
56 Vector x( size, 1.0 );
57 Vector b( size );
58 matrix_ptr->vectorProduct( x, b );
59 x = 0.0;
60 std::cout << "Vector b = " << b << '\n';
61
62 /***
63 * Setup solver of the linear system.
64 */
65 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
66 LinearSolver solver;
67 solver.setMatrix( matrix_ptr );
68 solver.setOmega( 0.0005 );
69
70 /***
71 * Setup monitor of the iterative solver.
72 */
73 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double >;
74 IterativeSolverMonitorType monitor;
75 TNL::Solvers::SolverMonitorThread mmonitorThread( monitor );
76 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
77 monitor.setVerbose( 1 );
78 monitor.setStage( "Jacobi stage:" );
79 TNL::Timer timer;
80 monitor.setTimer( timer );
81 timer.start();
82 solver.setSolverMonitor( monitor );
83 solver.setConvergenceResidue( 1.0e-6 );
84 solver.solve( b, x );
85 monitor.stopMainLoop();
86 std::cout << "Vector x = " << x << '\n';
87}
88
89int
90main( int argc, char* argv[] )
91{
92 std::cout << "Solving linear system on host:\n";
93 iterativeLinearSolverExample< TNL::Devices::Sequential >();
94
95#ifdef __CUDACC__
96 std::cout << "Solving linear system on CUDA device:\n";
97 iterativeLinearSolverExample< TNL::Devices::Cuda >();
98#endif
99}
TNL::Timer
Class for real time, CPU time and CPU cycles measuring.
Definition Timer.h:25
TNL::Timer::start
void start()
Starts timer.

The only changes are around the lines where we create an instance of TNL::Timer, connect it with the monitor using TNL::Solvers::SolverMonitor::setTimer and start the timer with 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:2.7485e-
ELA: 0.10015 Jacobi stage: ITER: 26717 RES: 0.011925
ELA: 0.20068 Jacobi stage: ITER: 55277 RES: 0.00014834
ELA: 0.30084 Jacobi stage: ITER: 85139 RES: 1.5113e-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:2.9244e-
ELA: 0.10014 Jacobi stage: ITER: 467 RES: 0.85961
ELA: 0.20026 Jacobi stage: ITER: 903 RES: 0.75805
ELA: 0.30037 Jacobi stage: ITER: 1211 RES: 0.69924
ELA: 0.40113 Jacobi stage: ITER: 1519 RES: 0.64863
ELA: 0.50125 Jacobi stage: ITER: 1827 RES: 0.60459
ELA: 0.60137 Jacobi stage: ITER: 2133 RES: 0.56583
ELA: 0.70149 Jacobi stage: ITER: 2441 RES: 0.53135
ELA: 0.80167 Jacobi stage: ITER: 2739 RES: 0.50154
ELA: 0.90178 Jacobi stage: ITER: 3045 RES: 0.47337
ELA: 1.0019 Jacobi stage: ITER: 3355 RES: 0.44766
ELA: 1.102 Jacobi stage: ITER: 3669 RES: 0.42343
ELA: 1.2031 Jacobi stage: ITER: 3999 RES: 0.40024
ELA: 1.3033 Jacobi stage: ITER: 4247 RES: 0.38386
ELA: 1.4071 Jacobi stage: ITER: 4483 RES: 0.36909
ELA: 1.5072 Jacobi stage: ITER: 4955 RES: 0.34167
ELA: 1.6096 Jacobi stage: ITER: 5227 RES: 0.32701
ELA: 1.7097 Jacobi stage: ITER: 5507 RES: 0.31268
ELA: 1.8099 Jacobi stage: ITER: 5793 RES: 0.2987
ELA: 1.91 Jacobi stage: ITER: 6035 RES: 0.28759
ELA: 2.0101 Jacobi stage: ITER: 6199 RES: 0.28026
ELA: 2.1102 Jacobi stage: ITER: 6413 RES: 0.27092
ELA: 2.2103 Jacobi stage: ITER: 6701 RES: 0.25899
ELA: 2.3104 Jacobi stage: ITER: 6803 RES: 0.25498
ELA: 2.4105 Jacobi stage: ITER: 7021 RES: 0.24639
ELA: 2.5106 Jacobi stage: ITER: 7251 RES: 0.23781
ELA: 2.6107 Jacobi stage: ITER: 7369 RES: 0.23342
ELA: 2.7109 Jacobi stage: ITER: 7581 RES: 0.22587
ELA: 2.811 Jacobi stage: ITER: 7825 RES: 0.21749
ELA: 2.9111 Jacobi stage: ITER: 7963 RES: 0.21296
ELA: 3.0112 Jacobi stage: ITER: 8177 RES: 0.20597
ELA: 3.1113 Jacobi stage: ITER: 8397 RES: 0.19909
ELA: 3.2114 Jacobi stage: ITER: 8681 RES: 0.19055
ELA: 3.3115 Jacobi stage: ITER: 8991 RES: 0.18171
ELA: 3.4116 Jacobi stage: ITER: 9301 RES: 0.17319
ELA: 3.5117 Jacobi stage: ITER: 9613 RES: 0.16506
ELA: 3.6119 Jacobi stage: ITER: 9927 RES: 0.15732
ELA: 3.712 Jacobi stage: ITER: 10233 RES: 0.15004
ELA: 3.8121 Jacobi stage: ITER: 10533 RES: 0.14328
ELA: 3.9122 Jacobi stage: ITER: 10847 RES: 0.13657
ELA: 4.0123 Jacobi stage: ITER: 11161 RES: 0.13009
ELA: 4.1145 Jacobi stage: ITER: 11477 RES: 0.12392
ELA: 4.2147 Jacobi stage: ITER: 11791 RES: 0.11812
ELA: 4.3153 Jacobi stage: ITER: 12087 RES: 0.11287
ELA: 4.4154 Jacobi stage: ITER: 12395 RES: 0.10765
ELA: 4.5156 Jacobi stage: ITER: 12703 RES: 0.10268
ELA: 4.6157 Jacobi stage: ITER: 13009 RES: 0.09793
ELA: 4.7158 Jacobi stage: ITER: 13323 RES: 0.093347
ELA: 4.8159 Jacobi stage: ITER: 13639 RES: 0.088924
ELA: 4.9161 Jacobi stage: ITER: 13947 RES: 0.084814
ELA: 5.0162 Jacobi stage: ITER: 14235 RES: 0.081144
ELA: 5.1164 Jacobi stage: ITER: 14543 RES: 0.077394
ELA: 5.2165 Jacobi stage: ITER: 14851 RES: 0.073818
ELA: 5.3166 Jacobi stage: ITER: 15155 RES: 0.07045
ELA: 5.4171 Jacobi stage: ITER: 15465 RES: 0.067153
ELA: 5.5173 Jacobi stage: ITER: 15781 RES: 0.063971
ELA: 5.6174 Jacobi stage: ITER: 16089 RES: 0.061015
ELA: 5.718 Jacobi stage: ITER: 16401 RES: 0.05816
ELA: 5.8202 Jacobi stage: ITER: 16711 RES: 0.055473
ELA: 5.9203 Jacobi stage: ITER: 17009 RES: 0.052974
ELA: 6.0205 Jacobi stage: ITER: 17317 RES: 0.050526
ELA: 6.1206 Jacobi stage: ITER: 17627 RES: 0.048192
ELA: 6.2207 Jacobi stage: ITER: 17935 RES: 0.045965
ELA: 6.3211 Jacobi stage: ITER: 18251 RES: 0.043787
ELA: 6.4212 Jacobi stage: ITER: 18551 RES: 0.041815
ELA: 6.5213 Jacobi stage: ITER: 18857 RES: 0.039883
ELA: 6.6214 Jacobi stage: ITER: 19159 RES: 0.038087
ELA: 6.7216 Jacobi stage: ITER: 19469 RES: 0.036305
ELA: 6.8241 Jacobi stage: ITER: 19785 RES: 0.034585
ELA: 6.9242 Jacobi stage: ITER: 20095 RES: 0.032986
ELA: 7.0251 Jacobi stage: ITER: 20403 RES: 0.031462
ELA: 7.1252 Jacobi stage: ITER: 20719 RES: 0.029972
ELA: 7.2254 Jacobi stage: ITER: 21035 RES: 0.028552
ELA: 7.3261 Jacobi stage: ITER: 21325 RES: 0.027299
ELA: 7.4262 Jacobi stage: ITER: 21623 RES: 0.026086
ELA: 7.5264 Jacobi stage: ITER: 21935 RES: 0.024865
ELA: 7.6265 Jacobi stage: ITER: 22243 RES: 0.023716
ELA: 7.7266 Jacobi stage: ITER: 22559 RES: 0.022593
ELA: 7.8267 Jacobi stage: ITER: 22867 RES: 0.021549
ELA: 7.9268 Jacobi stage: ITER: 23171 RES: 0.020566
ELA: 8.0269 Jacobi stage: ITER: 23479 RES: 0.019615
ELA: 8.1271 Jacobi stage: ITER: 23779 RES: 0.018732
ELA: 8.2271 Jacobi stage: ITER: 24087 RES: 0.017866
ELA: 8.3273 Jacobi stage: ITER: 24401 RES: 0.01702
ELA: 8.4274 Jacobi stage: ITER: 24717 RES: 0.016213
ELA: 8.5285 Jacobi stage: ITER: 25039 RES: 0.015436
ELA: 8.6286 Jacobi stage: ITER: 25351 RES: 0.014714
ELA: 8.7288 Jacobi stage: ITER: 25631 RES: 0.014094
ELA: 8.8289 Jacobi stage: ITER: 25919 RES: 0.013484
ELA: 8.929 Jacobi stage: ITER: 26215 RES: 0.012885
ELA: 9.0291 Jacobi stage: ITER: 26519 RES: 0.012297
ELA: 9.1299 Jacobi stage: ITER: 28077 RES: 0.0096769
ELA: 9.23 Jacobi stage: ITER: 30215 RES: 0.0069702
ELA: 9.3417 Jacobi stage: ITER: 31195 RES: 0.0059962
ELA: 9.4418 Jacobi stage: ITER: 32145 RES: 0.0051804
ELA: 9.544 Jacobi stage: ITER: 32975 RES: 0.0045618
ELA: 9.6441 Jacobi stage: ITER: 34785 RES: 0.0034535
ELA: 9.7465 Jacobi stage: ITER: 36299 RES: 0.0027377
ELA: 9.8466 Jacobi stage: ITER: 37213 RES: 0.0023784
ELA: 9.9467 Jacobi stage: ITER: 39059 RES: 0.0017918
ELA: 10.047 Jacobi stage: ITER: 40153 RES: 0.0015141
ELA: 10.15 Jacobi stage: ITER: 41115 RES: 0.0013066
ELA: 10.25 Jacobi stage: ITER: 41513 RES: 0.0012287
ELA: 10.356 Jacobi stage: ITER: 42107 RES: 0.0011219
ELA: 10.457 Jacobi stage: ITER: 43035 RES: 0.00097285
ELA: 10.557 Jacobi stage: ITER: 43455 RES: 0.00091207
ELA: 10.657 Jacobi stage: ITER: 43927 RES: 0.00084829
ELA: 10.757 Jacobi stage: ITER: 44403 RES: 0.00078848
ELA: 10.858 Jacobi stage: ITER: 44723 RES: 0.00075066
ELA: 10.959 Jacobi stage: ITER: 45403 RES: 0.00067621
ELA: 11.059 Jacobi stage: ITER: 46007 RES: 0.0006163
ELA: 11.159 Jacobi stage: ITER: 46625 RES: 0.00056031
ELA: 11.259 Jacobi stage: ITER: 47675 RES: 0.000477
ELA: 11.36 Jacobi stage: ITER: 48863 RES: 0.00039744
ELA: 11.46 Jacobi stage: ITER: 49479 RES: 0.00036156
ELA: 11.56 Jacobi stage: ITER: 49829 RES: 0.00034253
ELA: 11.66 Jacobi stage: ITER: 50909 RES: 0.00029017
ELA: 11.76 Jacobi stage: ITER: 51511 RES: 0.00026462
ELA: 11.86 Jacobi stage: ITER: 51867 RES: 0.00025054
ELA: 11.96 Jacobi stage: ITER: 52131 RES: 0.00024058
ELA: 12.06 Jacobi stage: ITER: 52511 RES: 0.00022694
ELA: 12.161 Jacobi stage: ITER: 52751 RES: 0.00021873
ELA: 12.263 Jacobi stage: ITER: 52971 RES: 0.00021146
ELA: 12.364 Jacobi stage: ITER: 53383 RES: 0.00019849
ELA: 12.464 Jacobi stage: ITER: 53629 RES: 0.00019108
ELA: 12.564 Jacobi stage: ITER: 54323 RES: 0.00017181
ELA: 12.665 Jacobi stage: ITER: 55473 RES: 0.00014394
ELA: 12.765 Jacobi stage: ITER: 56143 RES: 0.00012991
ELA: 12.865 Jacobi stage: ITER: 56369 RES: 0.00012544
ELA: 12.965 Jacobi stage: ITER: 56561 RES: 0.00012179
ELA: 13.066 Jacobi stage: ITER: 56893 RES: 0.00011574
ELA: 13.166 Jacobi stage: ITER: 57191 RES: 0.00011059
ELA: 13.266 Jacobi stage: ITER: 58103 RES: 9.6136e-05
ELA: 13.366 Jacobi stage: ITER: 58815 RES: 8.6176e-05
ELA: 13.466 Jacobi stage: ITER: 59049 RES: 8.3109e-05
ELA: 13.568 Jacobi stage: ITER: 59213 RES: 8.1041e-05
ELA: 13.668 Jacobi stage: ITER: 59431 RES: 7.8396e-05
ELA: 13.769 Jacobi stage: ITER: 60455 RES: 6.6986e-05
ELA: 13.87 Jacobi stage: ITER: 60809 RES: 6.3422e-05
ELA: 13.976 Jacobi stage: ITER: 62673 RES: 4.7632e-05
ELA: 14.076 Jacobi stage: ITER: 63725 RES: 4.0524e-05
ELA: 14.176 Jacobi stage: ITER: 64485 RES: 3.6059e-05
ELA: 14.28 Jacobi stage: ITER: 65019 RES: 3.323e-05
ELA: 14.38 Jacobi stage: ITER: 65489 RES: 3.0906e-05
ELA: 14.48 Jacobi stage: ITER: 65995 RES: 2.8604e-05
ELA: 14.581 Jacobi stage: ITER: 66539 RES: 2.6311e-05
ELA: 14.681 Jacobi stage: ITER: 66797 RES: 2.5281e-05
ELA: 14.785 Jacobi stage: ITER: 67097 RES: 2.4142e-05
ELA: 14.886 Jacobi stage: ITER: 67571 RES: 2.2454e-05
ELA: 14.998 Jacobi stage: ITER: 67879 RES: 2.1416e-05
ELA: 15.098 Jacobi stage: ITER: 68191 RES: 2.0414e-05
ELA: 15.199 Jacobi stage: ITER: 68407 RES: 1.9748e-05
ELA: 15.299 Jacobi stage: ITER: 68567 RES: 1.9269e-05
ELA: 15.399 Jacobi stage: ITER: 68711 RES: 1.8847e-05
ELA: 15.5 Jacobi stage: ITER: 68949 RES: 1.8165e-05
ELA: 15.608 Jacobi stage: ITER: 69219 RES: 1.7432e-05
ELA: 15.711 Jacobi stage: ITER: 69427 RES: 1.6884e-05
ELA: 15.813 Jacobi stage: ITER: 69705 RES: 1.6173e-05
ELA: 15.914 Jacobi stage: ITER: 70367 RES: 1.4614e-05
ELA: 16.014 Jacobi stage: ITER: 70735 RES: 1.3811e-05
ELA: 16.114 Jacobi stage: ITER: 71009 RES: 1.3238e-05
ELA: 16.214 Jacobi stage: ITER: 71347 RES: 1.2572e-05
ELA: 16.314 Jacobi stage: ITER: 71897 RES: 1.155e-05
ELA: 16.414 Jacobi stage: ITER: 72257 RES: 1.0929e-05
ELA: 16.515 Jacobi stage: ITER: 72647 RES: 1.0296e-05
ELA: 16.615 Jacobi stage: ITER: 72847 RES: 9.9848e-06
ELA: 16.716 Jacobi stage: ITER: 73065 RES: 9.653e-06
ELA: 16.817 Jacobi stage: ITER: 73271 RES: 9.3552e-06
ELA: 16.919 Jacobi stage: ITER: 73541 RES: 8.9724e-06
ELA: 17.02 Jacobi stage: ITER: 73743 RES: 8.7009e-06
ELA: 17.126 Jacobi stage: ITER: 74003 RES: 8.3603e-06
ELA: 17.226 Jacobi stage: ITER: 74383 RES: 7.8863e-06
ELA: 17.326 Jacobi stage: ITER: 74607 RES: 7.6196e-06
ELA: 17.427 Jacobi stage: ITER: 74979 RES: 7.1964e-06
ELA: 17.529 Jacobi stage: ITER: 75151 RES: 7.0088e-06
ELA: 17.629 Jacobi stage: ITER: 75371 RES: 6.7759e-06
ELA: 17.729 Jacobi stage: ITER: 75595 RES: 6.5467e-06
ELA: 17.829 Jacobi stage: ITER: 75839 RES: 6.3059e-06
ELA: 17.929 Jacobi stage: ITER: 76081 RES: 6.0739e-06
ELA: 18.029 Jacobi stage: ITER: 76303 RES: 5.8721e-06
ELA: 18.129 Jacobi stage: ITER: 76537 RES: 5.663e-06
ELA: 18.229 Jacobi stage: ITER: 76815 RES: 5.428e-06
ELA: 18.336 Jacobi stage: ITER: 77575 RES: 4.8299e-06
ELA: 18.436 Jacobi stage: ITER: 77961 RES: 4.5505e-06
ELA: 18.536 Jacobi stage: ITER: 79383 RES: 3.6587e-06
ELA: 18.636 Jacobi stage: ITER: 80985 RES: 2.8598e-06
ELA: 18.737 Jacobi stage: ITER: 81347 RES: 2.7059e-06
ELA: 18.837 Jacobi stage: ITER: 81711 RES: 2.5588e-06
ELA: 18.937 Jacobi stage: ITER: 82269 RES: 2.3479e-06
ELA: 19.037 Jacobi stage: ITER: 82627 RES: 2.2229e-06
ELA: 19.137 Jacobi stage: ITER: 82991 RES: 2.1021e-06
ELA: 19.237 Jacobi stage: ITER: 83241 RES: 2.0223e-06
ELA: 19.338 Jacobi stage: ITER: 83445 RES: 1.9599e-06
ELA: 19.438 Jacobi stage: ITER: 83743 RES: 1.8728e-06
ELA: 19.538 Jacobi stage: ITER: 83975 RES: 1.8072e-06
ELA: 19.638 Jacobi stage: ITER: 84201 RES: 1.745e-06
ELA: 19.738 Jacobi stage: ITER: 84563 RES: 1.6511e-06
ELA: 19.838 Jacobi stage: ITER: 84835 RES: 1.5836e-06
ELA: 19.938 Jacobi stage: ITER: 85171 RES: 1.5039e-06
ELA: 20.039 Jacobi stage: ITER: 85523 RES: 1.4248e-06
ELA: 20.146 Jacobi stage: ITER: 85905 RES: 1.3432e-06
ELA: 20.247 Jacobi stage: ITER: 86211 RES: 1.2819e-06
ELA: 20.347 Jacobi stage: ITER: 86499 RES: 1.2264e-06
ELA: 20.447 Jacobi stage: ITER: 87443 RES: 1.0609e-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/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Host.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/Linear/TFQMR.h>
7#include <TNL/Solvers/Linear/Preconditioners/Diagonal.h>
8
9template< typename Device >
10void
11iterativeLinearSolverExample()
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 */
22 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
23 using Vector = TNL::Containers::Vector< double, Device >;
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 {
31 const int rowIdx = row.getRowIndex();
32 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else {
41 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
42 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
43 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
44 }
45 };
46
47 /***
48 * Set the matrix elements.
49 */
50 matrix_ptr->forAllRows( f );
51 std::cout << *matrix_ptr << '\n';
52
53 /***
54 * Set the right-hand side vector.
55 */
56 Vector x( size, 1.0 );
57 Vector b( size );
58 matrix_ptr->vectorProduct( x, b );
59 x = 0.0;
60 std::cout << "Vector b = " << b << '\n';
61
62 /***
63 * Solve the linear system using diagonal (Jacobi) preconditioner.
64 */
65 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
66 using Preconditioner = TNL::Solvers::Linear::Preconditioners::Diagonal< MatrixType >;
67 auto preconditioner_ptr = std::make_shared< Preconditioner >();
68 preconditioner_ptr->update( matrix_ptr );
69 LinearSolver solver;
70 solver.setMatrix( matrix_ptr );
71 solver.setPreconditioner( preconditioner_ptr );
72 solver.setConvergenceResidue( 1.0e-6 );
73 solver.solve( b, x );
74 std::cout << "Vector x = " << x << '\n';
75}
76
77int
78main( int argc, char* argv[] )
79{
80 std::cout << "Solving linear system on host:\n";
81 iterativeLinearSolverExample< TNL::Devices::Sequential >();
82
83#ifdef __CUDACC__
84 std::cout << "Solving linear system on CUDA device:\n";
85 iterativeLinearSolverExample< TNL::Devices::Cuda >();
86#endif
87}
TNL::Solvers::Linear::Preconditioners::Diagonal
Diagonal (Jacobi) preconditioner for iterative solvers of linear systems.
Definition Diagonal.h:21

In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen in the solver setup. Similarly to the matrix of the linear system, the preconditioner needs to be passed to the solver as a std::shared_ptr. When the preconditioner object is created, we have to initialize it using the update method, which 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 using the setPreconditioner method.

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 provides the functions TNL::Solvers::getLinearSolver and TNL::Solvers::getPreconditioner for selecting the linear solver and preconditioner at runtime. The following example shows how to use these functions:

1#include <iostream>
2#include <memory>
3#include <TNL/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Host.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/LinearSolverTypeResolver.h>
7
8template< typename Device >
9void
10iterativeLinearSolverExample()
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 */
21 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
22 using Vector = TNL::Containers::Vector< double, Device >;
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 {
30 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
37 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
38 }
39 else {
40 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
41 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
42 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
43 }
44 };
45
46 /***
47 * Set the matrix elements.
48 */
49 matrix_ptr->forAllRows( f );
50 std::cout << *matrix_ptr << '\n';
51
52 /***
53 * Set the right-hand side vector.
54 */
55 Vector x( size, 1.0 );
56 Vector b( size );
57 matrix_ptr->vectorProduct( x, b );
58 x = 0.0;
59 std::cout << "Vector b = " << b << '\n';
60
61 /***
62 * Solve the linear system using diagonal (Jacobi) preconditioner.
63 */
64 auto solver_ptr = TNL::Solvers::getLinearSolver< MatrixType >( "tfqmr" );
65 auto preconditioner_ptr = TNL::Solvers::getPreconditioner< MatrixType >( "diagonal" );
66 preconditioner_ptr->update( matrix_ptr );
67 solver_ptr->setMatrix( matrix_ptr );
68 solver_ptr->setPreconditioner( preconditioner_ptr );
69 solver_ptr->setConvergenceResidue( 1.0e-6 );
70 solver_ptr->solve( b, x );
71 std::cout << "Vector x = " << x << '\n';
72}
73
74int
75main( int argc, char* argv[] )
76{
77 std::cout << "Solving linear system on host:\n";
78 iterativeLinearSolverExample< TNL::Devices::Sequential >();
79
80#ifdef __CUDACC__
81 std::cout << "Solving linear system on CUDA device:\n";
82 iterativeLinearSolverExample< TNL::Devices::Cuda >();
83#endif
84}
TNL::Solvers::getPreconditioner
std::shared_ptr< Linear::Preconditioners::Preconditioner< MatrixType > > getPreconditioner(const std::string &name)
Function returning shared pointer with linear preconditioner given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:144
TNL::Solvers::getLinearSolver
std::shared_ptr< Linear::LinearSolver< MatrixType > > getLinearSolver(const std::string &name)
Function returning shared pointer with linear solver given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:87

We still stay with the same problem and the only changes are in the solver setup. We first use TNL::Solvers::getLinearSolver to get a shared pointer holding the solver and then TNL::Solvers::getPreconditioner to get a shared pointer holding the preconditioner. 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 ]