Template Numerical Library version\ main:88903d64
Loading...
Searching...
No Matches
Linear solvers

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 */
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 */
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}
#define __cuda_callable__
Definition Macros.h:49
Vector extends Array with algebraic operations.
Definition Vector.h:36
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:57
void setMatrix(const MatrixPointer &matrix)
Set the matrix of the linear system.
Definition LinearSolver.h:120
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition TFQMR.h:21
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
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 */
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 */
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}
Iterative solver of linear systems based on the Jacobi method.
Definition Jacobi.h:21
A RAII wrapper for launching the SolverMonitor's main loop in a separate thread.
Definition SolverMonitor.h:133

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: 18679 RES: 0.040976
Jacobi stage: ITER: 37191 RES: 0.0023872
Jacobi stage: ITER: 55539 RES: 0.00014254
Jacobi stage: ITER: 73219 RES: 9.4302e-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: 4189 RES: 0.38748
Jacobi stage: ITER: 8675 RES: 0.19079
Jacobi stage: ITER: 13177 RES: 0.095435
Jacobi stage: ITER: 16379 RES: 0.058375
Jacobi stage: ITER: 17177 RES: 0.051625
Jacobi stage: ITER: 17999 RES: 0.045515
Jacobi stage: ITER: 18823 RES: 0.040104
Jacobi stage: ITER: 19649 RES: 0.035315
Jacobi stage: ITER: 20469 RES: 0.031135
Jacobi stage: ITER: 21293 RES: 0.027434
Jacobi stage: ITER: 22115 RES: 0.024187
Jacobi stage: ITER: 22937 RES: 0.021312
Jacobi stage: ITER: 23759 RES: 0.01879
Jacobi stage: ITER: 24577 RES: 0.016566
Jacobi stage: ITER: 25401 RES: 0.014596
Jacobi stage: ITER: 26223 RES: 0.012869
Jacobi stage: ITER: 27047 RES: 0.011339
Jacobi stage: ITER: 27883 RES: 0.0099727
Jacobi stage: ITER: 28707 RES: 0.0087871
Jacobi stage: ITER: 29531 RES: 0.0077424
Jacobi stage: ITER: 30353 RES: 0.0068219
Jacobi stage: ITER: 31175 RES: 0.0060146
Jacobi stage: ITER: 31999 RES: 0.0052996
Jacobi stage: ITER: 32821 RES: 0.0046695
Jacobi stage: ITER: 33643 RES: 0.0041144
Jacobi stage: ITER: 34467 RES: 0.0036275
Jacobi stage: ITER: 35291 RES: 0.0031962
Jacobi stage: ITER: 36113 RES: 0.0028162
Jacobi stage: ITER: 36935 RES: 0.0024829
Jacobi stage: ITER: 37759 RES: 0.0021878
Jacobi stage: ITER: 38581 RES: 0.0019277
Jacobi stage: ITER: 39403 RES: 0.0016995
Jacobi stage: ITER: 40227 RES: 0.0014975
Jacobi stage: ITER: 41049 RES: 0.0013195
Jacobi stage: ITER: 41871 RES: 0.0011633
Jacobi stage: ITER: 42693 RES: 0.001025
Jacobi stage: ITER: 43515 RES: 0.0009037
Jacobi stage: ITER: 44339 RES: 0.00079627
Jacobi stage: ITER: 45159 RES: 0.00070203
Jacobi stage: ITER: 45987 RES: 0.00061819
Jacobi stage: ITER: 46811 RES: 0.0005447
Jacobi stage: ITER: 47635 RES: 0.00047994
Jacobi stage: ITER: 48457 RES: 0.00042288
Jacobi stage: ITER: 49279 RES: 0.00037284
Jacobi stage: ITER: 50099 RES: 0.00032871
Jacobi stage: ITER: 50923 RES: 0.00028963
Jacobi stage: ITER: 51743 RES: 0.00025536
Jacobi stage: ITER: 52565 RES: 0.000225
Jacobi stage: ITER: 53387 RES: 0.00019837
Jacobi stage: ITER: 54209 RES: 0.00017479
Jacobi stage: ITER: 55031 RES: 0.0001541
Jacobi stage: ITER: 55853 RES: 0.00013578
Jacobi stage: ITER: 56675 RES: 0.00011971
Jacobi stage: ITER: 57495 RES: 0.00010555
Jacobi stage: ITER: 58319 RES: 9.2999e-05
Jacobi stage: ITER: 59141 RES: 8.1942e-05
Jacobi stage: ITER: 59963 RES: 7.2245e-05
Jacobi stage: ITER: 60795 RES: 6.3578e-05
Jacobi stage: ITER: 61619 RES: 5.6019e-05
Jacobi stage: ITER: 62441 RES: 4.936e-05
Jacobi stage: ITER: 63263 RES: 4.3518e-05
Jacobi stage: ITER: 64087 RES: 3.8344e-05
Jacobi stage: ITER: 64909 RES: 3.3786e-05
Jacobi stage: ITER: 65731 RES: 2.9787e-05
Jacobi stage: ITER: 66553 RES: 2.6246e-05
Jacobi stage: ITER: 67375 RES: 2.314e-05
Jacobi stage: ITER: 68197 RES: 2.0389e-05
Jacobi stage: ITER: 69019 RES: 1.7976e-05
Jacobi stage: ITER: 69843 RES: 1.5839e-05
Jacobi stage: ITER: 70663 RES: 1.3965e-05
Jacobi stage: ITER: 71487 RES: 1.2304e-05
Jacobi stage: ITER: 72307 RES: 1.0848e-05
Jacobi stage: ITER: 73131 RES: 9.5586e-06
Jacobi stage: ITER: 73951 RES: 8.4274e-06
Jacobi stage: ITER: 74775 RES: 7.4255e-06
Jacobi stage: ITER: 75599 RES: 6.5427e-06
Jacobi stage: ITER: 76419 RES: 5.7684e-06
Jacobi stage: ITER: 77243 RES: 5.0826e-06
Jacobi stage: ITER: 78063 RES: 4.4811e-06
Jacobi stage: ITER: 78885 RES: 3.9484e-06
Jacobi stage: ITER: 79707 RES: 3.4811e-06
Jacobi stage: ITER: 80531 RES: 3.0673e-06
Jacobi stage: ITER: 81353 RES: 2.7026e-06
Jacobi stage: ITER: 82179 RES: 2.3813e-06
Jacobi stage: ITER: 83005 RES: 2.0969e-06
Jacobi stage: ITER: 83831 RES: 1.8476e-06
Jacobi stage: ITER: 84655 RES: 1.628e-06
Jacobi stage: ITER: 85477 RES: 1.4344e-06
Jacobi stage: ITER: 86299 RES: 1.2647e-06
Jacobi stage: ITER: 87123 RES: 1.1143e-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 */
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 */
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}
Class for real time, CPU time and CPU cycles measuring.
Definition Timer.h:25
void start()
Starts timer.
Definition Timer.hpp:42

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:4.6595e-
ELA: 0.10052 Jacobi stage: ITER: 18935 RES: 0.03942
ELA: 0.20064 Jacobi stage: ITER: 37721 RES: 0.0021999
ELA: 0.30076 Jacobi stage: ITER: 56597 RES: 0.00012112
ELA: 0.40088 Jacobi stage: ITER: 75471 RES: 6.6726e-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:3.3703e-
ELA: 0.10014 Jacobi stage: ITER: 771 RES: 0.78627
ELA: 0.20025 Jacobi stage: ITER: 1591 RES: 0.6378
ELA: 0.30037 Jacobi stage: ITER: 2417 RES: 0.5339
ELA: 0.40049 Jacobi stage: ITER: 3241 RES: 0.45675
ELA: 0.50061 Jacobi stage: ITER: 4063 RES: 0.39592
ELA: 0.60073 Jacobi stage: ITER: 4883 RES: 0.34569
ELA: 0.70087 Jacobi stage: ITER: 5705 RES: 0.3029
ELA: 0.801 Jacobi stage: ITER: 6527 RES: 0.26621
ELA: 0.90114 Jacobi stage: ITER: 7351 RES: 0.23414
ELA: 1.0013 Jacobi stage: ITER: 8199 RES: 0.20533
ELA: 1.1014 Jacobi stage: ITER: 9023 RES: 0.18082
ELA: 1.2015 Jacobi stage: ITER: 9845 RES: 0.15927
ELA: 1.3017 Jacobi stage: ITER: 10667 RES: 0.1404
ELA: 1.4018 Jacobi stage: ITER: 11487 RES: 0.12377
ELA: 1.502 Jacobi stage: ITER: 12311 RES: 0.10905
ELA: 1.6021 Jacobi stage: ITER: 13133 RES: 0.096082
ELA: 1.7022 Jacobi stage: ITER: 13955 RES: 0.08471
ELA: 1.8023 Jacobi stage: ITER: 14779 RES: 0.074639
ELA: 1.9025 Jacobi stage: ITER: 15603 RES: 0.065765
ELA: 2.0025 Jacobi stage: ITER: 16423 RES: 0.057982
ELA: 2.1027 Jacobi stage: ITER: 17247 RES: 0.051088
ELA: 2.2028 Jacobi stage: ITER: 18071 RES: 0.045015
ELA: 2.3029 Jacobi stage: ITER: 18893 RES: 0.039663
ELA: 2.4031 Jacobi stage: ITER: 19715 RES: 0.034969
ELA: 2.5032 Jacobi stage: ITER: 20539 RES: 0.030812
ELA: 2.6033 Jacobi stage: ITER: 21361 RES: 0.027149
ELA: 2.7035 Jacobi stage: ITER: 22183 RES: 0.023936
ELA: 2.8036 Jacobi stage: ITER: 23007 RES: 0.02109
ELA: 2.9037 Jacobi stage: ITER: 23829 RES: 0.018583
ELA: 3.0039 Jacobi stage: ITER: 24651 RES: 0.016384
ELA: 3.104 Jacobi stage: ITER: 25475 RES: 0.014436
ELA: 3.2042 Jacobi stage: ITER: 26297 RES: 0.01272
ELA: 3.3043 Jacobi stage: ITER: 27119 RES: 0.011214
ELA: 3.4044 Jacobi stage: ITER: 27941 RES: 0.0098812
ELA: 3.5046 Jacobi stage: ITER: 28763 RES: 0.0087118
ELA: 3.6092 Jacobi stage: ITER: 29621 RES: 0.0076338
ELA: 3.7093 Jacobi stage: ITER: 30443 RES: 0.0067304
ELA: 3.8095 Jacobi stage: ITER: 31265 RES: 0.0059302
ELA: 3.9096 Jacobi stage: ITER: 32087 RES: 0.0052284
ELA: 4.0097 Jacobi stage: ITER: 32909 RES: 0.0046068
ELA: 4.1098 Jacobi stage: ITER: 33731 RES: 0.0040616
ELA: 4.2099 Jacobi stage: ITER: 34553 RES: 0.0035788
ELA: 4.3101 Jacobi stage: ITER: 35375 RES: 0.0031552
ELA: 4.4102 Jacobi stage: ITER: 36195 RES: 0.0027818
ELA: 4.5104 Jacobi stage: ITER: 37019 RES: 0.0024511
ELA: 4.6105 Jacobi stage: ITER: 37841 RES: 0.0021597
ELA: 4.7106 Jacobi stage: ITER: 38663 RES: 0.0019041
ELA: 4.8108 Jacobi stage: ITER: 39487 RES: 0.0016778
ELA: 4.9109 Jacobi stage: ITER: 40307 RES: 0.0014792
ELA: 5.0111 Jacobi stage: ITER: 41139 RES: 0.0013017
ELA: 5.1112 Jacobi stage: ITER: 41963 RES: 0.001147
ELA: 5.2113 Jacobi stage: ITER: 42783 RES: 0.0010112
ELA: 5.3115 Jacobi stage: ITER: 43607 RES: 0.00089102
ELA: 5.4116 Jacobi stage: ITER: 44427 RES: 0.00078558
ELA: 5.5124 Jacobi stage: ITER: 45255 RES: 0.00069176
ELA: 5.6125 Jacobi stage: ITER: 46079 RES: 0.00060952
ELA: 5.7127 Jacobi stage: ITER: 46901 RES: 0.00053705
ELA: 5.8128 Jacobi stage: ITER: 47723 RES: 0.0004735
ELA: 5.9129 Jacobi stage: ITER: 48545 RES: 0.00041721
ELA: 6.0131 Jacobi stage: ITER: 49369 RES: 0.00036761
ELA: 6.1132 Jacobi stage: ITER: 50191 RES: 0.0003241
ELA: 6.2133 Jacobi stage: ITER: 51013 RES: 0.00028557
ELA: 6.3135 Jacobi stage: ITER: 51835 RES: 0.00025178
ELA: 6.4136 Jacobi stage: ITER: 52657 RES: 0.00022184
ELA: 6.5137 Jacobi stage: ITER: 53479 RES: 0.00019559
ELA: 6.6139 Jacobi stage: ITER: 54301 RES: 0.00017234
ELA: 6.714 Jacobi stage: ITER: 55123 RES: 0.00015194
ELA: 6.8141 Jacobi stage: ITER: 55945 RES: 0.00013388
ELA: 6.9142 Jacobi stage: ITER: 56767 RES: 0.00011803
ELA: 7.0144 Jacobi stage: ITER: 57589 RES: 0.000104
ELA: 7.1145 Jacobi stage: ITER: 58411 RES: 9.1694e-05
ELA: 7.2146 Jacobi stage: ITER: 59233 RES: 8.0793e-05
ELA: 7.3147 Jacobi stage: ITER: 60055 RES: 7.1231e-05
ELA: 7.4148 Jacobi stage: ITER: 60877 RES: 6.2763e-05
ELA: 7.515 Jacobi stage: ITER: 61699 RES: 5.5335e-05
ELA: 7.6151 Jacobi stage: ITER: 62529 RES: 4.8697e-05
ELA: 7.7152 Jacobi stage: ITER: 63351 RES: 4.2934e-05
ELA: 7.8154 Jacobi stage: ITER: 64173 RES: 3.783e-05
ELA: 7.9155 Jacobi stage: ITER: 64995 RES: 3.3353e-05
ELA: 8.0157 Jacobi stage: ITER: 65815 RES: 2.9406e-05
ELA: 8.1158 Jacobi stage: ITER: 66639 RES: 2.591e-05
ELA: 8.2159 Jacobi stage: ITER: 67461 RES: 2.2829e-05
ELA: 8.316 Jacobi stage: ITER: 68283 RES: 2.0128e-05
ELA: 8.4161 Jacobi stage: ITER: 69105 RES: 1.7735e-05
ELA: 8.5163 Jacobi stage: ITER: 69927 RES: 1.5636e-05
ELA: 8.6164 Jacobi stage: ITER: 70751 RES: 1.3777e-05
ELA: 8.7165 Jacobi stage: ITER: 71573 RES: 1.2139e-05
ELA: 8.8167 Jacobi stage: ITER: 75975 RES: 6.1755e-06
ELA: 8.9168 Jacobi stage: ITER: 80855 RES: 2.9183e-06
ELA: 9.0175 Jacobi stage: ITER: 85737 RES: 1.3783e-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 */
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 */
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}
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 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
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 */
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}

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 ]