Template Numerical Library version\ main:7388705
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:57
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: 31345 RES: 0.0058578
Jacobi stage: ITER: 62597 RES: 4.8191e-05
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: 683 RES: 0.80624
Jacobi stage: ITER: 1375 RES: 0.67138
Jacobi stage: ITER: 2063 RES: 0.57447
Jacobi stage: ITER: 2729 RES: 0.5023
Jacobi stage: ITER: 3415 RES: 0.4429
Jacobi stage: ITER: 4103 RES: 0.39326
Jacobi stage: ITER: 4779 RES: 0.35159
Jacobi stage: ITER: 5447 RES: 0.31569
Jacobi stage: ITER: 6127 RES: 0.28345
Jacobi stage: ITER: 6607 RES: 0.2629
Jacobi stage: ITER: 7061 RES: 0.24486
Jacobi stage: ITER: 7521 RES: 0.22798
Jacobi stage: ITER: 7977 RES: 0.21244
Jacobi stage: ITER: 8415 RES: 0.1986
Jacobi stage: ITER: 8857 RES: 0.18545
Jacobi stage: ITER: 9195 RES: 0.17609
Jacobi stage: ITER: 9471 RES: 0.16876
Jacobi stage: ITER: 9747 RES: 0.16174
Jacobi stage: ITER: 10031 RES: 0.15483
Jacobi stage: ITER: 10311 RES: 0.1483
Jacobi stage: ITER: 10585 RES: 0.14214
Jacobi stage: ITER: 10863 RES: 0.13623
Jacobi stage: ITER: 11135 RES: 0.13065
Jacobi stage: ITER: 11545 RES: 0.12263
Jacobi stage: ITER: 11997 RES: 0.11441
Jacobi stage: ITER: 12449 RES: 0.10673
Jacobi stage: ITER: 12905 RES: 0.099508
Jacobi stage: ITER: 13365 RES: 0.092718
Jacobi stage: ITER: 13831 RES: 0.086339
Jacobi stage: ITER: 14295 RES: 0.080399
Jacobi stage: ITER: 14763 RES: 0.074822
Jacobi stage: ITER: 15219 RES: 0.069761
Jacobi stage: ITER: 15639 RES: 0.065402
Jacobi stage: ITER: 16069 RES: 0.061203
Jacobi stage: ITER: 16529 RES: 0.057028
Jacobi stage: ITER: 16991 RES: 0.053137
Jacobi stage: ITER: 17451 RES: 0.049512
Jacobi stage: ITER: 17919 RES: 0.046078
Jacobi stage: ITER: 18379 RES: 0.042935
Jacobi stage: ITER: 18851 RES: 0.039932
Jacobi stage: ITER: 19317 RES: 0.037162
Jacobi stage: ITER: 19757 RES: 0.034734
Jacobi stage: ITER: 20187 RES: 0.032524
Jacobi stage: ITER: 20645 RES: 0.030305
Jacobi stage: ITER: 21101 RES: 0.028255
Jacobi stage: ITER: 21567 RES: 0.026311
Jacobi stage: ITER: 22031 RES: 0.024501
Jacobi stage: ITER: 22497 RES: 0.022802
Jacobi stage: ITER: 22949 RES: 0.021272
Jacobi stage: ITER: 23405 RES: 0.019833
Jacobi stage: ITER: 23855 RES: 0.018515
Jacobi stage: ITER: 24299 RES: 0.017294
Jacobi stage: ITER: 24749 RES: 0.016134
Jacobi stage: ITER: 25209 RES: 0.015033
Jacobi stage: ITER: 25669 RES: 0.014008
Jacobi stage: ITER: 26127 RES: 0.01306
Jacobi stage: ITER: 26571 RES: 0.012199
Jacobi stage: ITER: 26999 RES: 0.011423
Jacobi stage: ITER: 27463 RES: 0.010637
Jacobi stage: ITER: 27927 RES: 0.0099055
Jacobi stage: ITER: 28383 RES: 0.0092354
Jacobi stage: ITER: 28849 RES: 0.0085948
Jacobi stage: ITER: 29319 RES: 0.0079987
Jacobi stage: ITER: 29783 RES: 0.0074484
Jacobi stage: ITER: 30229 RES: 0.0069531
Jacobi stage: ITER: 30639 RES: 0.0065307
Jacobi stage: ITER: 31089 RES: 0.0060927
Jacobi stage: ITER: 31537 RES: 0.0056876
Jacobi stage: ITER: 31985 RES: 0.0053093
Jacobi stage: ITER: 32443 RES: 0.0049502
Jacobi stage: ITER: 32903 RES: 0.0046125
Jacobi stage: ITER: 33371 RES: 0.0042926
Jacobi stage: ITER: 33831 RES: 0.0039997
Jacobi stage: ITER: 34291 RES: 0.0037269
Jacobi stage: ITER: 34733 RES: 0.0034812
Jacobi stage: ITER: 35175 RES: 0.0032537
Jacobi stage: ITER: 35635 RES: 0.0030317
Jacobi stage: ITER: 36091 RES: 0.0028266
Jacobi stage: ITER: 36555 RES: 0.0026322
Jacobi stage: ITER: 37023 RES: 0.0024496
Jacobi stage: ITER: 37493 RES: 0.0022783
Jacobi stage: ITER: 37957 RES: 0.0021216
Jacobi stage: ITER: 38403 RES: 0.0019817
Jacobi stage: ITER: 38827 RES: 0.0018568
Jacobi stage: ITER: 39269 RES: 0.0017343
Jacobi stage: ITER: 39719 RES: 0.001619
Jacobi stage: ITER: 40169 RES: 0.0015104
Jacobi stage: ITER: 40629 RES: 0.0014074
Jacobi stage: ITER: 41091 RES: 0.0013114
Jacobi stage: ITER: 41539 RES: 0.0012242
Jacobi stage: ITER: 41987 RES: 0.0011428
Jacobi stage: ITER: 42427 RES: 0.0010681
Jacobi stage: ITER: 42865 RES: 0.00099828
Jacobi stage: ITER: 43303 RES: 0.00093362
Jacobi stage: ITER: 43743 RES: 0.0008726
Jacobi stage: ITER: 44181 RES: 0.00081558
Jacobi stage: ITER: 44619 RES: 0.00076275
Jacobi stage: ITER: 45059 RES: 0.0007129
Jacobi stage: ITER: 45497 RES: 0.00066631
Jacobi stage: ITER: 45937 RES: 0.00062277
Jacobi stage: ITER: 46375 RES: 0.00058243
Jacobi stage: ITER: 46815 RES: 0.00054436
Jacobi stage: ITER: 47255 RES: 0.00050879
Jacobi stage: ITER: 47693 RES: 0.00047554
Jacobi stage: ITER: 48133 RES: 0.00044446
Jacobi stage: ITER: 48571 RES: 0.00041567
Jacobi stage: ITER: 49013 RES: 0.00038827
Jacobi stage: ITER: 49451 RES: 0.00036312
Jacobi stage: ITER: 49891 RES: 0.00033939
Jacobi stage: ITER: 50331 RES: 0.00031721
Jacobi stage: ITER: 50769 RES: 0.00029648
Jacobi stage: ITER: 51207 RES: 0.00027727
Jacobi stage: ITER: 51647 RES: 0.00025915
Jacobi stage: ITER: 52085 RES: 0.00024222
Jacobi stage: ITER: 52523 RES: 0.00022653
Jacobi stage: ITER: 52961 RES: 0.00021172
Jacobi stage: ITER: 53399 RES: 0.00019801
Jacobi stage: ITER: 53837 RES: 0.00018507
Jacobi stage: ITER: 54275 RES: 0.00017308
Jacobi stage: ITER: 54715 RES: 0.00016177
Jacobi stage: ITER: 55151 RES: 0.00015129
Jacobi stage: ITER: 55591 RES: 0.0001414
Jacobi stage: ITER: 56029 RES: 0.00013216
Jacobi stage: ITER: 56467 RES: 0.0001236
Jacobi stage: ITER: 56905 RES: 0.00011552
Jacobi stage: ITER: 57343 RES: 0.00010804
Jacobi stage: ITER: 57781 RES: 0.00010098
Jacobi stage: ITER: 58219 RES: 9.4438e-05
Jacobi stage: ITER: 58657 RES: 8.8266e-05
Jacobi stage: ITER: 59101 RES: 8.2447e-05
Jacobi stage: ITER: 59539 RES: 7.7107e-05
Jacobi stage: ITER: 59979 RES: 7.2068e-05
Jacobi stage: ITER: 60621 RES: 6.528e-05
Jacobi stage: ITER: 61059 RES: 6.1051e-05
Jacobi stage: ITER: 61497 RES: 5.7062e-05
Jacobi stage: ITER: 61935 RES: 5.3365e-05
Jacobi stage: ITER: 62373 RES: 4.9878e-05
Jacobi stage: ITER: 62811 RES: 4.6647e-05
Jacobi stage: ITER: 63249 RES: 4.3598e-05
Jacobi stage: ITER: 63687 RES: 4.0774e-05
Jacobi stage: ITER: 64123 RES: 3.8133e-05
Jacobi stage: ITER: 64563 RES: 3.5641e-05
Jacobi stage: ITER: 64999 RES: 3.3332e-05
Jacobi stage: ITER: 65437 RES: 3.1154e-05
Jacobi stage: ITER: 65875 RES: 2.9136e-05
Jacobi stage: ITER: 66313 RES: 2.7232e-05
Jacobi stage: ITER: 66751 RES: 2.5468e-05
Jacobi stage: ITER: 67187 RES: 2.3818e-05
Jacobi stage: ITER: 67625 RES: 2.2262e-05
Jacobi stage: ITER: 68063 RES: 2.0819e-05
Jacobi stage: ITER: 68501 RES: 1.9459e-05
Jacobi stage: ITER: 68939 RES: 1.8198e-05
Jacobi stage: ITER: 69375 RES: 1.702e-05
Jacobi stage: ITER: 69813 RES: 1.5907e-05
Jacobi stage: ITER: 70251 RES: 1.4877e-05
Jacobi stage: ITER: 70687 RES: 1.3913e-05
Jacobi stage: ITER: 71125 RES: 1.3004e-05
Jacobi stage: ITER: 71563 RES: 1.2162e-05
Jacobi stage: ITER: 71999 RES: 1.1374e-05
Jacobi stage: ITER: 72437 RES: 1.0631e-05
Jacobi stage: ITER: 72875 RES: 9.9419e-06
Jacobi stage: ITER: 73311 RES: 9.2979e-06
Jacobi stage: ITER: 73749 RES: 8.6903e-06
Jacobi stage: ITER: 74185 RES: 8.1273e-06
Jacobi stage: ITER: 74623 RES: 7.6009e-06
Jacobi stage: ITER: 75059 RES: 7.1085e-06
Jacobi stage: ITER: 75497 RES: 6.644e-06
Jacobi stage: ITER: 75935 RES: 6.2136e-06
Jacobi stage: ITER: 76371 RES: 5.8111e-06
Jacobi stage: ITER: 76809 RES: 5.4313e-06
Jacobi stage: ITER: 77247 RES: 5.0795e-06
Jacobi stage: ITER: 77683 RES: 4.7505e-06
Jacobi stage: ITER: 78121 RES: 4.44e-06
Jacobi stage: ITER: 78559 RES: 4.1524e-06
Jacobi stage: ITER: 78995 RES: 3.8834e-06
Jacobi stage: ITER: 79433 RES: 3.6296e-06
Jacobi stage: ITER: 79871 RES: 3.3945e-06
Jacobi stage: ITER: 80307 RES: 3.1746e-06
Jacobi stage: ITER: 80745 RES: 2.9672e-06
Jacobi stage: ITER: 81183 RES: 2.775e-06
Jacobi stage: ITER: 81619 RES: 2.5952e-06
Jacobi stage: ITER: 82057 RES: 2.4256e-06
Jacobi stage: ITER: 82495 RES: 2.2685e-06
Jacobi stage: ITER: 82931 RES: 2.1215e-06
Jacobi stage: ITER: 83369 RES: 1.9829e-06
Jacobi stage: ITER: 83807 RES: 1.8544e-06
Jacobi stage: ITER: 84243 RES: 1.7343e-06
Jacobi stage: ITER: 84723 RES: 1.611e-06
Jacobi stage: ITER: 85405 RES: 1.4504e-06
Jacobi stage: ITER: 86065 RES: 1.3105e-06
Jacobi stage: ITER: 86727 RES: 1.1842e-06
Jacobi stage: ITER: 87387 RES: 1.07e-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:3.3851e-
ELA: 0.10018 Jacobi stage: ITER: 30165 RES: 0.0070218
ELA: 0.20032 Jacobi stage: ITER: 63495 RES: 4.1969e-05
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.216e-0
ELA: 0.10016 Jacobi stage: ITER: 3119 RES: 0.46716
ELA: 0.20043 Jacobi stage: ITER: 3821 RES: 0.41245
ELA: 0.30072 Jacobi stage: ITER: 4519 RES: 0.36691
ELA: 0.40087 Jacobi stage: ITER: 5219 RES: 0.32743
ELA: 0.50103 Jacobi stage: ITER: 5897 RES: 0.29383
ELA: 0.60122 Jacobi stage: ITER: 6587 RES: 0.26373
ELA: 0.70134 Jacobi stage: ITER: 7289 RES: 0.23633
ELA: 0.80151 Jacobi stage: ITER: 7991 RES: 0.21204
ELA: 0.90167 Jacobi stage: ITER: 8671 RES: 0.19091
ELA: 1.0018 Jacobi stage: ITER: 9347 RES: 0.17202
ELA: 1.1019 Jacobi stage: ITER: 9811 RES: 0.16016
ELA: 1.202 Jacobi stage: ITER: 10263 RES: 0.1494
ELA: 1.3022 Jacobi stage: ITER: 10721 RES: 0.1392
ELA: 1.4025 Jacobi stage: ITER: 11177 RES: 0.12977
ELA: 1.5026 Jacobi stage: ITER: 11615 RES: 0.12136
ELA: 1.6034 Jacobi stage: ITER: 12065 RES: 0.11322
ELA: 1.7036 Jacobi stage: ITER: 12367 RES: 0.10812
ELA: 1.8037 Jacobi stage: ITER: 12647 RES: 0.10356
ELA: 1.9038 Jacobi stage: ITER: 12927 RES: 0.099202
ELA: 2.004 Jacobi stage: ITER: 13207 RES: 0.095026
ELA: 2.1041 Jacobi stage: ITER: 13489 RES: 0.090969
ELA: 2.2043 Jacobi stage: ITER: 13763 RES: 0.087246
ELA: 2.3044 Jacobi stage: ITER: 14039 RES: 0.083624
ELA: 2.4046 Jacobi stage: ITER: 14315 RES: 0.080153
ELA: 2.5047 Jacobi stage: ITER: 14755 RES: 0.074914
ELA: 2.6049 Jacobi stage: ITER: 15209 RES: 0.069846
ELA: 2.705 Jacobi stage: ITER: 15663 RES: 0.065161
ELA: 2.8051 Jacobi stage: ITER: 16117 RES: 0.060753
ELA: 2.9053 Jacobi stage: ITER: 16579 RES: 0.056609
ELA: 3.0054 Jacobi stage: ITER: 17047 RES: 0.052682
ELA: 3.1056 Jacobi stage: ITER: 17515 RES: 0.049028
ELA: 3.2057 Jacobi stage: ITER: 17979 RES: 0.045655
ELA: 3.3059 Jacobi stage: ITER: 18439 RES: 0.042541
ELA: 3.406 Jacobi stage: ITER: 18851 RES: 0.039932
ELA: 3.5062 Jacobi stage: ITER: 19289 RES: 0.037322
ELA: 3.6063 Jacobi stage: ITER: 19757 RES: 0.034734
ELA: 3.7065 Jacobi stage: ITER: 20221 RES: 0.032344
ELA: 3.8066 Jacobi stage: ITER: 20687 RES: 0.030119
ELA: 3.9067 Jacobi stage: ITER: 21151 RES: 0.028047
ELA: 4.0068 Jacobi stage: ITER: 21613 RES: 0.026118
ELA: 4.107 Jacobi stage: ITER: 22085 RES: 0.024291
ELA: 4.2071 Jacobi stage: ITER: 22551 RES: 0.02262
ELA: 4.3074 Jacobi stage: ITER: 22989 RES: 0.021142
ELA: 4.4076 Jacobi stage: ITER: 23423 RES: 0.019785
ELA: 4.5077 Jacobi stage: ITER: 23883 RES: 0.018435
ELA: 4.6078 Jacobi stage: ITER: 24345 RES: 0.017167
ELA: 4.708 Jacobi stage: ITER: 24815 RES: 0.015976
ELA: 4.8081 Jacobi stage: ITER: 25277 RES: 0.014877
ELA: 4.9082 Jacobi stage: ITER: 25739 RES: 0.013862
ELA: 5.0084 Jacobi stage: ITER: 26191 RES: 0.012932
ELA: 5.1092 Jacobi stage: ITER: 26647 RES: 0.012058
ELA: 5.2093 Jacobi stage: ITER: 27105 RES: 0.011235
ELA: 5.3095 Jacobi stage: ITER: 27555 RES: 0.010488
ELA: 5.4096 Jacobi stage: ITER: 28009 RES: 0.0097785
ELA: 5.5097 Jacobi stage: ITER: 28467 RES: 0.009117
ELA: 5.6099 Jacobi stage: ITER: 28931 RES: 0.0084899
ELA: 5.7103 Jacobi stage: ITER: 29385 RES: 0.0079156
ELA: 5.8104 Jacobi stage: ITER: 29827 RES: 0.0073983
ELA: 5.9105 Jacobi stage: ITER: 30257 RES: 0.0069233
ELA: 6.0107 Jacobi stage: ITER: 30725 RES: 0.0064431
ELA: 6.1114 Jacobi stage: ITER: 31189 RES: 0.0059998
ELA: 6.2116 Jacobi stage: ITER: 31651 RES: 0.0055905
ELA: 6.3117 Jacobi stage: ITER: 32119 RES: 0.0052028
ELA: 6.4118 Jacobi stage: ITER: 32589 RES: 0.0048389
ELA: 6.512 Jacobi stage: ITER: 33051 RES: 0.0045088
ELA: 6.6121 Jacobi stage: ITER: 33495 RES: 0.0042116
ELA: 6.7122 Jacobi stage: ITER: 33909 RES: 0.0039509
ELA: 6.8124 Jacobi stage: ITER: 34363 RES: 0.0036859
ELA: 6.9126 Jacobi stage: ITER: 34811 RES: 0.0034408
ELA: 7.0134 Jacobi stage: ITER: 35267 RES: 0.003208
ELA: 7.1136 Jacobi stage: ITER: 35725 RES: 0.0029892
ELA: 7.2137 Jacobi stage: ITER: 36187 RES: 0.0027853
ELA: 7.3138 Jacobi stage: ITER: 36655 RES: 0.0025921
ELA: 7.414 Jacobi stage: ITER: 37115 RES: 0.0024152
ELA: 7.5141 Jacobi stage: ITER: 37575 RES: 0.0022505
ELA: 7.6142 Jacobi stage: ITER: 38015 RES: 0.0021034
ELA: 7.7143 Jacobi stage: ITER: 38455 RES: 0.0019659
ELA: 7.8144 Jacobi stage: ITER: 38921 RES: 0.0018296
ELA: 7.9146 Jacobi stage: ITER: 39379 RES: 0.0017058
ELA: 8.0147 Jacobi stage: ITER: 39843 RES: 0.0015885
ELA: 8.1153 Jacobi stage: ITER: 40307 RES: 0.0014792
ELA: 8.2154 Jacobi stage: ITER: 40777 RES: 0.0013758
ELA: 8.3155 Jacobi stage: ITER: 41239 RES: 0.0012819
ELA: 8.4157 Jacobi stage: ITER: 41679 RES: 0.0011981
ELA: 8.5164 Jacobi stage: ITER: 42109 RES: 0.0011212
ELA: 8.6166 Jacobi stage: ITER: 42557 RES: 0.0010466
ELA: 8.7167 Jacobi stage: ITER: 43011 RES: 0.00097644
ELA: 8.8168 Jacobi stage: ITER: 43463 RES: 0.00091095
ELA: 8.917 Jacobi stage: ITER: 43919 RES: 0.00084933
ELA: 9.0171 Jacobi stage: ITER: 44389 RES: 0.00078993
ELA: 9.1172 Jacobi stage: ITER: 44827 RES: 0.00073876
ELA: 9.2174 Jacobi stage: ITER: 45279 RES: 0.00068921
ELA: 9.3175 Jacobi stage: ITER: 45719 RES: 0.00064417
ELA: 9.4177 Jacobi stage: ITER: 46157 RES: 0.00060207
ELA: 9.5178 Jacobi stage: ITER: 46597 RES: 0.00056273
ELA: 9.618 Jacobi stage: ITER: 47035 RES: 0.00052628
ELA: 9.7182 Jacobi stage: ITER: 47475 RES: 0.00049188
ELA: 9.8183 Jacobi stage: ITER: 47913 RES: 0.00045974
ELA: 9.9185 Jacobi stage: ITER: 48351 RES: 0.00042996
ELA: 10.019 Jacobi stage: ITER: 48791 RES: 0.00040186
ELA: 10.119 Jacobi stage: ITER: 49231 RES: 0.0003756
ELA: 10.219 Jacobi stage: ITER: 49669 RES: 0.00035105
ELA: 10.319 Jacobi stage: ITER: 50109 RES: 0.00032811
ELA: 10.419 Jacobi stage: ITER: 50547 RES: 0.00030686
ELA: 10.519 Jacobi stage: ITER: 50987 RES: 0.0002868
ELA: 10.62 Jacobi stage: ITER: 51427 RES: 0.00026806
ELA: 10.72 Jacobi stage: ITER: 51867 RES: 0.00025054
ELA: 10.82 Jacobi stage: ITER: 52307 RES: 0.00023417
ELA: 10.92 Jacobi stage: ITER: 52745 RES: 0.00021886
ELA: 11.02 Jacobi stage: ITER: 53185 RES: 0.00020456
ELA: 11.12 Jacobi stage: ITER: 53623 RES: 0.00019131
ELA: 11.22 Jacobi stage: ITER: 54063 RES: 0.00017881
ELA: 11.321 Jacobi stage: ITER: 54503 RES: 0.00016712
ELA: 11.421 Jacobi stage: ITER: 54941 RES: 0.0001562
ELA: 11.521 Jacobi stage: ITER: 55379 RES: 0.00014608
ELA: 11.621 Jacobi stage: ITER: 55817 RES: 0.00013654
ELA: 11.721 Jacobi stage: ITER: 56255 RES: 0.00012769
ELA: 11.821 Jacobi stage: ITER: 56695 RES: 0.00011935
ELA: 11.922 Jacobi stage: ITER: 57133 RES: 0.00011155
ELA: 12.022 Jacobi stage: ITER: 57571 RES: 0.00010432
ELA: 12.122 Jacobi stage: ITER: 58011 RES: 9.7504e-05
ELA: 12.222 Jacobi stage: ITER: 58449 RES: 9.1132e-05
ELA: 12.322 Jacobi stage: ITER: 58887 RES: 8.5229e-05
ELA: 12.422 Jacobi stage: ITER: 59327 RES: 7.9659e-05
ELA: 12.523 Jacobi stage: ITER: 59763 RES: 7.4499e-05
ELA: 12.623 Jacobi stage: ITER: 60203 RES: 6.963e-05
ELA: 12.723 Jacobi stage: ITER: 60641 RES: 6.508e-05
ELA: 12.823 Jacobi stage: ITER: 61079 RES: 6.0864e-05
ELA: 12.923 Jacobi stage: ITER: 61517 RES: 5.6887e-05
ELA: 13.023 Jacobi stage: ITER: 61955 RES: 5.3202e-05
ELA: 13.124 Jacobi stage: ITER: 62397 RES: 4.9694e-05
ELA: 13.224 Jacobi stage: ITER: 62835 RES: 4.6475e-05
ELA: 13.324 Jacobi stage: ITER: 63273 RES: 4.3438e-05
ELA: 13.456 Jacobi stage: ITER: 63851 RES: 3.976e-05
ELA: 13.556 Jacobi stage: ITER: 64289 RES: 3.7162e-05
ELA: 13.656 Jacobi stage: ITER: 64727 RES: 3.4754e-05
ELA: 13.757 Jacobi stage: ITER: 65165 RES: 3.2483e-05
ELA: 13.857 Jacobi stage: ITER: 65603 RES: 3.0379e-05
ELA: 13.957 Jacobi stage: ITER: 66041 RES: 2.8394e-05
ELA: 14.057 Jacobi stage: ITER: 66479 RES: 2.6554e-05
ELA: 14.157 Jacobi stage: ITER: 66917 RES: 2.4819e-05
ELA: 14.257 Jacobi stage: ITER: 67355 RES: 2.3211e-05
ELA: 14.357 Jacobi stage: ITER: 67791 RES: 2.1708e-05
ELA: 14.458 Jacobi stage: ITER: 68231 RES: 2.0289e-05
ELA: 14.558 Jacobi stage: ITER: 68667 RES: 1.8975e-05
ELA: 14.658 Jacobi stage: ITER: 69107 RES: 1.7735e-05
ELA: 14.758 Jacobi stage: ITER: 69543 RES: 1.6586e-05
ELA: 14.858 Jacobi stage: ITER: 69983 RES: 1.5502e-05
ELA: 14.958 Jacobi stage: ITER: 70419 RES: 1.4498e-05
ELA: 15.059 Jacobi stage: ITER: 70857 RES: 1.355e-05
ELA: 15.159 Jacobi stage: ITER: 71295 RES: 1.2673e-05
ELA: 15.259 Jacobi stage: ITER: 71731 RES: 1.1852e-05
ELA: 15.359 Jacobi stage: ITER: 72169 RES: 1.1077e-05
ELA: 15.459 Jacobi stage: ITER: 72607 RES: 1.036e-05
ELA: 15.559 Jacobi stage: ITER: 73043 RES: 9.6886e-06
ELA: 15.66 Jacobi stage: ITER: 73481 RES: 9.0555e-06
ELA: 15.76 Jacobi stage: ITER: 73919 RES: 8.4689e-06
ELA: 15.86 Jacobi stage: ITER: 74355 RES: 7.9203e-06
ELA: 15.96 Jacobi stage: ITER: 74793 RES: 7.4027e-06
ELA: 16.06 Jacobi stage: ITER: 75231 RES: 6.9232e-06
ELA: 16.16 Jacobi stage: ITER: 75667 RES: 6.4747e-06
ELA: 16.26 Jacobi stage: ITER: 76105 RES: 6.0516e-06
ELA: 16.361 Jacobi stage: ITER: 76543 RES: 5.6596e-06
ELA: 16.461 Jacobi stage: ITER: 76979 RES: 5.293e-06
ELA: 16.561 Jacobi stage: ITER: 77417 RES: 4.9471e-06
ELA: 16.661 Jacobi stage: ITER: 77855 RES: 4.6266e-06
ELA: 16.761 Jacobi stage: ITER: 78291 RES: 4.3269e-06
ELA: 16.861 Jacobi stage: ITER: 78729 RES: 4.0441e-06
ELA: 16.962 Jacobi stage: ITER: 79167 RES: 3.7822e-06
ELA: 17.062 Jacobi stage: ITER: 79603 RES: 3.5372e-06
ELA: 17.162 Jacobi stage: ITER: 80041 RES: 3.306e-06
ELA: 17.262 Jacobi stage: ITER: 80479 RES: 3.0919e-06
ELA: 17.362 Jacobi stage: ITER: 80915 RES: 2.8916e-06
ELA: 17.462 Jacobi stage: ITER: 81353 RES: 2.7026e-06
ELA: 17.563 Jacobi stage: ITER: 81791 RES: 2.5275e-06
ELA: 17.663 Jacobi stage: ITER: 82227 RES: 2.3638e-06
ELA: 17.763 Jacobi stage: ITER: 82665 RES: 2.2093e-06
ELA: 17.863 Jacobi stage: ITER: 83103 RES: 2.0662e-06
ELA: 17.963 Jacobi stage: ITER: 83539 RES: 1.9324e-06
ELA: 18.064 Jacobi stage: ITER: 83977 RES: 1.8061e-06
ELA: 18.164 Jacobi stage: ITER: 84415 RES: 1.6891e-06
ELA: 18.264 Jacobi stage: ITER: 84851 RES: 1.5797e-06
ELA: 18.364 Jacobi stage: ITER: 85289 RES: 1.4764e-06
ELA: 18.464 Jacobi stage: ITER: 85727 RES: 1.3808e-06
ELA: 18.564 Jacobi stage: ITER: 86163 RES: 1.2914e-06
ELA: 18.664 Jacobi stage: ITER: 86601 RES: 1.207e-06
ELA: 18.765 Jacobi stage: ITER: 87039 RES: 1.1288e-06
ELA: 18.865 Jacobi stage: ITER: 87477 RES: 1.055e-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 ]