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

Table of Contents

Introduction

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

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

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

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

Iterative solvers of linear systems

Basic setup

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

1#include <iostream>
2#include <memory>
3#include <TNL/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 iterativeLinearSolverExample()
10{
11 /***
12 * Set the following matrix (dots represent zero matrix elements):
13 *
14 * / 2.5 -1 . . . \
15 * | -1 2.5 -1 . . |
16 * | . -1 2.5 -1. . |
17 * | . . -1 2.5 -1 |
18 * \ . . . -1 2.5 /
19 */
20 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
21 using Vector = TNL::Containers::Vector< double, Device >;
22 const int size( 5 );
23 auto matrix_ptr = std::make_shared< MatrixType >();
24 matrix_ptr->setDimensions( size, size );
25 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
26
27 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
28 const int rowIdx = row.getRowIndex();
29 if( rowIdx == 0 )
30 {
31 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
32 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
33 }
34 else if( rowIdx == size - 1 )
35 {
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 {
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.
64 */
65 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
66 LinearSolver solver;
67 solver.setMatrix( matrix_ptr );
68 solver.setConvergenceResidue( 1.0e-6 );
69 solver.solve( b, x );
70 std::cout << "Vector x = " << x << std::endl;
71}
72
73int main( 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: CudaCallable.h:22
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition: Vector.h:40
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:127
TNL::Solvers::Linear::TFQMR
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition: TFQMR.h:26
std::endl
T endl(T... args)

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

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

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

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

The result looks as follows:

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

Setup with a solver monitor

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

1#include <iostream>
2#include <memory>
3#include <chrono>
4#include <thread>
5#include <TNL/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 iterativeLinearSolverExample()
12{
13 /***
14 * Set the following matrix (dots represent zero matrix elements):
15 *
16 * / 2.5 -1 . . . \
17 * | -1 2.5 -1 . . |
18 * | . -1 2.5 -1. . |
19 * | . . -1 2.5 -1 |
20 * \ . . . -1 2.5 /
21 */
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 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 )
32 {
33 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
34 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
35 }
36 else if( rowIdx == size - 1 )
37 {
38 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
39 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
40 }
41 else
42 {
43 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
44 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
45 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
46 }
47 };
48
49 /***
50 * Set the matrix elements.
51 */
52 matrix_ptr->forAllRows( f );
53 std::cout << *matrix_ptr << std::endl;
54
55 /***
56 * Set the right-hand side vector.
57 */
58 Vector x( size, 1.0 );
59 Vector b( size );
60 matrix_ptr->vectorProduct( x, b );
61 x = 0.0;
62 std::cout << "Vector b = " << b << std::endl;
63
64 /***
65 * 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 monitorThread(monitor);
78 monitor.setRefreshRate(10); // refresh rate in milliseconds
79 monitor.setVerbose(1);
80 monitor.setStage( "Jacobi stage:" );
81 solver.setSolverMonitor(monitor);
82 solver.setConvergenceResidue( 1.0e-6 );
83 solver.solve( b, x );
84 monitor.stopMainLoop();
85 std::cout << "Vector x = " << x << std::endl;
86}
87
88int main( 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:26
TNL::Solvers::SolverMonitorThread
A RAII wrapper for launching the SolverMonitor's main loop in a separate thread.
Definition: SolverMonitor.h:137

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

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Jacobi stage: ITER: 26437 RES: 0.012449
Jacobi stage: ITER: 66687 RES: 2.5719e-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: 5579 RES: 0.30912
Jacobi stage: ITER: 9419 RES: 0.17012
Jacobi stage: ITER: 10227 RES: 0.15023
Jacobi stage: ITER: 11115 RES: 0.13105
Jacobi stage: ITER: 12019 RES: 0.11405
Jacobi stage: ITER: 12937 RES: 0.09902
Jacobi stage: ITER: 13851 RES: 0.086074
Jacobi stage: ITER: 14757 RES: 0.074868
Jacobi stage: ITER: 15663 RES: 0.065161
Jacobi stage: ITER: 16569 RES: 0.056678
Jacobi stage: ITER: 17475 RES: 0.04933
Jacobi stage: ITER: 18383 RES: 0.042908
Jacobi stage: ITER: 19289 RES: 0.037322
Jacobi stage: ITER: 20197 RES: 0.032464
Jacobi stage: ITER: 21103 RES: 0.028255
Jacobi stage: ITER: 22011 RES: 0.024577
Jacobi stage: ITER: 22919 RES: 0.021377
Jacobi stage: ITER: 23825 RES: 0.018594
Jacobi stage: ITER: 24733 RES: 0.016174
Jacobi stage: ITER: 25641 RES: 0.014068
Jacobi stage: ITER: 26547 RES: 0.012244
Jacobi stage: ITER: 27455 RES: 0.01065
Jacobi stage: ITER: 28361 RES: 0.0092638
Jacobi stage: ITER: 29269 RES: 0.0080579
Jacobi stage: ITER: 30175 RES: 0.0070132
Jacobi stage: ITER: 31083 RES: 0.0061002
Jacobi stage: ITER: 31991 RES: 0.0053061
Jacobi stage: ITER: 32895 RES: 0.0046182
Jacobi stage: ITER: 33803 RES: 0.004017
Jacobi stage: ITER: 34711 RES: 0.003494
Jacobi stage: ITER: 35617 RES: 0.0030392
Jacobi stage: ITER: 36523 RES: 0.0026452
Jacobi stage: ITER: 37431 RES: 0.0023008
Jacobi stage: ITER: 38337 RES: 0.0020013
Jacobi stage: ITER: 39245 RES: 0.0017408
Jacobi stage: ITER: 40151 RES: 0.0015151
Jacobi stage: ITER: 41059 RES: 0.0013178
Jacobi stage: ITER: 41967 RES: 0.0011463
Jacobi stage: ITER: 42873 RES: 0.00099706
Jacobi stage: ITER: 43779 RES: 0.00086779
Jacobi stage: ITER: 44685 RES: 0.00075482
Jacobi stage: ITER: 45593 RES: 0.00065656
Jacobi stage: ITER: 46499 RES: 0.00057144
Jacobi stage: ITER: 47405 RES: 0.00049705
Jacobi stage: ITER: 48311 RES: 0.00043261
Jacobi stage: ITER: 49215 RES: 0.00037652
Jacobi stage: ITER: 49955 RES: 0.00033607
Jacobi stage: ITER: 50861 RES: 0.00029232
Jacobi stage: ITER: 51775 RES: 0.00025411
Jacobi stage: ITER: 52689 RES: 0.00022076
Jacobi stage: ITER: 53603 RES: 0.0001919
Jacobi stage: ITER: 54515 RES: 0.00016681
Jacobi stage: ITER: 55427 RES: 0.00014501
Jacobi stage: ITER: 56339 RES: 0.00012605
Jacobi stage: ITER: 57251 RES: 0.00010958
Jacobi stage: ITER: 58163 RES: 9.5254e-05
Jacobi stage: ITER: 59075 RES: 8.2803e-05
Jacobi stage: ITER: 59987 RES: 7.1979e-05
Jacobi stage: ITER: 60899 RES: 6.257e-05
Jacobi stage: ITER: 61819 RES: 5.4325e-05
Jacobi stage: ITER: 62731 RES: 4.7224e-05
Jacobi stage: ITER: 63645 RES: 4.1026e-05
Jacobi stage: ITER: 64557 RES: 3.5663e-05
Jacobi stage: ITER: 65469 RES: 3.1001e-05
Jacobi stage: ITER: 66379 RES: 2.6965e-05
Jacobi stage: ITER: 67287 RES: 2.3455e-05
Jacobi stage: ITER: 68195 RES: 2.0402e-05
Jacobi stage: ITER: 69101 RES: 1.7746e-05
Jacobi stage: ITER: 70009 RES: 1.5436e-05
Jacobi stage: ITER: 70917 RES: 1.3426e-05
Jacobi stage: ITER: 71823 RES: 1.1685e-05
Jacobi stage: ITER: 72731 RES: 1.0164e-05
Jacobi stage: ITER: 73637 RES: 8.8411e-06
Jacobi stage: ITER: 74545 RES: 7.6901e-06
Jacobi stage: ITER: 75453 RES: 6.689e-06
Jacobi stage: ITER: 76359 RES: 5.8218e-06
Jacobi stage: ITER: 77267 RES: 5.0639e-06
Jacobi stage: ITER: 78175 RES: 4.4047e-06
Jacobi stage: ITER: 79087 RES: 3.8289e-06
Jacobi stage: ITER: 79999 RES: 3.3284e-06
Jacobi stage: ITER: 80911 RES: 2.8934e-06
Jacobi stage: ITER: 81823 RES: 2.5151e-06
Jacobi stage: ITER: 82735 RES: 2.1864e-06
Jacobi stage: ITER: 83639 RES: 1.9029e-06
Jacobi stage: ITER: 84551 RES: 1.6542e-06
Jacobi stage: ITER: 85463 RES: 1.438e-06
Jacobi stage: ITER: 86375 RES: 1.25e-06
Jacobi stage: ITER: 87293 RES: 1.0853e-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 iterativeLinearSolverExample()
13{
14 /***
15 * Set the following matrix (dots represent zero matrix elements):
16 *
17 * / 2.5 -1 . . . \
18 * | -1 2.5 -1 . . |
19 * | . -1 2.5 -1. . |
20 * | . . -1 2.5 -1 |
21 * \ . . . -1 2.5 /
22 */
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 const int rowIdx = row.getRowIndex();
32 if( rowIdx == 0 )
33 {
34 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
35 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
36 }
37 else if( rowIdx == size - 1 )
38 {
39 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
40 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
41 }
42 else
43 {
44 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
45 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
46 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
47 }
48 };
49
50 /***
51 * Set the matrix elements.
52 */
53 matrix_ptr->forAllRows( f );
54 std::cout << *matrix_ptr << std::endl;
55
56 /***
57 * Set the right-hand side vector.
58 */
59 Vector x( size, 1.0 );
60 Vector b( size );
61 matrix_ptr->vectorProduct( x, b );
62 x = 0.0;
63 std::cout << "Vector b = " << b << std::endl;
64
65 /***
66 * Setup solver of the linear system.
67 */
68 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
69 LinearSolver solver;
70 solver.setMatrix( matrix_ptr );
71 solver.setOmega( 0.0005 );
72
73 /***
74 * Setup monitor of the iterative solver.
75 */
76 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
77 IterativeSolverMonitorType monitor;
78 TNL::Solvers::SolverMonitorThread mmonitorThread(monitor);
79 monitor.setRefreshRate(10); // refresh rate in milliseconds
80 monitor.setVerbose(1);
81 monitor.setStage( "Jacobi stage:" );
82 TNL::Timer timer;
83 monitor.setTimer( timer );
84 timer.start();
85 solver.setSolverMonitor(monitor);
86 solver.setConvergenceResidue( 1.0e-6 );
87 solver.solve( b, x );
88 monitor.stopMainLoop();
89 std::cout << "Vector x = " << x << std::endl;
90}
91
92int main( 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:28
TNL::Timer::start
void start()
Starts timer.
Definition: Timer.hpp:49

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

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
ELA:3.16e-06
ELA: 0.10034 Jacobi stage: ITER: 40421 RES: 0.0014531
ELA: 0.20044 Jacobi stage: ITER: 80941 RES: 2.8792e-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:9.34e-06
ELA: 0.1002 Jacobi stage: ITER: 807 RES: 0.77838
ELA: 0.2003 Jacobi stage: ITER: 1663 RES: 0.62731
ELA: 0.30039 Jacobi stage: ITER: 2555 RES: 0.51972
ELA: 0.40049 Jacobi stage: ITER: 3475 RES: 0.43822
ELA: 0.50059 Jacobi stage: ITER: 4393 RES: 0.37451
ELA: 0.60069 Jacobi stage: ITER: 5299 RES: 0.32325
ELA: 0.70079 Jacobi stage: ITER: 6205 RES: 0.27991
ELA: 0.80089 Jacobi stage: ITER: 7111 RES: 0.24304
ELA: 0.90099 Jacobi stage: ITER: 8017 RES: 0.21113
ELA: 1.0011 Jacobi stage: ITER: 8923 RES: 0.18363
ELA: 1.1012 Jacobi stage: ITER: 9831 RES: 0.15967
ELA: 1.2013 Jacobi stage: ITER: 10739 RES: 0.13885
ELA: 1.3014 Jacobi stage: ITER: 11645 RES: 0.12076
ELA: 1.4015 Jacobi stage: ITER: 12551 RES: 0.1051
ELA: 1.5016 Jacobi stage: ITER: 13459 RES: 0.091417
ELA: 1.6017 Jacobi stage: ITER: 14367 RES: 0.079515
ELA: 1.7018 Jacobi stage: ITER: 15275 RES: 0.069163
ELA: 1.8019 Jacobi stage: ITER: 16179 RES: 0.060196
ELA: 1.902 Jacobi stage: ITER: 17087 RES: 0.05236
ELA: 2.0021 Jacobi stage: ITER: 17995 RES: 0.045543
ELA: 2.1022 Jacobi stage: ITER: 18901 RES: 0.039614
ELA: 2.2023 Jacobi stage: ITER: 19807 RES: 0.034457
ELA: 2.3024 Jacobi stage: ITER: 20715 RES: 0.02999
ELA: 2.4025 Jacobi stage: ITER: 21621 RES: 0.026086
ELA: 2.5026 Jacobi stage: ITER: 22603 RES: 0.02244
ELA: 2.6027 Jacobi stage: ITER: 23511 RES: 0.019519
ELA: 2.7028 Jacobi stage: ITER: 24417 RES: 0.016978
ELA: 2.8029 Jacobi stage: ITER: 25323 RES: 0.014777
ELA: 2.903 Jacobi stage: ITER: 26231 RES: 0.012853
ELA: 3.003 Jacobi stage: ITER: 27137 RES: 0.01118
ELA: 3.1031 Jacobi stage: ITER: 28043 RES: 0.0097306
ELA: 3.2032 Jacobi stage: ITER: 28951 RES: 0.0084638
ELA: 3.3033 Jacobi stage: ITER: 29857 RES: 0.007362
ELA: 3.4034 Jacobi stage: ITER: 30763 RES: 0.0064075
ELA: 3.5035 Jacobi stage: ITER: 31671 RES: 0.0055734
ELA: 3.6037 Jacobi stage: ITER: 32589 RES: 0.0048389
ELA: 3.7038 Jacobi stage: ITER: 33495 RES: 0.0042116
ELA: 3.8039 Jacobi stage: ITER: 34403 RES: 0.0036633
ELA: 3.904 Jacobi stage: ITER: 35309 RES: 0.0031864
ELA: 4.0041 Jacobi stage: ITER: 36215 RES: 0.0027733
ELA: 4.1042 Jacobi stage: ITER: 37121 RES: 0.0024123
ELA: 4.2043 Jacobi stage: ITER: 38027 RES: 0.0020995
ELA: 4.3043 Jacobi stage: ITER: 38935 RES: 0.0018262
ELA: 4.4044 Jacobi stage: ITER: 39841 RES: 0.0015885
ELA: 4.5045 Jacobi stage: ITER: 40589 RES: 0.0014161
ELA: 4.6046 Jacobi stage: ITER: 41483 RES: 0.0012347
ELA: 4.7047 Jacobi stage: ITER: 42395 RES: 0.0010733
ELA: 4.8048 Jacobi stage: ITER: 43309 RES: 0.00093247
ELA: 4.9049 Jacobi stage: ITER: 44223 RES: 0.00081058
ELA: 5.005 Jacobi stage: ITER: 45135 RES: 0.00070463
ELA: 5.1051 Jacobi stage: ITER: 46047 RES: 0.00061252
ELA: 5.2052 Jacobi stage: ITER: 46959 RES: 0.00053245
ELA: 5.3053 Jacobi stage: ITER: 47871 RES: 0.00046285
ELA: 5.4054 Jacobi stage: ITER: 48783 RES: 0.00040235
ELA: 5.5054 Jacobi stage: ITER: 49697 RES: 0.00034954
ELA: 5.6055 Jacobi stage: ITER: 50609 RES: 0.00030385
ELA: 5.7056 Jacobi stage: ITER: 51521 RES: 0.00026414
ELA: 5.8057 Jacobi stage: ITER: 52429 RES: 0.00022975
ELA: 5.9058 Jacobi stage: ITER: 53341 RES: 0.00019972
ELA: 6.0059 Jacobi stage: ITER: 54255 RES: 0.00017361
ELA: 6.106 Jacobi stage: ITER: 55167 RES: 0.00015092
ELA: 6.2061 Jacobi stage: ITER: 56079 RES: 0.00013119
ELA: 6.3062 Jacobi stage: ITER: 56991 RES: 0.00011404
ELA: 6.4063 Jacobi stage: ITER: 57899 RES: 9.9196e-05
ELA: 6.5064 Jacobi stage: ITER: 58805 RES: 8.6282e-05
ELA: 6.6065 Jacobi stage: ITER: 59711 RES: 7.5096e-05
ELA: 6.7066 Jacobi stage: ITER: 60619 RES: 6.532e-05
ELA: 6.8067 Jacobi stage: ITER: 61527 RES: 5.6817e-05
ELA: 6.9068 Jacobi stage: ITER: 62435 RES: 4.942e-05
ELA: 7.0069 Jacobi stage: ITER: 63339 RES: 4.3013e-05
ELA: 7.107 Jacobi stage: ITER: 64247 RES: 3.7414e-05
ELA: 7.2071 Jacobi stage: ITER: 65155 RES: 3.2543e-05
ELA: 7.3072 Jacobi stage: ITER: 66063 RES: 2.8307e-05
ELA: 7.4073 Jacobi stage: ITER: 66969 RES: 2.4622e-05
ELA: 7.5074 Jacobi stage: ITER: 67875 RES: 2.1429e-05
ELA: 7.6075 Jacobi stage: ITER: 68783 RES: 1.864e-05
ELA: 7.7076 Jacobi stage: ITER: 69695 RES: 1.6203e-05
ELA: 7.8077 Jacobi stage: ITER: 70607 RES: 1.4085e-05
ELA: 7.9078 Jacobi stage: ITER: 71519 RES: 1.2244e-05
ELA: 8.0079 Jacobi stage: ITER: 72431 RES: 1.0644e-05
ELA: 8.108 Jacobi stage: ITER: 73343 RES: 9.2523e-06
ELA: 8.2081 Jacobi stage: ITER: 74303 RES: 7.9838e-06
ELA: 8.3082 Jacobi stage: ITER: 75215 RES: 6.9402e-06
ELA: 8.4083 Jacobi stage: ITER: 76127 RES: 6.033e-06
ELA: 8.5083 Jacobi stage: ITER: 77039 RES: 5.2444e-06
ELA: 8.6085 Jacobi stage: ITER: 77951 RES: 4.5589e-06
ELA: 8.7088 Jacobi stage: ITER: 78865 RES: 3.9605e-06
ELA: 8.8089 Jacobi stage: ITER: 83031 RES: 2.0892e-06
ELA: 8.909 Jacobi stage: ITER: 87267 RES: 1.0899e-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 iterativeLinearSolverExample()
11{
12 /***
13 * Set the following matrix (dots represent zero matrix elements):
14 *
15 * / 2.5 -1 . . . \
16 * | -1 2.5 -1 . . |
17 * | . -1 2.5 -1. . |
18 * | . . -1 2.5 -1 |
19 * \ . . . -1 2.5 /
20 */
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 const int rowIdx = row.getRowIndex();
30 if( rowIdx == 0 )
31 {
32 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
33 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
34 }
35 else if( rowIdx == size - 1 )
36 {
37 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else
41 {
42 row.setElement( 0, rowIdx-1, -1.0 ); // element below the diagonal
43 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
44 row.setElement( 2, rowIdx+1, -1.0 ); // element above the diagonal
45 }
46 };
47
48 /***
49 * Set the matrix elements.
50 */
51 matrix_ptr->forAllRows( f );
52 std::cout << *matrix_ptr << std::endl;
53
54 /***
55 * Set the right-hand side vector.
56 */
57 Vector x( size, 1.0 );
58 Vector b( size );
59 matrix_ptr->vectorProduct( x, b );
60 x = 0.0;
61 std::cout << "Vector b = " << b << std::endl;
62
63 /***
64 * Solve the linear system using diagonal (Jacobi) preconditioner.
65 */
66 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
67 using Preconditioner = TNL::Solvers::Linear::Preconditioners::Diagonal< MatrixType >;
68 auto preconditioner_ptr = std::make_shared< Preconditioner >();
69 preconditioner_ptr->update( matrix_ptr );
70 LinearSolver solver;
71 solver.setMatrix( matrix_ptr );
72 solver.setPreconditioner( preconditioner_ptr );
73 solver.setConvergenceResidue( 1.0e-6 );
74 solver.solve( b, x );
75 std::cout << "Vector x = " << x << std::endl;
76}
77
78int main( 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:29

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

The result looks as follows:

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

Choosing the solver and preconditioner type at runtime

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

1#include <iostream>
2#include <memory>
3#include <TNL/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 iterativeLinearSolverExample()
10{
11 /***
12 * Set the following matrix (dots represent zero matrix elements):
13 *
14 * / 2.5 -1 . . . \
15 * | -1 2.5 -1 . . |
16 * | . -1 2.5 -1. . |
17 * | . . -1 2.5 -1 |
18 * \ . . . -1 2.5 /
19 */
20 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
21 using Vector = TNL::Containers::Vector< double, Device >;
22 const int size( 5 );
23 auto matrix_ptr = std::make_shared< MatrixType >();
24 matrix_ptr->setDimensions( size, size );
25 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
26
27 auto f = [=] __cuda_callable__ ( typename MatrixType::RowView& row ) mutable {
28 const int rowIdx = row.getRowIndex();
29 if( rowIdx == 0 )
30 {
31 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
32 row.setElement( 1, rowIdx+1, -1 ); // element above the diagonal
33 }
34 else if( rowIdx == size - 1 )
35 {
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 {
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 auto solver_ptr = TNL::Solvers::getLinearSolver< MatrixType >( "tfqmr" );
66 auto preconditioner_ptr = TNL::Solvers::getPreconditioner< MatrixType >( "diagonal");
67 preconditioner_ptr->update( matrix_ptr );
68 solver_ptr->setMatrix( matrix_ptr );
69 solver_ptr->setPreconditioner( preconditioner_ptr );
70 solver_ptr->setConvergenceResidue( 1.0e-6 );
71 solver_ptr->solve( b, x );
72 std::cout << "Vector x = " << x << std::endl;
73}
74
75int main( 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}

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

The result looks as follows:

Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Generated on Sat Apr 1 2023 11:40:57 for Template Numerical Library by doxygen 1.9.5