Template Numerical Library version\ main:4e58ea6
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 << std::endl;
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 << std::endl;
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 << std::endl;
70}
71
72int
73main( int argc, char* argv[] )
74{
75 std::cout << "Solving linear system on host: " << std::endl;
76 iterativeLinearSolverExample< TNL::Devices::Sequential >();
77
78#ifdef __CUDACC__
79 std::cout << "Solving linear system on CUDA device: " << std::endl;
80 iterativeLinearSolverExample< TNL::Devices::Cuda >();
81#endif
82}
__cuda_callable__
#define __cuda_callable__
Definition Macros.h:49
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition Vector.h:36
TNL::Matrices::SparseMatrix
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:57
TNL::Solvers::Linear::LinearSolver::setMatrix
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::endl
T endl(T... args)
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 <chrono>
4#include <thread>
5#include <TNL/Matrices/SparseMatrix.h>
6#include <TNL/Devices/Sequential.h>
7#include <TNL/Devices/Cuda.h>
8#include <TNL/Solvers/Linear/Jacobi.h>
9
10template< typename Device >
11void
12iterativeLinearSolverExample()
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 */
23 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
24 using Vector = TNL::Containers::Vector< double, Device >;
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 {
32 const int rowIdx = row.getRowIndex();
33 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
39 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
40 }
41 else {
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 * Setup solver of the linear system.
65 */
66 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
67 LinearSolver solver;
68 solver.setMatrix( matrix_ptr );
69 solver.setOmega( 0.0005 );
70
71 /***
72 * Setup monitor of the iterative solver.
73 */
74 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
75 IterativeSolverMonitorType monitor;
76 TNL::Solvers::SolverMonitorThread monitorThread( monitor );
77 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
78 monitor.setVerbose( 1 );
79 monitor.setStage( "Jacobi stage:" );
80 solver.setSolverMonitor( monitor );
81 solver.setConvergenceResidue( 1.0e-6 );
82 solver.solve( b, x );
83 monitor.stopMainLoop();
84 std::cout << "Vector x = " << x << std::endl;
85}
86
87int
88main( int argc, char* argv[] )
89{
90 std::cout << "Solving linear system on host: " << std::endl;
91 iterativeLinearSolverExample< TNL::Devices::Sequential >();
92
93#ifdef __CUDACC__
94 std::cout << "Solving linear system on CUDA device: " << std::endl;
95 iterativeLinearSolverExample< TNL::Devices::Cuda >();
96#endif
97}
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:133

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: 27033 RES: 0.01136
Jacobi stage: ITER: 53807 RES: 0.00018586
Jacobi stage: ITER: 83467 RES: 1.9539e-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: 3995 RES: 0.40051
Jacobi stage: ITER: 4443 RES: 0.37154
Jacobi stage: ITER: 4769 RES: 0.35205
Jacobi stage: ITER: 5093 RES: 0.33403
Jacobi stage: ITER: 5415 RES: 0.31731
Jacobi stage: ITER: 5739 RES: 0.30136
Jacobi stage: ITER: 6063 RES: 0.28632
Jacobi stage: ITER: 6385 RES: 0.27211
Jacobi stage: ITER: 6709 RES: 0.25867
Jacobi stage: ITER: 7031 RES: 0.24608
Jacobi stage: ITER: 7355 RES: 0.234
Jacobi stage: ITER: 7679 RES: 0.22253
Jacobi stage: ITER: 8001 RES: 0.21165
Jacobi stage: ITER: 8323 RES: 0.20144
Jacobi stage: ITER: 8647 RES: 0.19161
Jacobi stage: ITER: 8971 RES: 0.18228
Jacobi stage: ITER: 9293 RES: 0.1734
Jacobi stage: ITER: 9617 RES: 0.16496
Jacobi stage: ITER: 9939 RES: 0.15703
Jacobi stage: ITER: 10263 RES: 0.1494
Jacobi stage: ITER: 10587 RES: 0.14214
Jacobi stage: ITER: 10911 RES: 0.13523
Jacobi stage: ITER: 11235 RES: 0.12866
Jacobi stage: ITER: 11539 RES: 0.12279
Jacobi stage: ITER: 11843 RES: 0.11718
Jacobi stage: ITER: 12171 RES: 0.11142
Jacobi stage: ITER: 12499 RES: 0.10595
Jacobi stage: ITER: 12827 RES: 0.10074
Jacobi stage: ITER: 13155 RES: 0.095788
Jacobi stage: ITER: 13449 RES: 0.091529
Jacobi stage: ITER: 13779 RES: 0.087032
Jacobi stage: ITER: 14109 RES: 0.082704
Jacobi stage: ITER: 14437 RES: 0.07864
Jacobi stage: ITER: 14765 RES: 0.074776
Jacobi stage: ITER: 15095 RES: 0.071102
Jacobi stage: ITER: 15411 RES: 0.067733
Jacobi stage: ITER: 15707 RES: 0.064722
Jacobi stage: ITER: 16011 RES: 0.06177
Jacobi stage: ITER: 16335 RES: 0.058771
Jacobi stage: ITER: 16663 RES: 0.055883
Jacobi stage: ITER: 16987 RES: 0.05317
Jacobi stage: ITER: 17311 RES: 0.050589
Jacobi stage: ITER: 17631 RES: 0.048162
Jacobi stage: ITER: 17959 RES: 0.045796
Jacobi stage: ITER: 18287 RES: 0.043546
Jacobi stage: ITER: 18617 RES: 0.041381
Jacobi stage: ITER: 18911 RES: 0.039566
Jacobi stage: ITER: 19227 RES: 0.037691
Jacobi stage: ITER: 19539 RES: 0.035927
Jacobi stage: ITER: 19867 RES: 0.034162
Jacobi stage: ITER: 20193 RES: 0.032484
Jacobi stage: ITER: 20523 RES: 0.030888
Jacobi stage: ITER: 20831 RES: 0.02946
Jacobi stage: ITER: 21155 RES: 0.02803
Jacobi stage: ITER: 21483 RES: 0.026653
Jacobi stage: ITER: 21813 RES: 0.025328
Jacobi stage: ITER: 22127 RES: 0.024143
Jacobi stage: ITER: 22427 RES: 0.023055
Jacobi stage: ITER: 22739 RES: 0.021977
Jacobi stage: ITER: 23063 RES: 0.02091
Jacobi stage: ITER: 23395 RES: 0.01987
Jacobi stage: ITER: 23727 RES: 0.018882
Jacobi stage: ITER: 24041 RES: 0.017988
Jacobi stage: ITER: 24363 RES: 0.017125
Jacobi stage: ITER: 24695 RES: 0.016273
Jacobi stage: ITER: 25023 RES: 0.015474
Jacobi stage: ITER: 25353 RES: 0.014704
Jacobi stage: ITER: 25617 RES: 0.01412
Jacobi stage: ITER: 25913 RES: 0.013493
Jacobi stage: ITER: 26239 RES: 0.012837
Jacobi stage: ITER: 26567 RES: 0.012207
Jacobi stage: ITER: 26897 RES: 0.0116
Jacobi stage: ITER: 27219 RES: 0.011043
Jacobi stage: ITER: 27543 RES: 0.010507
Jacobi stage: ITER: 27867 RES: 0.0099972
Jacobi stage: ITER: 28193 RES: 0.009506
Jacobi stage: ITER: 28523 RES: 0.0090389
Jacobi stage: ITER: 28841 RES: 0.0086054
Jacobi stage: ITER: 29139 RES: 0.0082229
Jacobi stage: ITER: 29463 RES: 0.0078237
Jacobi stage: ITER: 29793 RES: 0.0074347
Jacobi stage: ITER: 30123 RES: 0.0070694
Jacobi stage: ITER: 30453 RES: 0.006718
Jacobi stage: ITER: 30767 RES: 0.0064036
Jacobi stage: ITER: 31089 RES: 0.0060927
Jacobi stage: ITER: 31417 RES: 0.0057934
Jacobi stage: ITER: 31745 RES: 0.0055087
Jacobi stage: ITER: 32075 RES: 0.0052381
Jacobi stage: ITER: 32399 RES: 0.0049838
Jacobi stage: ITER: 32679 RES: 0.004774
Jacobi stage: ITER: 32979 RES: 0.004559
Jacobi stage: ITER: 33307 RES: 0.004335
Jacobi stage: ITER: 33635 RES: 0.004122
Jacobi stage: ITER: 33963 RES: 0.0039194
Jacobi stage: ITER: 34283 RES: 0.0037315
Jacobi stage: ITER: 34601 RES: 0.0035525
Jacobi stage: ITER: 34925 RES: 0.00338
Jacobi stage: ITER: 35255 RES: 0.0032139
Jacobi stage: ITER: 35571 RES: 0.0030617
Jacobi stage: ITER: 35887 RES: 0.0029166
Jacobi stage: ITER: 36191 RES: 0.0027835
Jacobi stage: ITER: 36519 RES: 0.0026468
Jacobi stage: ITER: 36847 RES: 0.0025167
Jacobi stage: ITER: 37171 RES: 0.0023946
Jacobi stage: ITER: 37499 RES: 0.0022769
Jacobi stage: ITER: 37825 RES: 0.002165
Jacobi stage: ITER: 38151 RES: 0.0020599
Jacobi stage: ITER: 38477 RES: 0.0019587
Jacobi stage: ITER: 38803 RES: 0.0018636
Jacobi stage: ITER: 39129 RES: 0.001772
Jacobi stage: ITER: 39455 RES: 0.001686
Jacobi stage: ITER: 39783 RES: 0.0016032
Jacobi stage: ITER: 40107 RES: 0.0015253
Jacobi stage: ITER: 40449 RES: 0.0014468
Jacobi stage: ITER: 40789 RES: 0.0013732
Jacobi stage: ITER: 41115 RES: 0.0013066
Jacobi stage: ITER: 41441 RES: 0.0012424
Jacobi stage: ITER: 41767 RES: 0.001182
Jacobi stage: ITER: 42095 RES: 0.001124
Jacobi stage: ITER: 42419 RES: 0.0010694
Jacobi stage: ITER: 42747 RES: 0.0010169
Jacobi stage: ITER: 43073 RES: 0.00096689
Jacobi stage: ITER: 43399 RES: 0.00091995
Jacobi stage: ITER: 43725 RES: 0.00087475
Jacobi stage: ITER: 44051 RES: 0.00083228
Jacobi stage: ITER: 44377 RES: 0.00079139
Jacobi stage: ITER: 44703 RES: 0.00075297
Jacobi stage: ITER: 45029 RES: 0.00071597
Jacobi stage: ITER: 45355 RES: 0.00068121
Jacobi stage: ITER: 45683 RES: 0.00064774
Jacobi stage: ITER: 46007 RES: 0.0006163
Jacobi stage: ITER: 46335 RES: 0.00058602
Jacobi stage: ITER: 46661 RES: 0.00055722
Jacobi stage: ITER: 46987 RES: 0.00053017
Jacobi stage: ITER: 47313 RES: 0.00050412
Jacobi stage: ITER: 47639 RES: 0.00047965
Jacobi stage: ITER: 47967 RES: 0.00045608
Jacobi stage: ITER: 48291 RES: 0.00043394
Jacobi stage: ITER: 48619 RES: 0.00041262
Jacobi stage: ITER: 48945 RES: 0.00039234
Jacobi stage: ITER: 49271 RES: 0.0003733
Jacobi stage: ITER: 49597 RES: 0.00035495
Jacobi stage: ITER: 49923 RES: 0.00033772
Jacobi stage: ITER: 50249 RES: 0.00032113
Jacobi stage: ITER: 50575 RES: 0.00030554
Jacobi stage: ITER: 50901 RES: 0.00029053
Jacobi stage: ITER: 51227 RES: 0.00027642
Jacobi stage: ITER: 51555 RES: 0.00026284
Jacobi stage: ITER: 51879 RES: 0.00025008
Jacobi stage: ITER: 52207 RES: 0.00023779
Jacobi stage: ITER: 52531 RES: 0.00022625
Jacobi stage: ITER: 52859 RES: 0.00021513
Jacobi stage: ITER: 53185 RES: 0.00020456
Jacobi stage: ITER: 53511 RES: 0.00019463
Jacobi stage: ITER: 53839 RES: 0.00018507
Jacobi stage: ITER: 54163 RES: 0.00017608
Jacobi stage: ITER: 54491 RES: 0.00016743
Jacobi stage: ITER: 54817 RES: 0.0001592
Jacobi stage: ITER: 55143 RES: 0.00015148
Jacobi stage: ITER: 55469 RES: 0.00014403
Jacobi stage: ITER: 55795 RES: 0.00013704
Jacobi stage: ITER: 56121 RES: 0.00013031
Jacobi stage: ITER: 56447 RES: 0.00012398
Jacobi stage: ITER: 56775 RES: 0.00011789
Jacobi stage: ITER: 57099 RES: 0.00011217
Jacobi stage: ITER: 57427 RES: 0.00010665
Jacobi stage: ITER: 57753 RES: 0.00010141
Jacobi stage: ITER: 58079 RES: 9.6491e-05
Jacobi stage: ITER: 58405 RES: 9.175e-05
Jacobi stage: ITER: 58731 RES: 8.7296e-05
Jacobi stage: ITER: 59059 RES: 8.3006e-05
Jacobi stage: ITER: 59383 RES: 7.8977e-05
Jacobi stage: ITER: 59711 RES: 7.5096e-05
Jacobi stage: ITER: 60037 RES: 7.1406e-05
Jacobi stage: ITER: 60363 RES: 6.794e-05
Jacobi stage: ITER: 60689 RES: 6.4602e-05
Jacobi stage: ITER: 61015 RES: 6.1465e-05
Jacobi stage: ITER: 61343 RES: 5.8445e-05
Jacobi stage: ITER: 61669 RES: 5.5574e-05
Jacobi stage: ITER: 61995 RES: 5.2876e-05
Jacobi stage: ITER: 62321 RES: 5.0278e-05
Jacobi stage: ITER: 62647 RES: 4.7837e-05
Jacobi stage: ITER: 62975 RES: 4.5486e-05
Jacobi stage: ITER: 63301 RES: 4.3252e-05
Jacobi stage: ITER: 63627 RES: 4.1152e-05
Jacobi stage: ITER: 63953 RES: 3.913e-05
Jacobi stage: ITER: 64279 RES: 3.723e-05
Jacobi stage: ITER: 64607 RES: 3.5401e-05
Jacobi stage: ITER: 64931 RES: 3.3682e-05
Jacobi stage: ITER: 65259 RES: 3.2027e-05
Jacobi stage: ITER: 65585 RES: 3.0454e-05
Jacobi stage: ITER: 65911 RES: 2.8975e-05
Jacobi stage: ITER: 66237 RES: 2.7552e-05
Jacobi stage: ITER: 66563 RES: 2.6214e-05
Jacobi stage: ITER: 66891 RES: 2.4926e-05
Jacobi stage: ITER: 67215 RES: 2.3716e-05
Jacobi stage: ITER: 67543 RES: 2.2551e-05
Jacobi stage: ITER: 67869 RES: 2.1443e-05
Jacobi stage: ITER: 68195 RES: 2.0402e-05
Jacobi stage: ITER: 68523 RES: 1.9399e-05
Jacobi stage: ITER: 68849 RES: 1.8446e-05
Jacobi stage: ITER: 69175 RES: 1.7551e-05
Jacobi stage: ITER: 69501 RES: 1.6688e-05
Jacobi stage: ITER: 69827 RES: 1.5878e-05
Jacobi stage: ITER: 70155 RES: 1.5098e-05
Jacobi stage: ITER: 70481 RES: 1.4356e-05
Jacobi stage: ITER: 70807 RES: 1.3659e-05
Jacobi stage: ITER: 71133 RES: 1.2988e-05
Jacobi stage: ITER: 71459 RES: 1.2357e-05
Jacobi stage: ITER: 71787 RES: 1.175e-05
Jacobi stage: ITER: 72111 RES: 1.118e-05
Jacobi stage: ITER: 72439 RES: 1.0631e-05
Jacobi stage: ITER: 72765 RES: 1.0108e-05
Jacobi stage: ITER: 73091 RES: 9.6175e-06
Jacobi stage: ITER: 73417 RES: 9.1449e-06
Jacobi stage: ITER: 73743 RES: 8.7009e-06
Jacobi stage: ITER: 74069 RES: 8.2734e-06
Jacobi stage: ITER: 74395 RES: 7.8718e-06
Jacobi stage: ITER: 74723 RES: 7.485e-06
Jacobi stage: ITER: 75047 RES: 7.1216e-06
Jacobi stage: ITER: 75375 RES: 6.7717e-06
Jacobi stage: ITER: 75701 RES: 6.439e-06
Jacobi stage: ITER: 76027 RES: 6.1264e-06
Jacobi stage: ITER: 76353 RES: 5.8254e-06
Jacobi stage: ITER: 76679 RES: 5.5426e-06
Jacobi stage: ITER: 77005 RES: 5.2702e-06
Jacobi stage: ITER: 77331 RES: 5.0144e-06
Jacobi stage: ITER: 77659 RES: 4.768e-06
Jacobi stage: ITER: 77983 RES: 4.5365e-06
Jacobi stage: ITER: 78311 RES: 4.3136e-06
Jacobi stage: ITER: 78637 RES: 4.1017e-06
Jacobi stage: ITER: 78963 RES: 3.9026e-06
Jacobi stage: ITER: 79291 RES: 3.7108e-06
Jacobi stage: ITER: 79615 RES: 3.5307e-06
Jacobi stage: ITER: 79943 RES: 3.3572e-06
Jacobi stage: ITER: 80269 RES: 3.1922e-06
Jacobi stage: ITER: 80595 RES: 3.0373e-06
Jacobi stage: ITER: 80921 RES: 2.888e-06
Jacobi stage: ITER: 81247 RES: 2.7478e-06
Jacobi stage: ITER: 81573 RES: 2.6128e-06
Jacobi stage: ITER: 81899 RES: 2.486e-06
Jacobi stage: ITER: 82227 RES: 2.3638e-06
Jacobi stage: ITER: 82551 RES: 2.2491e-06
Jacobi stage: ITER: 82879 RES: 2.1385e-06
Jacobi stage: ITER: 83205 RES: 2.0335e-06
Jacobi stage: ITER: 83531 RES: 1.9348e-06
Jacobi stage: ITER: 83857 RES: 1.8397e-06
Jacobi stage: ITER: 84183 RES: 1.7504e-06
Jacobi stage: ITER: 84509 RES: 1.6644e-06
Jacobi stage: ITER: 84835 RES: 1.5836e-06
Jacobi stage: ITER: 85163 RES: 1.5058e-06
Jacobi stage: ITER: 85487 RES: 1.4327e-06
Jacobi stage: ITER: 85815 RES: 1.3623e-06
Jacobi stage: ITER: 86141 RES: 1.2953e-06
Jacobi stage: ITER: 86467 RES: 1.2325e-06
Jacobi stage: ITER: 86795 RES: 1.1719e-06
Jacobi stage: ITER: 87119 RES: 1.115e-06
Jacobi stage: ITER: 87447 RES: 1.0602e-06
Jacobi stage: ITER: 87773 RES: 1.0081e-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/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
13iterativeLinearSolverExample()
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 */
24 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
25 using Vector = TNL::Containers::Vector< double, Device >;
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 {
33 const int rowIdx = row.getRowIndex();
34 if( rowIdx == 0 ) {
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 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
40 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
41 }
42 else {
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 * Setup solver of the linear system.
66 */
67 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
68 LinearSolver solver;
69 solver.setMatrix( matrix_ptr );
70 solver.setOmega( 0.0005 );
71
72 /***
73 * Setup monitor of the iterative solver.
74 */
75 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
76 IterativeSolverMonitorType monitor;
77 TNL::Solvers::SolverMonitorThread mmonitorThread( monitor );
78 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
79 monitor.setVerbose( 1 );
80 monitor.setStage( "Jacobi stage:" );
81 TNL::Timer timer;
82 monitor.setTimer( timer );
83 timer.start();
84 solver.setSolverMonitor( monitor );
85 solver.setConvergenceResidue( 1.0e-6 );
86 solver.solve( b, x );
87 monitor.stopMainLoop();
88 std::cout << "Vector x = " << x << std::endl;
89}
90
91int
92main( int argc, char* argv[] )
93{
94 std::cout << "Solving linear system on host: " << std::endl;
95 iterativeLinearSolverExample< TNL::Devices::Sequential >();
96
97#ifdef __CUDACC__
98 std::cout << "Solving linear system on CUDA device: " << std::endl;
99 iterativeLinearSolverExample< TNL::Devices::Cuda >();
100#endif
101}
TNL::Timer
Class for real time, CPU time and CPU cycles measuring.
Definition Timer.h:25
TNL::Timer::start
void start()
Starts timer.
Definition Timer.hpp:42

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:1.3022e-
ELA: 0.10011 Jacobi stage: ITER: 26835 RES: 0.011707
ELA: 0.20022 Jacobi stage: ITER: 54557 RES: 0.00016569
ELA: 0.30033 Jacobi stage: ITER: 83925 RES: 1.8206e-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.5538e-
ELA: 0.10013 Jacobi stage: ITER: 383 RES: 0.88222
ELA: 0.20025 Jacobi stage: ITER: 721 RES: 0.79705
ELA: 0.30038 Jacobi stage: ITER: 1043 RES: 0.73017
ELA: 0.4005 Jacobi stage: ITER: 1367 RES: 0.67269
ELA: 0.50062 Jacobi stage: ITER: 1691 RES: 0.62332
ELA: 0.60074 Jacobi stage: ITER: 2015 RES: 0.58037
ELA: 0.70086 Jacobi stage: ITER: 2337 RES: 0.54257
ELA: 0.80096 Jacobi stage: ITER: 2661 RES: 0.50891
ELA: 0.90109 Jacobi stage: ITER: 2983 RES: 0.47901
ELA: 1.0012 Jacobi stage: ITER: 3307 RES: 0.45152
ELA: 1.1014 Jacobi stage: ITER: 3631 RES: 0.42639
ELA: 1.2015 Jacobi stage: ITER: 3953 RES: 0.40325
ELA: 1.3016 Jacobi stage: ITER: 4275 RES: 0.38181
ELA: 1.4017 Jacobi stage: ITER: 4599 RES: 0.3621
ELA: 1.5018 Jacobi stage: ITER: 4923 RES: 0.34345
ELA: 1.6019 Jacobi stage: ITER: 5245 RES: 0.32596
ELA: 1.7021 Jacobi stage: ITER: 5569 RES: 0.30951
ELA: 1.8022 Jacobi stage: ITER: 5891 RES: 0.2942
ELA: 1.9023 Jacobi stage: ITER: 6215 RES: 0.27956
ELA: 2.0024 Jacobi stage: ITER: 6539 RES: 0.26571
ELA: 2.1025 Jacobi stage: ITER: 6863 RES: 0.2526
ELA: 2.2026 Jacobi stage: ITER: 7187 RES: 0.24018
ELA: 2.3027 Jacobi stage: ITER: 7511 RES: 0.2284
ELA: 2.4029 Jacobi stage: ITER: 7797 RES: 0.21844
ELA: 2.503 Jacobi stage: ITER: 8119 RES: 0.20789
ELA: 2.6031 Jacobi stage: ITER: 8447 RES: 0.19762
ELA: 2.7032 Jacobi stage: ITER: 8775 RES: 0.18787
ELA: 2.8033 Jacobi stage: ITER: 9101 RES: 0.17861
ELA: 2.9035 Jacobi stage: ITER: 9407 RES: 0.17043
ELA: 3.0036 Jacobi stage: ITER: 9721 RES: 0.16234
ELA: 3.1037 Jacobi stage: ITER: 10051 RES: 0.15435
ELA: 3.2038 Jacobi stage: ITER: 10379 RES: 0.14676
ELA: 3.3039 Jacobi stage: ITER: 10707 RES: 0.13954
ELA: 3.404 Jacobi stage: ITER: 11033 RES: 0.13268
ELA: 3.5041 Jacobi stage: ITER: 11351 RES: 0.12639
ELA: 3.6043 Jacobi stage: ITER: 11651 RES: 0.12069
ELA: 3.7044 Jacobi stage: ITER: 11971 RES: 0.1149
ELA: 3.8045 Jacobi stage: ITER: 12275 RES: 0.10965
ELA: 3.9046 Jacobi stage: ITER: 12599 RES: 0.10433
ELA: 4.0047 Jacobi stage: ITER: 12919 RES: 0.099324
ELA: 4.1048 Jacobi stage: ITER: 13249 RES: 0.094385
ELA: 4.2049 Jacobi stage: ITER: 13567 RES: 0.089913
ELA: 4.305 Jacobi stage: ITER: 13887 RES: 0.0856
ELA: 4.4052 Jacobi stage: ITER: 14215 RES: 0.081393
ELA: 4.5053 Jacobi stage: ITER: 14547 RES: 0.077346
ELA: 4.6054 Jacobi stage: ITER: 14843 RES: 0.073908
ELA: 4.7055 Jacobi stage: ITER: 15163 RES: 0.070363
ELA: 4.8056 Jacobi stage: ITER: 15469 RES: 0.067112
ELA: 4.9057 Jacobi stage: ITER: 15795 RES: 0.063853
ELA: 5.0058 Jacobi stage: ITER: 16121 RES: 0.060716
ELA: 5.1059 Jacobi stage: ITER: 16451 RES: 0.057733
ELA: 5.206 Jacobi stage: ITER: 16763 RES: 0.055031
ELA: 5.3061 Jacobi stage: ITER: 17085 RES: 0.05236
ELA: 5.4062 Jacobi stage: ITER: 17413 RES: 0.049787
ELA: 5.5063 Jacobi stage: ITER: 17739 RES: 0.04737
ELA: 5.6064 Jacobi stage: ITER: 18071 RES: 0.045015
ELA: 5.7066 Jacobi stage: ITER: 18369 RES: 0.042987
ELA: 5.8067 Jacobi stage: ITER: 18667 RES: 0.041077
ELA: 5.9068 Jacobi stage: ITER: 18993 RES: 0.039058
ELA: 6.007 Jacobi stage: ITER: 19319 RES: 0.037162
ELA: 6.1071 Jacobi stage: ITER: 19647 RES: 0.035336
ELA: 6.2072 Jacobi stage: ITER: 19965 RES: 0.033641
ELA: 6.3073 Jacobi stage: ITER: 20287 RES: 0.032028
ELA: 6.4074 Jacobi stage: ITER: 20613 RES: 0.030454
ELA: 6.5075 Jacobi stage: ITER: 20941 RES: 0.028958
ELA: 6.6076 Jacobi stage: ITER: 21273 RES: 0.027518
ELA: 6.7077 Jacobi stage: ITER: 21537 RES: 0.026425
ELA: 6.8078 Jacobi stage: ITER: 21835 RES: 0.02525
ELA: 6.908 Jacobi stage: ITER: 22159 RES: 0.024024
ELA: 7.0081 Jacobi stage: ITER: 22483 RES: 0.022858
ELA: 7.1082 Jacobi stage: ITER: 22811 RES: 0.021735
ELA: 7.2083 Jacobi stage: ITER: 23143 RES: 0.020654
ELA: 7.3084 Jacobi stage: ITER: 23459 RES: 0.019676
ELA: 7.4085 Jacobi stage: ITER: 23781 RES: 0.01872
ELA: 7.5086 Jacobi stage: ITER: 24107 RES: 0.017812
ELA: 7.6088 Jacobi stage: ITER: 24437 RES: 0.016926
ELA: 7.7089 Jacobi stage: ITER: 24769 RES: 0.016084
ELA: 7.809 Jacobi stage: ITER: 25063 RES: 0.015379
ELA: 7.9091 Jacobi stage: ITER: 25371 RES: 0.014668
ELA: 8.0092 Jacobi stage: ITER: 25701 RES: 0.013939
ELA: 8.1093 Jacobi stage: ITER: 26031 RES: 0.013254
ELA: 8.2094 Jacobi stage: ITER: 26361 RES: 0.012595
ELA: 8.3095 Jacobi stage: ITER: 26679 RES: 0.011999
ELA: 8.4096 Jacobi stage: ITER: 26999 RES: 0.011423
ELA: 8.5097 Jacobi stage: ITER: 27325 RES: 0.010862
ELA: 8.6098 Jacobi stage: ITER: 27651 RES: 0.010334
ELA: 8.71 Jacobi stage: ITER: 27979 RES: 0.0098267
ELA: 8.8101 Jacobi stage: ITER: 28313 RES: 0.0093324
ELA: 8.9102 Jacobi stage: ITER: 28595 RES: 0.0089395
ELA: 9.0103 Jacobi stage: ITER: 28887 RES: 0.0085474
ELA: 9.1104 Jacobi stage: ITER: 29207 RES: 0.0081375
ELA: 9.2105 Jacobi stage: ITER: 29535 RES: 0.0077376
ELA: 9.3106 Jacobi stage: ITER: 29863 RES: 0.0073575
ELA: 9.4108 Jacobi stage: ITER: 30195 RES: 0.0069917
ELA: 9.5109 Jacobi stage: ITER: 30507 RES: 0.0066645
ELA: 9.611 Jacobi stage: ITER: 30825 RES: 0.0063449
ELA: 9.7128 Jacobi stage: ITER: 31159 RES: 0.0060294
ELA: 9.8129 Jacobi stage: ITER: 31481 RES: 0.0057367
ELA: 9.913 Jacobi stage: ITER: 31785 RES: 0.005475
ELA: 10.015 Jacobi stage: ITER: 32095 RES: 0.005222
ELA: 10.115 Jacobi stage: ITER: 32423 RES: 0.0049654
ELA: 10.215 Jacobi stage: ITER: 32751 RES: 0.0047214
ELA: 10.315 Jacobi stage: ITER: 33077 RES: 0.0044895
ELA: 10.415 Jacobi stage: ITER: 33403 RES: 0.0042715
ELA: 10.515 Jacobi stage: ITER: 33729 RES: 0.0040616
ELA: 10.616 Jacobi stage: ITER: 34055 RES: 0.0038644
ELA: 10.716 Jacobi stage: ITER: 34381 RES: 0.0036746
ELA: 10.816 Jacobi stage: ITER: 34707 RES: 0.0034962
ELA: 10.916 Jacobi stage: ITER: 35035 RES: 0.0033244
ELA: 11.016 Jacobi stage: ITER: 35361 RES: 0.0031611
ELA: 11.116 Jacobi stage: ITER: 35687 RES: 0.0030076
ELA: 11.216 Jacobi stage: ITER: 36013 RES: 0.0028598
ELA: 11.322 Jacobi stage: ITER: 36351 RES: 0.002716
ELA: 11.422 Jacobi stage: ITER: 36643 RES: 0.0025968
ELA: 11.522 Jacobi stage: ITER: 36969 RES: 0.0024693
ELA: 11.622 Jacobi stage: ITER: 37295 RES: 0.0023494
ELA: 11.722 Jacobi stage: ITER: 37621 RES: 0.0022339
ELA: 11.822 Jacobi stage: ITER: 37947 RES: 0.0021255
ELA: 11.922 Jacobi stage: ITER: 38273 RES: 0.0020211
ELA: 12.022 Jacobi stage: ITER: 38599 RES: 0.0019229
ELA: 12.122 Jacobi stage: ITER: 38927 RES: 0.0018285
ELA: 12.222 Jacobi stage: ITER: 39251 RES: 0.0017397
ELA: 12.323 Jacobi stage: ITER: 39579 RES: 0.0016542
ELA: 12.423 Jacobi stage: ITER: 39905 RES: 0.0015729
ELA: 12.523 Jacobi stage: ITER: 40231 RES: 0.0014966
ELA: 12.623 Jacobi stage: ITER: 40557 RES: 0.001423
ELA: 12.723 Jacobi stage: ITER: 40883 RES: 0.001354
ELA: 12.823 Jacobi stage: ITER: 41209 RES: 0.0012874
ELA: 12.923 Jacobi stage: ITER: 41535 RES: 0.0012249
ELA: 13.023 Jacobi stage: ITER: 41863 RES: 0.0011647
ELA: 13.123 Jacobi stage: ITER: 42187 RES: 0.0011082
ELA: 13.223 Jacobi stage: ITER: 42515 RES: 0.0010537
ELA: 13.323 Jacobi stage: ITER: 42841 RES: 0.001002
ELA: 13.424 Jacobi stage: ITER: 43167 RES: 0.00095332
ELA: 13.524 Jacobi stage: ITER: 43493 RES: 0.00090648
ELA: 13.624 Jacobi stage: ITER: 43819 RES: 0.00086248
ELA: 13.724 Jacobi stage: ITER: 44147 RES: 0.0008201
ELA: 13.824 Jacobi stage: ITER: 44471 RES: 0.00078028
ELA: 13.924 Jacobi stage: ITER: 44799 RES: 0.00074195
ELA: 14.024 Jacobi stage: ITER: 45123 RES: 0.00070549
ELA: 14.124 Jacobi stage: ITER: 45451 RES: 0.00067124
ELA: 14.224 Jacobi stage: ITER: 45777 RES: 0.00063826
ELA: 14.325 Jacobi stage: ITER: 46103 RES: 0.00060727
ELA: 14.425 Jacobi stage: ITER: 46429 RES: 0.00057744
ELA: 14.525 Jacobi stage: ITER: 46755 RES: 0.0005494
ELA: 14.625 Jacobi stage: ITER: 47083 RES: 0.00052241
ELA: 14.725 Jacobi stage: ITER: 47407 RES: 0.00049705
ELA: 14.825 Jacobi stage: ITER: 47735 RES: 0.00047263
ELA: 14.925 Jacobi stage: ITER: 48061 RES: 0.0004494
ELA: 15.025 Jacobi stage: ITER: 48387 RES: 0.00042759
ELA: 15.125 Jacobi stage: ITER: 48713 RES: 0.00040658
ELA: 15.225 Jacobi stage: ITER: 49039 RES: 0.00038684
ELA: 15.326 Jacobi stage: ITER: 49365 RES: 0.00036783
ELA: 15.426 Jacobi stage: ITER: 49691 RES: 0.00034997
ELA: 15.526 Jacobi stage: ITER: 50019 RES: 0.00033278
ELA: 15.626 Jacobi stage: ITER: 50343 RES: 0.00031662
ELA: 15.726 Jacobi stage: ITER: 50671 RES: 0.00030107
ELA: 15.826 Jacobi stage: ITER: 50997 RES: 0.00028627
ELA: 15.926 Jacobi stage: ITER: 51323 RES: 0.00027238
ELA: 16.026 Jacobi stage: ITER: 51649 RES: 0.00025899
ELA: 16.126 Jacobi stage: ITER: 51975 RES: 0.00024642
ELA: 16.227 Jacobi stage: ITER: 52301 RES: 0.00023431
ELA: 16.327 Jacobi stage: ITER: 52627 RES: 0.00022294
ELA: 16.427 Jacobi stage: ITER: 52955 RES: 0.00021198
ELA: 16.527 Jacobi stage: ITER: 53279 RES: 0.00020169
ELA: 16.627 Jacobi stage: ITER: 53607 RES: 0.00019178
ELA: 16.727 Jacobi stage: ITER: 53933 RES: 0.00018236
ELA: 16.827 Jacobi stage: ITER: 54259 RES: 0.0001735
ELA: 16.927 Jacobi stage: ITER: 54585 RES: 0.00016498
ELA: 17.027 Jacobi stage: ITER: 54911 RES: 0.00015697
ELA: 17.127 Jacobi stage: ITER: 55237 RES: 0.00014926
ELA: 17.227 Jacobi stage: ITER: 55563 RES: 0.00014201
ELA: 17.328 Jacobi stage: ITER: 55891 RES: 0.00013503
ELA: 17.428 Jacobi stage: ITER: 56215 RES: 0.00012848
ELA: 17.528 Jacobi stage: ITER: 56543 RES: 0.00012217
ELA: 17.628 Jacobi stage: ITER: 56867 RES: 0.00011623
ELA: 17.728 Jacobi stage: ITER: 57195 RES: 0.00011052
ELA: 17.828 Jacobi stage: ITER: 57521 RES: 0.00010509
ELA: 17.928 Jacobi stage: ITER: 57847 RES: 9.9991e-05
ELA: 18.028 Jacobi stage: ITER: 58173 RES: 9.5078e-05
ELA: 18.128 Jacobi stage: ITER: 58499 RES: 9.0462e-05
ELA: 18.228 Jacobi stage: ITER: 58825 RES: 8.6018e-05
ELA: 18.328 Jacobi stage: ITER: 59151 RES: 8.1842e-05
ELA: 18.428 Jacobi stage: ITER: 59479 RES: 7.7821e-05
ELA: 18.528 Jacobi stage: ITER: 59803 RES: 7.4042e-05
ELA: 18.629 Jacobi stage: ITER: 60131 RES: 7.0405e-05
ELA: 18.729 Jacobi stage: ITER: 60455 RES: 6.6986e-05
ELA: 18.829 Jacobi stage: ITER: 60783 RES: 6.3695e-05
ELA: 18.929 Jacobi stage: ITER: 61109 RES: 6.0566e-05
ELA: 19.029 Jacobi stage: ITER: 61435 RES: 5.7625e-05
ELA: 19.129 Jacobi stage: ITER: 61761 RES: 5.4794e-05
ELA: 19.229 Jacobi stage: ITER: 62087 RES: 5.2134e-05
ELA: 19.329 Jacobi stage: ITER: 62415 RES: 4.9572e-05
ELA: 19.429 Jacobi stage: ITER: 62739 RES: 4.7166e-05
ELA: 19.53 Jacobi stage: ITER: 63067 RES: 4.4848e-05
ELA: 19.63 Jacobi stage: ITER: 63393 RES: 4.2645e-05
ELA: 19.73 Jacobi stage: ITER: 63719 RES: 4.0574e-05
ELA: 19.83 Jacobi stage: ITER: 64045 RES: 3.8581e-05
ELA: 19.93 Jacobi stage: ITER: 64371 RES: 3.6708e-05
ELA: 20.03 Jacobi stage: ITER: 64699 RES: 3.4904e-05
ELA: 20.13 Jacobi stage: ITER: 65025 RES: 3.3189e-05
ELA: 20.231 Jacobi stage: ITER: 65351 RES: 3.1578e-05
ELA: 20.331 Jacobi stage: ITER: 65677 RES: 3.0026e-05
ELA: 20.431 Jacobi stage: ITER: 66003 RES: 2.8569e-05
ELA: 20.531 Jacobi stage: ITER: 66331 RES: 2.7165e-05
ELA: 20.631 Jacobi stage: ITER: 66657 RES: 2.583e-05
ELA: 20.731 Jacobi stage: ITER: 66983 RES: 2.4576e-05
ELA: 20.831 Jacobi stage: ITER: 67309 RES: 2.3369e-05
ELA: 20.932 Jacobi stage: ITER: 67635 RES: 2.2234e-05
ELA: 21.032 Jacobi stage: ITER: 67963 RES: 2.1142e-05
ELA: 21.132 Jacobi stage: ITER: 68289 RES: 2.0103e-05
ELA: 21.232 Jacobi stage: ITER: 68615 RES: 1.9127e-05
ELA: 21.332 Jacobi stage: ITER: 68941 RES: 1.8187e-05
ELA: 21.432 Jacobi stage: ITER: 69267 RES: 1.7304e-05
ELA: 21.532 Jacobi stage: ITER: 69595 RES: 1.6454e-05
ELA: 21.632 Jacobi stage: ITER: 69919 RES: 1.5655e-05
ELA: 21.733 Jacobi stage: ITER: 70247 RES: 1.4886e-05
ELA: 21.833 Jacobi stage: ITER: 70573 RES: 1.4155e-05
ELA: 21.933 Jacobi stage: ITER: 70899 RES: 1.3467e-05
ELA: 22.033 Jacobi stage: ITER: 71225 RES: 1.2806e-05
ELA: 22.133 Jacobi stage: ITER: 71551 RES: 1.2184e-05
ELA: 22.233 Jacobi stage: ITER: 71879 RES: 1.1585e-05
ELA: 22.333 Jacobi stage: ITER: 72203 RES: 1.1023e-05
ELA: 22.433 Jacobi stage: ITER: 72531 RES: 1.0481e-05
ELA: 22.534 Jacobi stage: ITER: 72857 RES: 9.9664e-06
ELA: 22.634 Jacobi stage: ITER: 73183 RES: 9.4825e-06
ELA: 22.734 Jacobi stage: ITER: 73511 RES: 9.0166e-06
ELA: 22.834 Jacobi stage: ITER: 73835 RES: 8.5789e-06
ELA: 22.934 Jacobi stage: ITER: 74163 RES: 8.1573e-06
ELA: 23.034 Jacobi stage: ITER: 74489 RES: 7.7566e-06
ELA: 23.135 Jacobi stage: ITER: 74815 RES: 7.38e-06
ELA: 23.235 Jacobi stage: ITER: 75143 RES: 7.0174e-06
ELA: 23.335 Jacobi stage: ITER: 75469 RES: 6.6726e-06
ELA: 23.435 Jacobi stage: ITER: 75795 RES: 6.3486e-06
ELA: 23.535 Jacobi stage: ITER: 76121 RES: 6.0367e-06
ELA: 23.635 Jacobi stage: ITER: 76447 RES: 5.7436e-06
ELA: 23.735 Jacobi stage: ITER: 76775 RES: 5.4614e-06
ELA: 23.835 Jacobi stage: ITER: 77099 RES: 5.1963e-06
ELA: 23.935 Jacobi stage: ITER: 77427 RES: 4.941e-06
ELA: 24.036 Jacobi stage: ITER: 77753 RES: 4.6982e-06
ELA: 24.136 Jacobi stage: ITER: 78079 RES: 4.4701e-06
ELA: 24.236 Jacobi stage: ITER: 78405 RES: 4.2505e-06
ELA: 24.336 Jacobi stage: ITER: 78731 RES: 4.0441e-06
ELA: 24.436 Jacobi stage: ITER: 79057 RES: 3.8454e-06
ELA: 24.536 Jacobi stage: ITER: 79383 RES: 3.6587e-06
ELA: 24.636 Jacobi stage: ITER: 79711 RES: 3.479e-06
ELA: 24.736 Jacobi stage: ITER: 80037 RES: 3.308e-06
ELA: 24.836 Jacobi stage: ITER: 80363 RES: 3.1474e-06
ELA: 24.937 Jacobi stage: ITER: 80689 RES: 2.9928e-06
ELA: 25.037 Jacobi stage: ITER: 81015 RES: 2.8475e-06
ELA: 25.137 Jacobi stage: ITER: 81341 RES: 2.7076e-06
ELA: 25.237 Jacobi stage: ITER: 81667 RES: 2.5761e-06
ELA: 25.337 Jacobi stage: ITER: 81995 RES: 2.4496e-06
ELA: 25.437 Jacobi stage: ITER: 82321 RES: 2.3292e-06
ELA: 25.537 Jacobi stage: ITER: 82647 RES: 2.2161e-06
ELA: 25.638 Jacobi stage: ITER: 82973 RES: 2.1072e-06
ELA: 25.738 Jacobi stage: ITER: 83299 RES: 2.0049e-06
ELA: 25.838 Jacobi stage: ITER: 83627 RES: 1.9064e-06
ELA: 25.938 Jacobi stage: ITER: 84051 RES: 1.7862e-06
ELA: 26.038 Jacobi stage: ITER: 84545 RES: 1.6552e-06
ELA: 26.138 Jacobi stage: ITER: 85035 RES: 1.5357e-06
ELA: 26.238 Jacobi stage: ITER: 85527 RES: 1.4239e-06
ELA: 26.339 Jacobi stage: ITER: 86035 RES: 1.317e-06
ELA: 26.439 Jacobi stage: ITER: 86527 RES: 1.2211e-06
ELA: 26.539 Jacobi stage: ITER: 87017 RES: 1.1323e-06
ELA: 26.639 Jacobi stage: ITER: 87507 RES: 1.0505e-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 << std::endl;
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 << std::endl;
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 << std::endl;
75}
76
77int
78main( int argc, char* argv[] )
79{
80 std::cout << "Solving linear system on host: " << std::endl;
81 iterativeLinearSolverExample< TNL::Devices::Sequential >();
82
83#ifdef __CUDACC__
84 std::cout << "Solving linear system on CUDA device: " << std::endl;
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 << std::endl;
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 << std::endl;
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 << std::endl;
72}
73
74int
75main( int argc, char* argv[] )
76{
77 std::cout << "Solving linear system on host: " << std::endl;
78 iterativeLinearSolverExample< TNL::Devices::Sequential >();
79
80#ifdef __CUDACC__
81 std::cout << "Solving linear system on CUDA device: " << std::endl;
82 iterativeLinearSolverExample< TNL::Devices::Cuda >();
83#endif
84}
TNL::Solvers::getLinearSolver
std::shared_ptr< Linear::LinearSolver< MatrixType > > getLinearSolver(std::string name)
Function returning shared pointer with linear solver given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:84
TNL::Solvers::getPreconditioner
std::shared_ptr< Linear::Preconditioners::Preconditioner< MatrixType > > getPreconditioner(std::string name)
Function returning shared pointer with linear preconditioner given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:138

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 ]
Generated on Wed Dec 11 2024 22:03:19 for Template Numerical Library by doxygen 1.12.0