Template Numerical Library version\ main:40802e2e
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: 29847 RES: 0.0073756
Jacobi stage: ITER: 60183 RES: 6.9844e-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: 5353 RES: 0.32037
Jacobi stage: ITER: 9891 RES: 0.1582
Jacobi stage: ITER: 10275 RES: 0.14912
Jacobi stage: ITER: 10603 RES: 0.14179
Jacobi stage: ITER: 10927 RES: 0.1349
Jacobi stage: ITER: 11255 RES: 0.12826
Jacobi stage: ITER: 11579 RES: 0.12203
Jacobi stage: ITER: 11907 RES: 0.11603
Jacobi stage: ITER: 12231 RES: 0.1104
Jacobi stage: ITER: 12559 RES: 0.10497
Jacobi stage: ITER: 12885 RES: 0.099814
Jacobi stage: ITER: 13199 RES: 0.095142
Jacobi stage: ITER: 13511 RES: 0.09069
Jacobi stage: ITER: 13837 RES: 0.086233
Jacobi stage: ITER: 14163 RES: 0.082046
Jacobi stage: ITER: 14491 RES: 0.078015
Jacobi stage: ITER: 14817 RES: 0.074181
Jacobi stage: ITER: 15131 RES: 0.07071
Jacobi stage: ITER: 15461 RES: 0.067194
Jacobi stage: ITER: 15791 RES: 0.063893
Jacobi stage: ITER: 16067 RES: 0.061241
Jacobi stage: ITER: 16375 RES: 0.058411
Jacobi stage: ITER: 16703 RES: 0.055541
Jacobi stage: ITER: 17035 RES: 0.052779
Jacobi stage: ITER: 17365 RES: 0.050155
Jacobi stage: ITER: 17689 RES: 0.04772
Jacobi stage: ITER: 18003 RES: 0.045487
Jacobi stage: ITER: 18327 RES: 0.043279
Jacobi stage: ITER: 18651 RES: 0.041178
Jacobi stage: ITER: 18979 RES: 0.039155
Jacobi stage: ITER: 19311 RES: 0.037208
Jacobi stage: ITER: 19631 RES: 0.035423
Jacobi stage: ITER: 19939 RES: 0.033786
Jacobi stage: ITER: 20267 RES: 0.032126
Jacobi stage: ITER: 20595 RES: 0.030548
Jacobi stage: ITER: 20907 RES: 0.029118
Jacobi stage: ITER: 21219 RES: 0.027756
Jacobi stage: ITER: 21547 RES: 0.026392
Jacobi stage: ITER: 21875 RES: 0.025095
Jacobi stage: ITER: 22207 RES: 0.023848
Jacobi stage: ITER: 22539 RES: 0.022662
Jacobi stage: ITER: 22843 RES: 0.021628
Jacobi stage: ITER: 23161 RES: 0.020591
Jacobi stage: ITER: 23491 RES: 0.019579
Jacobi stage: ITER: 23817 RES: 0.018617
Jacobi stage: ITER: 24127 RES: 0.017757
Jacobi stage: ITER: 24439 RES: 0.016926
Jacobi stage: ITER: 24765 RES: 0.016094
Jacobi stage: ITER: 25097 RES: 0.015294
Jacobi stage: ITER: 25429 RES: 0.014534
Jacobi stage: ITER: 25739 RES: 0.013862
Jacobi stage: ITER: 26051 RES: 0.013214
Jacobi stage: ITER: 26383 RES: 0.012557
Jacobi stage: ITER: 26713 RES: 0.011932
Jacobi stage: ITER: 27047 RES: 0.011339
Jacobi stage: ITER: 27345 RES: 0.010828
Jacobi stage: ITER: 27665 RES: 0.010309
Jacobi stage: ITER: 27991 RES: 0.0098086
Jacobi stage: ITER: 28323 RES: 0.0093209
Jacobi stage: ITER: 28659 RES: 0.0088521
Jacobi stage: ITER: 28957 RES: 0.0084534
Jacobi stage: ITER: 29273 RES: 0.0080529
Jacobi stage: ITER: 29603 RES: 0.0076572
Jacobi stage: ITER: 29935 RES: 0.0072766
Jacobi stage: ITER: 30247 RES: 0.0069361
Jacobi stage: ITER: 30561 RES: 0.0066074
Jacobi stage: ITER: 30891 RES: 0.0062828
Jacobi stage: ITER: 31223 RES: 0.0059704
Jacobi stage: ITER: 31553 RES: 0.0056736
Jacobi stage: ITER: 31887 RES: 0.0053915
Jacobi stage: ITER: 32209 RES: 0.0051298
Jacobi stage: ITER: 32499 RES: 0.0049078
Jacobi stage: ITER: 32821 RES: 0.0046695
Jacobi stage: ITER: 33149 RES: 0.0044401
Jacobi stage: ITER: 33479 RES: 0.0042219
Jacobi stage: ITER: 33809 RES: 0.004012
Jacobi stage: ITER: 34127 RES: 0.0038219
Jacobi stage: ITER: 34439 RES: 0.0036431
Jacobi stage: ITER: 34751 RES: 0.0034726
Jacobi stage: ITER: 35079 RES: 0.003302
Jacobi stage: ITER: 35411 RES: 0.0031378
Jacobi stage: ITER: 35743 RES: 0.0029818
Jacobi stage: ITER: 35991 RES: 0.0028704
Jacobi stage: ITER: 36299 RES: 0.0027377
Jacobi stage: ITER: 36625 RES: 0.0026032
Jacobi stage: ITER: 36957 RES: 0.0024738
Jacobi stage: ITER: 37279 RES: 0.0023552
Jacobi stage: ITER: 37597 RES: 0.0022422
Jacobi stage: ITER: 37931 RES: 0.0021307
Jacobi stage: ITER: 38223 RES: 0.0020373
Jacobi stage: ITER: 38553 RES: 0.001936
Jacobi stage: ITER: 38883 RES: 0.0018409
Jacobi stage: ITER: 39211 RES: 0.0017504
Jacobi stage: ITER: 39539 RES: 0.0016644
Jacobi stage: ITER: 39869 RES: 0.0015817
Jacobi stage: ITER: 40199 RES: 0.0015039
Jacobi stage: ITER: 40527 RES: 0.00143
Jacobi stage: ITER: 40855 RES: 0.0013598
Jacobi stage: ITER: 41185 RES: 0.0012922
Jacobi stage: ITER: 41515 RES: 0.0012287
Jacobi stage: ITER: 41843 RES: 0.0011683
Jacobi stage: ITER: 42171 RES: 0.0011109
Jacobi stage: ITER: 42501 RES: 0.0010557
Jacobi stage: ITER: 42831 RES: 0.0010038
Jacobi stage: ITER: 43159 RES: 0.0009545
Jacobi stage: ITER: 43487 RES: 0.0009076
Jacobi stage: ITER: 43815 RES: 0.00086301
Jacobi stage: ITER: 44145 RES: 0.0008201
Jacobi stage: ITER: 44475 RES: 0.00077981
Jacobi stage: ITER: 44803 RES: 0.00074149
Jacobi stage: ITER: 45131 RES: 0.00070506
Jacobi stage: ITER: 45461 RES: 0.00067001
Jacobi stage: ITER: 45791 RES: 0.00063709
Jacobi stage: ITER: 46119 RES: 0.00060578
Jacobi stage: ITER: 46447 RES: 0.00057602
Jacobi stage: ITER: 46775 RES: 0.00054772
Jacobi stage: ITER: 47105 RES: 0.00052049
Jacobi stage: ITER: 47435 RES: 0.00049491
Jacobi stage: ITER: 47763 RES: 0.0004706
Jacobi stage: ITER: 48091 RES: 0.00044748
Jacobi stage: ITER: 48421 RES: 0.00042523
Jacobi stage: ITER: 48751 RES: 0.00040434
Jacobi stage: ITER: 49079 RES: 0.00038447
Jacobi stage: ITER: 49407 RES: 0.00036558
Jacobi stage: ITER: 49737 RES: 0.0003474
Jacobi stage: ITER: 50067 RES: 0.00033033
Jacobi stage: ITER: 50395 RES: 0.0003141
Jacobi stage: ITER: 50723 RES: 0.00029867
Jacobi stage: ITER: 51053 RES: 0.00028382
Jacobi stage: ITER: 51381 RES: 0.00026988
Jacobi stage: ITER: 51711 RES: 0.00025662
Jacobi stage: ITER: 52039 RES: 0.00024401
Jacobi stage: ITER: 52367 RES: 0.00023202
Jacobi stage: ITER: 52697 RES: 0.00022048
Jacobi stage: ITER: 53027 RES: 0.00020965
Jacobi stage: ITER: 53355 RES: 0.00019935
Jacobi stage: ITER: 53683 RES: 0.00018956
Jacobi stage: ITER: 54013 RES: 0.00018013
Jacobi stage: ITER: 54343 RES: 0.00017128
Jacobi stage: ITER: 54671 RES: 0.00016287
Jacobi stage: ITER: 54999 RES: 0.00015486
Jacobi stage: ITER: 55327 RES: 0.00014725
Jacobi stage: ITER: 55657 RES: 0.00013993
Jacobi stage: ITER: 55987 RES: 0.00013306
Jacobi stage: ITER: 56315 RES: 0.00012652
Jacobi stage: ITER: 56643 RES: 0.0001203
Jacobi stage: ITER: 56973 RES: 0.00011432
Jacobi stage: ITER: 57301 RES: 0.00010871
Jacobi stage: ITER: 57631 RES: 0.00010336
Jacobi stage: ITER: 57959 RES: 9.8286e-05
Jacobi stage: ITER: 58287 RES: 9.3457e-05
Jacobi stage: ITER: 58617 RES: 8.881e-05
Jacobi stage: ITER: 58945 RES: 8.4447e-05
Jacobi stage: ITER: 59275 RES: 8.0298e-05
Jacobi stage: ITER: 59603 RES: 7.6352e-05
Jacobi stage: ITER: 59931 RES: 7.2601e-05
Jacobi stage: ITER: 60261 RES: 6.8991e-05
Jacobi stage: ITER: 60591 RES: 6.5602e-05
Jacobi stage: ITER: 60919 RES: 6.2378e-05
Jacobi stage: ITER: 61247 RES: 5.9314e-05
Jacobi stage: ITER: 61577 RES: 5.6365e-05
Jacobi stage: ITER: 61907 RES: 5.3595e-05
Jacobi stage: ITER: 62235 RES: 5.0962e-05
Jacobi stage: ITER: 62563 RES: 4.8458e-05
Jacobi stage: ITER: 62893 RES: 4.6049e-05
Jacobi stage: ITER: 63233 RES: 4.3706e-05
Jacobi stage: ITER: 63563 RES: 4.1558e-05
Jacobi stage: ITER: 63891 RES: 3.9516e-05
Jacobi stage: ITER: 64221 RES: 3.7552e-05
Jacobi stage: ITER: 64549 RES: 3.5707e-05
Jacobi stage: ITER: 64879 RES: 3.3952e-05
Jacobi stage: ITER: 65207 RES: 3.2284e-05
Jacobi stage: ITER: 65535 RES: 3.0698e-05
Jacobi stage: ITER: 65865 RES: 2.9172e-05
Jacobi stage: ITER: 66193 RES: 2.7738e-05
Jacobi stage: ITER: 66523 RES: 2.6376e-05
Jacobi stage: ITER: 66851 RES: 2.508e-05
Jacobi stage: ITER: 67179 RES: 2.3847e-05
Jacobi stage: ITER: 67509 RES: 2.2662e-05
Jacobi stage: ITER: 67837 RES: 2.1548e-05
Jacobi stage: ITER: 68167 RES: 2.049e-05
Jacobi stage: ITER: 68495 RES: 1.9483e-05
Jacobi stage: ITER: 68823 RES: 1.8526e-05
Jacobi stage: ITER: 69153 RES: 1.7605e-05
Jacobi stage: ITER: 69481 RES: 1.674e-05
Jacobi stage: ITER: 69811 RES: 1.5917e-05
Jacobi stage: ITER: 70139 RES: 1.5135e-05
Jacobi stage: ITER: 70469 RES: 1.4383e-05
Jacobi stage: ITER: 70797 RES: 1.3676e-05
Jacobi stage: ITER: 71127 RES: 1.3004e-05
Jacobi stage: ITER: 71455 RES: 1.2365e-05
Jacobi stage: ITER: 71783 RES: 1.1757e-05
Jacobi stage: ITER: 72113 RES: 1.1173e-05
Jacobi stage: ITER: 72443 RES: 1.0624e-05
Jacobi stage: ITER: 72771 RES: 1.0102e-05
Jacobi stage: ITER: 73099 RES: 9.6057e-06
Jacobi stage: ITER: 73429 RES: 9.1281e-06
Jacobi stage: ITER: 73759 RES: 8.6796e-06
Jacobi stage: ITER: 74087 RES: 8.2531e-06
Jacobi stage: ITER: 74417 RES: 7.8428e-06
Jacobi stage: ITER: 74745 RES: 7.4575e-06
Jacobi stage: ITER: 75075 RES: 7.0911e-06
Jacobi stage: ITER: 75403 RES: 6.7427e-06
Jacobi stage: ITER: 75731 RES: 6.4114e-06
Jacobi stage: ITER: 76061 RES: 6.0926e-06
Jacobi stage: ITER: 76389 RES: 5.7933e-06
Jacobi stage: ITER: 76719 RES: 5.5086e-06
Jacobi stage: ITER: 77047 RES: 5.238e-06
Jacobi stage: ITER: 77375 RES: 4.9806e-06
Jacobi stage: ITER: 77705 RES: 4.733e-06
Jacobi stage: ITER: 78033 RES: 4.5004e-06
Jacobi stage: ITER: 78363 RES: 4.2793e-06
Jacobi stage: ITER: 78691 RES: 4.0691e-06
Jacobi stage: ITER: 79019 RES: 3.8691e-06
Jacobi stage: ITER: 79349 RES: 3.6768e-06
Jacobi stage: ITER: 79677 RES: 3.4961e-06
Jacobi stage: ITER: 80007 RES: 3.3243e-06
Jacobi stage: ITER: 80335 RES: 3.161e-06
Jacobi stage: ITER: 80663 RES: 3.0057e-06
Jacobi stage: ITER: 80993 RES: 2.8563e-06
Jacobi stage: ITER: 81323 RES: 2.7159e-06
Jacobi stage: ITER: 81651 RES: 2.5825e-06
Jacobi stage: ITER: 81979 RES: 2.4556e-06
Jacobi stage: ITER: 82309 RES: 2.3335e-06
Jacobi stage: ITER: 82639 RES: 2.2189e-06
Jacobi stage: ITER: 82967 RES: 2.1098e-06
Jacobi stage: ITER: 83295 RES: 2.0062e-06
Jacobi stage: ITER: 83625 RES: 1.9064e-06
Jacobi stage: ITER: 83953 RES: 1.8128e-06
Jacobi stage: ITER: 84283 RES: 1.7237e-06
Jacobi stage: ITER: 84611 RES: 1.639e-06
Jacobi stage: ITER: 84939 RES: 1.5585e-06
Jacobi stage: ITER: 85269 RES: 1.481e-06
Jacobi stage: ITER: 85599 RES: 1.4082e-06
Jacobi stage: ITER: 85927 RES: 1.339e-06
Jacobi stage: ITER: 86255 RES: 1.2732e-06
Jacobi stage: ITER: 86583 RES: 1.2107e-06
Jacobi stage: ITER: 86913 RES: 1.1505e-06
Jacobi stage: ITER: 87243 RES: 1.094e-06
Jacobi stage: ITER: 87571 RES: 1.0402e-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.7405e-
ELA: 0.10011 Jacobi stage: ITER: 30387 RES: 0.0067843
ELA: 0.20021 Jacobi stage: ITER: 60743 RES: 6.4088e-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.3695e-
ELA: 0.10012 Jacobi stage: ITER: 381 RES: 0.88222
ELA: 0.20023 Jacobi stage: ITER: 709 RES: 0.79978
ELA: 0.30034 Jacobi stage: ITER: 1035 RES: 0.73171
ELA: 0.40046 Jacobi stage: ITER: 1361 RES: 0.67335
ELA: 0.50057 Jacobi stage: ITER: 1687 RES: 0.62389
ELA: 0.60068 Jacobi stage: ITER: 2015 RES: 0.58037
ELA: 0.70079 Jacobi stage: ITER: 2339 RES: 0.54257
ELA: 0.8009 Jacobi stage: ITER: 2667 RES: 0.50852
ELA: 0.90101 Jacobi stage: ITER: 2991 RES: 0.4783
ELA: 1.0011 Jacobi stage: ITER: 3319 RES: 0.45055
ELA: 1.1012 Jacobi stage: ITER: 3663 RES: 0.42402
ELA: 1.2013 Jacobi stage: ITER: 3987 RES: 0.40106
ELA: 1.3014 Jacobi stage: ITER: 4315 RES: 0.37953
ELA: 1.4015 Jacobi stage: ITER: 4641 RES: 0.3595
ELA: 1.5016 Jacobi stage: ITER: 4967 RES: 0.34101
ELA: 1.6017 Jacobi stage: ITER: 5281 RES: 0.32408
ELA: 1.7018 Jacobi stage: ITER: 5611 RES: 0.30755
ELA: 1.8019 Jacobi stage: ITER: 5941 RES: 0.29179
ELA: 1.902 Jacobi stage: ITER: 6271 RES: 0.27711
ELA: 2.0021 Jacobi stage: ITER: 6547 RES: 0.26538
ELA: 2.1021 Jacobi stage: ITER: 6855 RES: 0.25292
ELA: 2.2022 Jacobi stage: ITER: 7187 RES: 0.24018
ELA: 2.3023 Jacobi stage: ITER: 7517 RES: 0.22812
ELA: 2.4024 Jacobi stage: ITER: 7849 RES: 0.21669
ELA: 2.5025 Jacobi stage: ITER: 8175 RES: 0.2061
ELA: 2.6026 Jacobi stage: ITER: 8483 RES: 0.19652
ELA: 2.7027 Jacobi stage: ITER: 8811 RES: 0.18683
ELA: 2.8029 Jacobi stage: ITER: 9139 RES: 0.17762
ELA: 2.903 Jacobi stage: ITER: 9471 RES: 0.16876
ELA: 3.0031 Jacobi stage: ITER: 9793 RES: 0.16055
ELA: 3.1031 Jacobi stage: ITER: 10103 RES: 0.15312
ELA: 3.2032 Jacobi stage: ITER: 10427 RES: 0.14568
ELA: 3.3033 Jacobi stage: ITER: 10759 RES: 0.13843
ELA: 3.4034 Jacobi stage: ITER: 11091 RES: 0.13154
ELA: 3.5035 Jacobi stage: ITER: 11401 RES: 0.12538
ELA: 3.6035 Jacobi stage: ITER: 11715 RES: 0.11951
ELA: 3.7036 Jacobi stage: ITER: 12041 RES: 0.11363
ELA: 3.8037 Jacobi stage: ITER: 12375 RES: 0.10798
ELA: 3.9038 Jacobi stage: ITER: 12705 RES: 0.10261
ELA: 4.0039 Jacobi stage: ITER: 13015 RES: 0.09787
ELA: 4.104 Jacobi stage: ITER: 13331 RES: 0.093232
ELA: 4.204 Jacobi stage: ITER: 13659 RES: 0.088651
ELA: 4.3041 Jacobi stage: ITER: 13993 RES: 0.084191
ELA: 4.4042 Jacobi stage: ITER: 14303 RES: 0.080301
ELA: 4.5043 Jacobi stage: ITER: 14613 RES: 0.076543
ELA: 4.6044 Jacobi stage: ITER: 14941 RES: 0.072782
ELA: 4.7045 Jacobi stage: ITER: 15269 RES: 0.069206
ELA: 4.8046 Jacobi stage: ITER: 15603 RES: 0.065765
ELA: 4.9047 Jacobi stage: ITER: 15937 RES: 0.062457
ELA: 5.0048 Jacobi stage: ITER: 16235 RES: 0.05968
ELA: 5.1049 Jacobi stage: ITER: 16559 RES: 0.056783
ELA: 5.205 Jacobi stage: ITER: 16891 RES: 0.05396
ELA: 5.305 Jacobi stage: ITER: 17223 RES: 0.051277
ELA: 5.4051 Jacobi stage: ITER: 17543 RES: 0.048818
ELA: 5.5052 Jacobi stage: ITER: 17851 RES: 0.046562
ELA: 5.6053 Jacobi stage: ITER: 18179 RES: 0.044274
ELA: 5.7053 Jacobi stage: ITER: 18511 RES: 0.042073
ELA: 5.8054 Jacobi stage: ITER: 18843 RES: 0.039981
ELA: 5.9055 Jacobi stage: ITER: 19173 RES: 0.037993
ELA: 6.0057 Jacobi stage: ITER: 19465 RES: 0.036327
ELA: 6.1057 Jacobi stage: ITER: 19793 RES: 0.034542
ELA: 6.2058 Jacobi stage: ITER: 20127 RES: 0.032825
ELA: 6.3059 Jacobi stage: ITER: 20461 RES: 0.031174
ELA: 6.406 Jacobi stage: ITER: 20761 RES: 0.02977
ELA: 6.5061 Jacobi stage: ITER: 21087 RES: 0.028324
ELA: 6.6062 Jacobi stage: ITER: 21419 RES: 0.026916
ELA: 6.7063 Jacobi stage: ITER: 21747 RES: 0.025594
ELA: 6.8064 Jacobi stage: ITER: 22079 RES: 0.024321
ELA: 6.9065 Jacobi stage: ITER: 22413 RES: 0.023098
ELA: 7.0065 Jacobi stage: ITER: 22721 RES: 0.022031
ELA: 7.1066 Jacobi stage: ITER: 23023 RES: 0.021038
ELA: 7.2067 Jacobi stage: ITER: 23351 RES: 0.020005
ELA: 7.3068 Jacobi stage: ITER: 23679 RES: 0.019022
ELA: 7.4069 Jacobi stage: ITER: 24011 RES: 0.018076
ELA: 7.507 Jacobi stage: ITER: 24345 RES: 0.017167
ELA: 7.6071 Jacobi stage: ITER: 24651 RES: 0.016384
ELA: 7.7072 Jacobi stage: ITER: 24959 RES: 0.015627
ELA: 7.8073 Jacobi stage: ITER: 25289 RES: 0.01485
ELA: 7.9073 Jacobi stage: ITER: 25621 RES: 0.014111
ELA: 8.0074 Jacobi stage: ITER: 25953 RES: 0.01341
ELA: 8.1075 Jacobi stage: ITER: 26259 RES: 0.012798
ELA: 8.2076 Jacobi stage: ITER: 26519 RES: 0.012297
ELA: 8.3077 Jacobi stage: ITER: 26839 RES: 0.011707
ELA: 8.4078 Jacobi stage: ITER: 27173 RES: 0.011118
ELA: 8.5079 Jacobi stage: ITER: 27507 RES: 0.010566
ELA: 8.608 Jacobi stage: ITER: 27821 RES: 0.010065
ELA: 8.708 Jacobi stage: ITER: 28151 RES: 0.0095705
ELA: 8.8081 Jacobi stage: ITER: 28459 RES: 0.0091282
ELA: 8.9082 Jacobi stage: ITER: 28775 RES: 0.0086958
ELA: 9.0083 Jacobi stage: ITER: 29103 RES: 0.0082685
ELA: 9.1084 Jacobi stage: ITER: 29433 RES: 0.0078574
ELA: 9.2085 Jacobi stage: ITER: 29763 RES: 0.0074714
ELA: 9.3085 Jacobi stage: ITER: 30091 RES: 0.0071043
ELA: 9.4086 Jacobi stage: ITER: 30419 RES: 0.0067552
ELA: 9.5087 Jacobi stage: ITER: 30749 RES: 0.0064194
ELA: 9.6088 Jacobi stage: ITER: 31077 RES: 0.006104
ELA: 9.7089 Jacobi stage: ITER: 31407 RES: 0.005804
ELA: 9.8089 Jacobi stage: ITER: 31735 RES: 0.0055189
ELA: 9.909 Jacobi stage: ITER: 32063 RES: 0.0052477
ELA: 10.009 Jacobi stage: ITER: 32393 RES: 0.0049868
ELA: 10.109 Jacobi stage: ITER: 32723 RES: 0.0047418
ELA: 10.209 Jacobi stage: ITER: 33051 RES: 0.0045088
ELA: 10.309 Jacobi stage: ITER: 33379 RES: 0.0042873
ELA: 10.409 Jacobi stage: ITER: 33709 RES: 0.0040741
ELA: 10.509 Jacobi stage: ITER: 34039 RES: 0.003874
ELA: 10.61 Jacobi stage: ITER: 34367 RES: 0.0036836
ELA: 10.71 Jacobi stage: ITER: 34695 RES: 0.0035026
ELA: 10.81 Jacobi stage: ITER: 35025 RES: 0.0033285
ELA: 10.91 Jacobi stage: ITER: 35353 RES: 0.0031649
ELA: 11.01 Jacobi stage: ITER: 35683 RES: 0.0030094
ELA: 11.11 Jacobi stage: ITER: 36011 RES: 0.0028616
ELA: 11.21 Jacobi stage: ITER: 36339 RES: 0.002721
ELA: 11.31 Jacobi stage: ITER: 36669 RES: 0.0025857
ELA: 11.41 Jacobi stage: ITER: 36999 RES: 0.0024587
ELA: 11.51 Jacobi stage: ITER: 37327 RES: 0.0023379
ELA: 11.61 Jacobi stage: ITER: 37655 RES: 0.002223
ELA: 11.711 Jacobi stage: ITER: 37985 RES: 0.0021125
ELA: 11.811 Jacobi stage: ITER: 38315 RES: 0.0020087
ELA: 11.911 Jacobi stage: ITER: 38643 RES: 0.00191
ELA: 12.011 Jacobi stage: ITER: 38971 RES: 0.0018161
ELA: 12.111 Jacobi stage: ITER: 39301 RES: 0.0017258
ELA: 12.211 Jacobi stage: ITER: 39629 RES: 0.001641
ELA: 12.311 Jacobi stage: ITER: 39959 RES: 0.0015604
ELA: 12.412 Jacobi stage: ITER: 40287 RES: 0.0014838
ELA: 12.512 Jacobi stage: ITER: 40617 RES: 0.00141
ELA: 12.612 Jacobi stage: ITER: 40947 RES: 0.0013407
ELA: 12.712 Jacobi stage: ITER: 41275 RES: 0.0012748
ELA: 12.812 Jacobi stage: ITER: 41603 RES: 0.0012122
ELA: 12.912 Jacobi stage: ITER: 41933 RES: 0.0011519
ELA: 13.012 Jacobi stage: ITER: 42261 RES: 0.0010953
ELA: 13.112 Jacobi stage: ITER: 42591 RES: 0.0010415
ELA: 13.212 Jacobi stage: ITER: 42919 RES: 0.00099034
ELA: 13.313 Jacobi stage: ITER: 43249 RES: 0.0009411
ELA: 13.413 Jacobi stage: ITER: 43577 RES: 0.00089486
ELA: 13.513 Jacobi stage: ITER: 43907 RES: 0.0008509
ELA: 13.613 Jacobi stage: ITER: 44235 RES: 0.00080909
ELA: 13.713 Jacobi stage: ITER: 44563 RES: 0.00076934
ELA: 13.813 Jacobi stage: ITER: 44891 RES: 0.00073154
ELA: 13.913 Jacobi stage: ITER: 45221 RES: 0.00069517
ELA: 14.013 Jacobi stage: ITER: 45551 RES: 0.00066101
ELA: 14.113 Jacobi stage: ITER: 45879 RES: 0.00062853
ELA: 14.213 Jacobi stage: ITER: 46207 RES: 0.00059765
ELA: 14.313 Jacobi stage: ITER: 46535 RES: 0.00056829
ELA: 14.413 Jacobi stage: ITER: 46865 RES: 0.00054003
ELA: 14.513 Jacobi stage: ITER: 47195 RES: 0.0005135
ELA: 14.614 Jacobi stage: ITER: 47523 RES: 0.00048827
ELA: 14.714 Jacobi stage: ITER: 47851 RES: 0.00046428
ELA: 14.814 Jacobi stage: ITER: 48179 RES: 0.00044147
ELA: 14.914 Jacobi stage: ITER: 48509 RES: 0.00041952
ELA: 15.014 Jacobi stage: ITER: 48839 RES: 0.00039891
ELA: 15.114 Jacobi stage: ITER: 49167 RES: 0.00037931
ELA: 15.214 Jacobi stage: ITER: 49495 RES: 0.00036067
ELA: 15.314 Jacobi stage: ITER: 49823 RES: 0.00034295
ELA: 15.414 Jacobi stage: ITER: 50153 RES: 0.0003259
ELA: 15.514 Jacobi stage: ITER: 50483 RES: 0.00030989
ELA: 15.614 Jacobi stage: ITER: 50811 RES: 0.00029466
ELA: 15.714 Jacobi stage: ITER: 51139 RES: 0.00028018
ELA: 15.814 Jacobi stage: ITER: 51467 RES: 0.00026642
ELA: 15.914 Jacobi stage: ITER: 51797 RES: 0.00025317
ELA: 16.015 Jacobi stage: ITER: 52125 RES: 0.00024073
ELA: 16.115 Jacobi stage: ITER: 52455 RES: 0.0002289
ELA: 16.215 Jacobi stage: ITER: 52783 RES: 0.00021766
ELA: 16.315 Jacobi stage: ITER: 53111 RES: 0.00020696
ELA: 16.415 Jacobi stage: ITER: 53441 RES: 0.00019667
ELA: 16.515 Jacobi stage: ITER: 53769 RES: 0.00018701
ELA: 16.615 Jacobi stage: ITER: 54099 RES: 0.00017782
ELA: 16.715 Jacobi stage: ITER: 54427 RES: 0.00016908
ELA: 16.815 Jacobi stage: ITER: 54755 RES: 0.00016078
ELA: 16.915 Jacobi stage: ITER: 55085 RES: 0.00015278
ELA: 17.015 Jacobi stage: ITER: 55415 RES: 0.00014528
ELA: 17.115 Jacobi stage: ITER: 55743 RES: 0.00013814
ELA: 17.215 Jacobi stage: ITER: 56071 RES: 0.00013135
ELA: 17.315 Jacobi stage: ITER: 56399 RES: 0.0001249
ELA: 17.416 Jacobi stage: ITER: 56729 RES: 0.00011869
ELA: 17.516 Jacobi stage: ITER: 57057 RES: 0.00011286
ELA: 17.616 Jacobi stage: ITER: 57387 RES: 0.00010731
ELA: 17.716 Jacobi stage: ITER: 57715 RES: 0.00010204
ELA: 17.816 Jacobi stage: ITER: 58043 RES: 9.7026e-05
ELA: 17.916 Jacobi stage: ITER: 58373 RES: 9.2202e-05
ELA: 18.016 Jacobi stage: ITER: 58701 RES: 8.7672e-05
ELA: 18.116 Jacobi stage: ITER: 59031 RES: 8.3364e-05
ELA: 18.216 Jacobi stage: ITER: 59359 RES: 7.9268e-05
ELA: 18.316 Jacobi stage: ITER: 59687 RES: 7.5374e-05
ELA: 18.416 Jacobi stage: ITER: 60017 RES: 7.1626e-05
ELA: 18.516 Jacobi stage: ITER: 60345 RES: 6.8107e-05
ELA: 18.616 Jacobi stage: ITER: 60675 RES: 6.4761e-05
ELA: 18.716 Jacobi stage: ITER: 61003 RES: 6.1579e-05
ELA: 18.817 Jacobi stage: ITER: 61331 RES: 5.8553e-05
ELA: 18.917 Jacobi stage: ITER: 61661 RES: 5.5642e-05
ELA: 19.017 Jacobi stage: ITER: 61989 RES: 5.2908e-05
ELA: 19.117 Jacobi stage: ITER: 62319 RES: 5.0309e-05
ELA: 19.217 Jacobi stage: ITER: 62647 RES: 4.7837e-05
ELA: 19.317 Jacobi stage: ITER: 62975 RES: 4.5486e-05
ELA: 19.417 Jacobi stage: ITER: 63305 RES: 4.3225e-05
ELA: 19.517 Jacobi stage: ITER: 63633 RES: 4.1101e-05
ELA: 19.617 Jacobi stage: ITER: 63963 RES: 3.9082e-05
ELA: 19.717 Jacobi stage: ITER: 64291 RES: 3.7162e-05
ELA: 19.817 Jacobi stage: ITER: 64619 RES: 3.5336e-05
ELA: 19.917 Jacobi stage: ITER: 64949 RES: 3.3579e-05
ELA: 20.018 Jacobi stage: ITER: 65279 RES: 3.1929e-05
ELA: 20.118 Jacobi stage: ITER: 65607 RES: 3.036e-05
ELA: 20.218 Jacobi stage: ITER: 65935 RES: 2.8869e-05
ELA: 20.318 Jacobi stage: ITER: 66265 RES: 2.7433e-05
ELA: 20.418 Jacobi stage: ITER: 66595 RES: 2.6085e-05
ELA: 20.518 Jacobi stage: ITER: 66923 RES: 2.4804e-05
ELA: 20.618 Jacobi stage: ITER: 67251 RES: 2.3585e-05
ELA: 20.718 Jacobi stage: ITER: 67581 RES: 2.2412e-05
ELA: 20.819 Jacobi stage: ITER: 67909 RES: 2.1311e-05
ELA: 20.919 Jacobi stage: ITER: 68239 RES: 2.0264e-05
ELA: 21.019 Jacobi stage: ITER: 68567 RES: 1.9269e-05
ELA: 21.119 Jacobi stage: ITER: 68895 RES: 1.8322e-05
ELA: 21.219 Jacobi stage: ITER: 69225 RES: 1.7411e-05
ELA: 21.319 Jacobi stage: ITER: 69553 RES: 1.6555e-05
ELA: 21.419 Jacobi stage: ITER: 69883 RES: 1.5742e-05
ELA: 21.519 Jacobi stage: ITER: 70211 RES: 1.4969e-05
ELA: 21.619 Jacobi stage: ITER: 70539 RES: 1.4233e-05
ELA: 21.719 Jacobi stage: ITER: 70869 RES: 1.3525e-05
ELA: 21.819 Jacobi stage: ITER: 71199 RES: 1.2861e-05
ELA: 21.92 Jacobi stage: ITER: 71527 RES: 1.2229e-05
ELA: 22.02 Jacobi stage: ITER: 71855 RES: 1.1628e-05
ELA: 22.12 Jacobi stage: ITER: 72185 RES: 1.105e-05
ELA: 22.22 Jacobi stage: ITER: 72515 RES: 1.0507e-05
ELA: 22.32 Jacobi stage: ITER: 72843 RES: 9.9909e-06
ELA: 22.42 Jacobi stage: ITER: 73171 RES: 9.5e-06
ELA: 22.52 Jacobi stage: ITER: 73501 RES: 9.0277e-06
ELA: 22.62 Jacobi stage: ITER: 73831 RES: 8.5841e-06
ELA: 22.721 Jacobi stage: ITER: 74159 RES: 8.1624e-06
ELA: 22.821 Jacobi stage: ITER: 74487 RES: 7.7613e-06
ELA: 22.921 Jacobi stage: ITER: 74817 RES: 7.3754e-06
ELA: 23.021 Jacobi stage: ITER: 75145 RES: 7.0131e-06
ELA: 23.121 Jacobi stage: ITER: 75475 RES: 6.6685e-06
ELA: 23.221 Jacobi stage: ITER: 75803 RES: 6.3408e-06
ELA: 23.321 Jacobi stage: ITER: 76133 RES: 6.0256e-06
ELA: 23.421 Jacobi stage: ITER: 76461 RES: 5.7295e-06
ELA: 23.521 Jacobi stage: ITER: 76791 RES: 5.448e-06
ELA: 23.621 Jacobi stage: ITER: 77119 RES: 5.1804e-06
ELA: 23.722 Jacobi stage: ITER: 77447 RES: 4.9258e-06
ELA: 23.822 Jacobi stage: ITER: 77775 RES: 4.6838e-06
ELA: 23.922 Jacobi stage: ITER: 78105 RES: 4.4509e-06
ELA: 24.022 Jacobi stage: ITER: 78543 RES: 4.1626e-06
ELA: 24.122 Jacobi stage: ITER: 79039 RES: 3.8573e-06
ELA: 24.222 Jacobi stage: ITER: 79533 RES: 3.5743e-06
ELA: 24.322 Jacobi stage: ITER: 80027 RES: 3.3141e-06
ELA: 24.422 Jacobi stage: ITER: 80521 RES: 3.071e-06
ELA: 24.522 Jacobi stage: ITER: 81015 RES: 2.8475e-06
ELA: 24.622 Jacobi stage: ITER: 81511 RES: 2.6386e-06
ELA: 24.723 Jacobi stage: ITER: 82003 RES: 2.4466e-06
ELA: 24.823 Jacobi stage: ITER: 82499 RES: 2.2671e-06
ELA: 24.923 Jacobi stage: ITER: 82993 RES: 2.1008e-06
ELA: 25.023 Jacobi stage: ITER: 83487 RES: 1.9479e-06
ELA: 25.123 Jacobi stage: ITER: 83983 RES: 1.805e-06
ELA: 25.223 Jacobi stage: ITER: 84477 RES: 1.6726e-06
ELA: 25.323 Jacobi stage: ITER: 84971 RES: 1.5508e-06
ELA: 25.423 Jacobi stage: ITER: 85465 RES: 1.4371e-06
ELA: 25.524 Jacobi stage: ITER: 85979 RES: 1.3284e-06
ELA: 25.624 Jacobi stage: ITER: 86473 RES: 1.2309e-06
ELA: 25.724 Jacobi stage: ITER: 86967 RES: 1.1413e-06
ELA: 25.824 Jacobi stage: ITER: 87461 RES: 1.0576e-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 Thu Dec 5 2024 15:23:19 for Template Numerical Library by doxygen 1.12.0