Template Numerical Library version\ main:481315e2
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 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:39
TNL::Matrices::SparseMatrix
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:60
TNL::Solvers::Linear::LinearSolver::setMatrix
void setMatrix(const MatrixPointer &matrix)
Set the matrix of the linear system.
Definition LinearSolver.h:123
TNL::Solvers::Linear::TFQMR
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition TFQMR.h:24
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:24
TNL::Solvers::SolverMonitorThread
A RAII wrapper for launching the SolverMonitor's main loop in a separate thread.
Definition SolverMonitor.h:136

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: 17325 RES: 0.050464
Jacobi stage: ITER: 33927 RES: 0.0039412
Jacobi stage: ITER: 50747 RES: 0.00029757
Jacobi stage: ITER: 67343 RES: 2.3254e-05
Jacobi stage: ITER: 82695 RES: 2.1985e-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: 2371 RES: 0.53907
Jacobi stage: ITER: 4731 RES: 0.35436
Jacobi stage: ITER: 7359 RES: 0.23385
Jacobi stage: ITER: 9539 RES: 0.16701
Jacobi stage: ITER: 11687 RES: 0.12002
Jacobi stage: ITER: 14137 RES: 0.082349
Jacobi stage: ITER: 16335 RES: 0.058771
Jacobi stage: ITER: 18775 RES: 0.040401
Jacobi stage: ITER: 20983 RES: 0.028781
Jacobi stage: ITER: 23531 RES: 0.019459
Jacobi stage: ITER: 25635 RES: 0.014085
Jacobi stage: ITER: 28025 RES: 0.0097545
Jacobi stage: ITER: 30095 RES: 0.0070999
Jacobi stage: ITER: 31491 RES: 0.0057296
Jacobi stage: ITER: 33133 RES: 0.004451
Jacobi stage: ITER: 34267 RES: 0.0037406
Jacobi stage: ITER: 35943 RES: 0.0028916
Jacobi stage: ITER: 37567 RES: 0.0022532
Jacobi stage: ITER: 38501 RES: 0.0019515
Jacobi stage: ITER: 40159 RES: 0.0015123
Jacobi stage: ITER: 41759 RES: 0.0011835
Jacobi stage: ITER: 43423 RES: 0.00091657
Jacobi stage: ITER: 44607 RES: 0.00076415
Jacobi stage: ITER: 46161 RES: 0.0006017
Jacobi stage: ITER: 47335 RES: 0.00050257
Jacobi stage: ITER: 48889 RES: 0.00039573
Jacobi stage: ITER: 50183 RES: 0.0003245
Jacobi stage: ITER: 51607 RES: 0.00026075
Jacobi stage: ITER: 53147 RES: 0.00020582
Jacobi stage: ITER: 54729 RES: 0.00016137
Jacobi stage: ITER: 56889 RES: 0.00011581
Jacobi stage: ITER: 58871 RES: 8.5438e-05
Jacobi stage: ITER: 59385 RES: 7.8928e-05
Jacobi stage: ITER: 59883 RES: 7.3138e-05
Jacobi stage: ITER: 61715 RES: 5.5199e-05
Jacobi stage: ITER: 63519 RES: 4.184e-05
Jacobi stage: ITER: 64691 RES: 3.4947e-05
Jacobi stage: ITER: 65189 RES: 3.2364e-05
Jacobi stage: ITER: 65675 RES: 3.0045e-05
Jacobi stage: ITER: 66395 RES: 2.6899e-05
Jacobi stage: ITER: 66893 RES: 2.4911e-05
Jacobi stage: ITER: 67381 RES: 2.3112e-05
Jacobi stage: ITER: 67891 RES: 2.1377e-05
Jacobi stage: ITER: 68381 RES: 1.9821e-05
Jacobi stage: ITER: 68861 RES: 1.8412e-05
Jacobi stage: ITER: 69349 RES: 1.7082e-05
Jacobi stage: ITER: 69963 RES: 1.555e-05
Jacobi stage: ITER: 70483 RES: 1.4356e-05
Jacobi stage: ITER: 71197 RES: 1.2861e-05
Jacobi stage: ITER: 72181 RES: 1.1057e-05
Jacobi stage: ITER: 73339 RES: 9.258e-06
Jacobi stage: ITER: 74947 RES: 7.2319e-06
Jacobi stage: ITER: 76805 RES: 5.4347e-06
Jacobi stage: ITER: 78785 RES: 4.0095e-06
Jacobi stage: ITER: 79727 RES: 3.4704e-06
Jacobi stage: ITER: 80385 RES: 3.1359e-06
Jacobi stage: ITER: 81171 RES: 2.7801e-06
Jacobi stage: ITER: 81919 RES: 2.4783e-06
Jacobi stage: ITER: 83487 RES: 1.9479e-06
Jacobi stage: ITER: 85655 RES: 1.3962e-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:48

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:0.001140
ELA: 0.10446 Jacobi stage: ITER: 16071 RES: 0.061203
ELA: 0.2078 Jacobi stage: ITER: 33635 RES: 0.004122
ELA: 0.31113 Jacobi stage: ITER: 47799 RES: 0.000468
ELA: 0.41125 Jacobi stage: ITER: 57053 RES: 0.00011293
ELA: 0.51137 Jacobi stage: ITER: 65047 RES: 3.3087e-05
ELA: 0.6115 Jacobi stage: ITER: 72949 RES: 9.8265e-06
ELA: 0.71162 Jacobi stage: ITER: 82167 RES: 2.3857e-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: 0.0127
ELA: 0.11281 Jacobi stage: ITER: 333 RES: 0.89565
ELA: 0.21292 Jacobi stage: ITER: 1323 RES: 0.67998
ELA: 0.31304 Jacobi stage: ITER: 2343 RES: 0.54213
ELA: 0.41327 Jacobi stage: ITER: 2777 RES: 0.49772
ELA: 0.51341 Jacobi stage: ITER: 3115 RES: 0.4675
ELA: 0.61353 Jacobi stage: ITER: 3465 RES: 0.43884
ELA: 0.71365 Jacobi stage: ITER: 3915 RES: 0.406
ELA: 0.81377 Jacobi stage: ITER: 4267 RES: 0.38258
ELA: 0.9139 Jacobi stage: ITER: 4627 RES: 0.36044
ELA: 1.014 Jacobi stage: ITER: 4985 RES: 0.33991
ELA: 1.1141 Jacobi stage: ITER: 5385 RES: 0.31873
ELA: 1.2143 Jacobi stage: ITER: 5761 RES: 0.30022
ELA: 1.3144 Jacobi stage: ITER: 6179 RES: 0.28114
ELA: 1.4145 Jacobi stage: ITER: 6555 RES: 0.26505
ELA: 1.5146 Jacobi stage: ITER: 6937 RES: 0.24963
ELA: 1.6148 Jacobi stage: ITER: 7321 RES: 0.23516
ELA: 1.7149 Jacobi stage: ITER: 7693 RES: 0.22198
ELA: 1.815 Jacobi stage: ITER: 8059 RES: 0.20983
ELA: 1.926 Jacobi stage: ITER: 8415 RES: 0.1986
ELA: 2.0261 Jacobi stage: ITER: 8733 RES: 0.18903
ELA: 2.1294 Jacobi stage: ITER: 9247 RES: 0.17469
ELA: 2.2295 Jacobi stage: ITER: 9747 RES: 0.16174
ELA: 2.3296 Jacobi stage: ITER: 10079 RES: 0.15369
ELA: 2.4298 Jacobi stage: ITER: 10415 RES: 0.14595
ELA: 2.5299 Jacobi stage: ITER: 10815 RES: 0.13724
ELA: 2.63 Jacobi stage: ITER: 11315 RES: 0.12709
ELA: 2.7301 Jacobi stage: ITER: 11799 RES: 0.11798
ELA: 2.8302 Jacobi stage: ITER: 12289 RES: 0.10939
ELA: 2.9303 Jacobi stage: ITER: 12749 RES: 0.10192
ELA: 3.0305 Jacobi stage: ITER: 13239 RES: 0.09456
ELA: 3.1306 Jacobi stage: ITER: 13737 RES: 0.087568
ELA: 3.2307 Jacobi stage: ITER: 14225 RES: 0.081244
ELA: 3.3308 Jacobi stage: ITER: 14715 RES: 0.075376
ELA: 3.4309 Jacobi stage: ITER: 15203 RES: 0.069932
ELA: 3.531 Jacobi stage: ITER: 15695 RES: 0.064842
ELA: 3.6312 Jacobi stage: ITER: 16845 RES: 0.054326
ELA: 3.7313 Jacobi stage: ITER: 18135 RES: 0.044574
ELA: 3.8314 Jacobi stage: ITER: 19139 RES: 0.038204
ELA: 3.9315 Jacobi stage: ITER: 19483 RES: 0.036238
ELA: 4.0316 Jacobi stage: ITER: 20209 RES: 0.032404
ELA: 4.1317 Jacobi stage: ITER: 21055 RES: 0.028464
ELA: 4.2319 Jacobi stage: ITER: 21715 RES: 0.02572
ELA: 4.332 Jacobi stage: ITER: 22295 RES: 0.023528
ELA: 4.4321 Jacobi stage: ITER: 24193 RES: 0.017572
ELA: 4.5322 Jacobi stage: ITER: 25185 RES: 0.015089
ELA: 4.6323 Jacobi stage: ITER: 26709 RES: 0.01194
ELA: 4.7325 Jacobi stage: ITER: 29591 RES: 0.0076714
ELA: 4.8326 Jacobi stage: ITER: 32161 RES: 0.0051677
ELA: 4.9327 Jacobi stage: ITER: 33103 RES: 0.0044729
ELA: 5.0328 Jacobi stage: ITER: 36195 RES: 0.0027818
ELA: 5.1329 Jacobi stage: ITER: 40787 RES: 0.0013741
ELA: 5.233 Jacobi stage: ITER: 44827 RES: 0.00073876
ELA: 5.3331 Jacobi stage: ITER: 47359 RES: 0.00050073
ELA: 5.4333 Jacobi stage: ITER: 49651 RES: 0.00035191
ELA: 5.5335 Jacobi stage: ITER: 53935 RES: 0.00018236
ELA: 5.6336 Jacobi stage: ITER: 58495 RES: 9.0518e-05
ELA: 5.7337 Jacobi stage: ITER: 60927 RES: 6.2302e-05
ELA: 5.8338 Jacobi stage: ITER: 65343 RES: 3.1617e-05
ELA: 5.9339 Jacobi stage: ITER: 70095 RES: 1.5238e-05
ELA: 6.034 Jacobi stage: ITER: 74547 RES: 7.6901e-06
ELA: 6.1342 Jacobi stage: ITER: 77771 RES: 4.6867e-06
ELA: 6.2343 Jacobi stage: ITER: 82155 RES: 2.3901e-06
ELA: 6.3344 Jacobi stage: ITER: 85531 RES: 1.423e-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:24

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 Thu Sep 7 2023 14:21:33 for Template Numerical Library by doxygen 1.9.8