Linear solvers
Table of Contents
Introduction
Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the following iterative methods:
- Stationary methods
- Krylov subspace methods
- Conjugate gradient method, CG (TNL::Solvers::Linear::CG)
- Biconjugate gradient stabilized method, BICGStab (TNL::Solvers::Linear::BICGStab)
- Biconjugate gradient stabilized method, BICGStab(l) (TNL::Solvers::Linear::BICGStabL)
- Transpose-free quasi-minimal residual method, TFQMR (TNL::Solvers::Linear::TFQMR)
- Generalized minimal residual method, GMRES (TNL::Solvers::Linear::GMRES) with various methods of orthogonalization:
- Classical Gramm-Schmidt, CGS
- Classical Gramm-Schmidt with reorthogonalization, CGSR
- Modified Gramm-Schmidt, MGS
- Modified Gramm-Schmidt with reorthogonalization, MGSR
- 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:
- Diagonal or Jacobi (TNL::Solvers::Linear::Preconditioners::Diagonal)
- ILU (Incomplete LU) - CPU only currently
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 );
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
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 );
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;
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 );
70}
71
72int
73main( int argc, char* argv[] )
74{
76 iterativeLinearSolverExample< TNL::Devices::Sequential >();
77
78#ifdef __CUDACC__
80 iterativeLinearSolverExample< TNL::Devices::Cuda >();
81#endif
82}
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:55
virtual 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
In this example we solve a linear system \( A \vec x = \vec b \) where
\[A = \left( \begin{array}{cccc} 2.5 & -1 & & & \\ -1 & 2.5 & -1 & & \\ & -1 & 2.5 & -1 & \\ & & -1 & 2.5 & -1 \\ & & & -1 & 2.5 \\ \end{array} \right) \]
The right-hand side vector \(\vec b \) is set to \(( 1.5, 0.5, 0.5, 0.5, 1.5 )^T \) so that the exact solution is \( \vec x = ( 1, 1, 1, 1, 1 )^T \). The elements of the matrix \( A \) are set using the method TNL::Matrices::SparseMatrix::forAllRows. In this example, we use the sparse matrix but any other matrix type can be used as well (see the namespace TNL::Matrices). Next we set the solution vector \( \vec x = ( 1, 1, 1, 1, 1 )^T \) and multiply it with matrix \( A \) to get the right-hand side vector \( \vec b \). Finally, we reset the vector \( \vec x \) to zero.
To solve the linear system in the example, we use the TFQMR solver. Other solvers can be used as well (see the namespace TNL::Solvers::Linear). The solver needs only one template parameter which is the matrix type. Next we create an instance of the solver and set the matrix of the linear system. Note that the matrix is passed to the solver as a std::shared_ptr. Then we set the stopping criterion for the iterative method in terms of the relative residue, i.e. \( \lVert \vec b - A \vec x \rVert / \lVert b \rVert \). The solver is executed by calling the TNL::Solvers::Linear::LinearSolver::solve method which accepts the right-hand side vector \( \vec b \) and the solution vector \( \vec x \).
The result looks as follows:
Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Setup with a solver monitor
Solution of large linear systems may take a lot of time. In such situations, it is useful to be able to monitor the convergence of the solver or the solver status in general. For this purpose, TNL provides a solver monitor which can show various metrics in real time, such as current number of iterations, current residue of the approximate solution, etc. The solver monitor in TNL runs in a separate thread and it refreshes the status of the solver with a configurable refresh rate (once per 500 ms by default). The use of the solver monitor is demonstrated in the following example.
1#include <iostream>
2#include <memory>
3#include <TNL/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Sequential.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/Linear/Jacobi.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 );
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
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 );
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;
60
61 /***
62 * Setup solver of the linear system.
63 */
65 LinearSolver solver;
66 solver.setMatrix( matrix_ptr );
67 solver.setOmega( 0.0005 );
68
69 /***
70 * Setup monitor of the iterative solver.
71 */
73 IterativeSolverMonitorType monitor;
74 TNL::Solvers::SolverMonitorThread monitorThread( monitor );
75 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
76 monitor.setVerbose( 1 );
77 monitor.setStage( "Jacobi stage:" );
78 solver.setSolverMonitor( monitor );
79 solver.setConvergenceResidue( 1.0e-6 );
80 solver.solve( b, x );
81 monitor.stopMainLoop();
83}
84
85int
86main( int argc, char* argv[] )
87{
89 iterativeLinearSolverExample< TNL::Devices::Sequential >();
90
91#ifdef __CUDACC__
93 iterativeLinearSolverExample< TNL::Devices::Cuda >();
94#endif
95}
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:142
First, we set up the same linear system as in the previous example, we create an instance of the Jacobi solver and we pass the matrix of the linear system to the solver. Then, we set the relaxation parameter \( \omega \) of the Jacobi solver to 0.0005. The reason is to artificially slow down the convergence, because we want to see some iterations in this example. Next, we create an instance of the solver monitor and a special thread for the monitor (an instance of the TNL::Solvers::SolverMonitorThread class). We use the following methods to configure the solver monitor:
- TNL::Solvers::SolverMonitor::setRefreshRate sets the refresh rate of the monitor to 10 milliseconds.
- TNL::Solvers::IterativeSolverMonitor::setVerbose sets verbosity of the monitor to 1.
- TNL::Solvers::IterativeSolverMonitor::setStage sets a name of the solver stage. The monitor stages serve for distinguishing between different stages of more complex solvers (for example when the linear system solver is embedded into a time-dependent PDE solver).
Next, we call TNL::Solvers::IterativeSolver::setSolverMonitor to connect the solver with the monitor and we set the convergence criterion based on the relative residue. Finally, we start the solver using the TNL::Solvers::Linear::Jacobi::solve method and when the solver finishes, we stop the monitor using TNL::Solvers::SolverMonitor::stopMainLoop.
The result looks as follows:
Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Jacobi stage: ITER: 79741 RES: 3.4619e-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: 7339 RES: 0.23458
Jacobi stage: ITER: 15003 RES: 0.072114
Jacobi stage: ITER: 21987 RES: 0.024667
Jacobi stage: ITER: 22579 RES: 0.022523
Jacobi stage: ITER: 22987 RES: 0.021155
Jacobi stage: ITER: 23407 RES: 0.019833
Jacobi stage: ITER: 23819 RES: 0.018617
Jacobi stage: ITER: 24211 RES: 0.017529
Jacobi stage: ITER: 24609 RES: 0.016485
Jacobi stage: ITER: 25027 RES: 0.015464
Jacobi stage: ITER: 25437 RES: 0.014516
Jacobi stage: ITER: 25847 RES: 0.013634
Jacobi stage: ITER: 26255 RES: 0.012806
Jacobi stage: ITER: 26663 RES: 0.012028
Jacobi stage: ITER: 27071 RES: 0.011297
Jacobi stage: ITER: 27491 RES: 0.010592
Jacobi stage: ITER: 27903 RES: 0.0099421
Jacobi stage: ITER: 28295 RES: 0.0093611
Jacobi stage: ITER: 28675 RES: 0.0088304
Jacobi stage: ITER: 29081 RES: 0.0082939
Jacobi stage: ITER: 29499 RES: 0.0077806
Jacobi stage: ITER: 29909 RES: 0.0073034
Jacobi stage: ITER: 30109 RES: 0.0070825
Jacobi stage: ITER: 30289 RES: 0.0068893
Jacobi stage: ITER: 30469 RES: 0.0067015
Jacobi stage: ITER: 30649 RES: 0.0065187
Jacobi stage: ITER: 30827 RES: 0.0063449
Jacobi stage: ITER: 31007 RES: 0.0061718
Jacobi stage: ITER: 31183 RES: 0.0060072
Jacobi stage: ITER: 31365 RES: 0.0058398
Jacobi stage: ITER: 31539 RES: 0.0056876
Jacobi stage: ITER: 31719 RES: 0.0055325
Jacobi stage: ITER: 31899 RES: 0.0053816
Jacobi stage: ITER: 32081 RES: 0.0052316
Jacobi stage: ITER: 32351 RES: 0.0050206
Jacobi stage: ITER: 32753 RES: 0.0047185
Jacobi stage: ITER: 33165 RES: 0.0044292
Jacobi stage: ITER: 33569 RES: 0.0041627
Jacobi stage: ITER: 33989 RES: 0.0039026
Jacobi stage: ITER: 34401 RES: 0.0036633
Jacobi stage: ITER: 34795 RES: 0.0034492
Jacobi stage: ITER: 35175 RES: 0.0032537
Jacobi stage: ITER: 35561 RES: 0.0030654
Jacobi stage: ITER: 35947 RES: 0.0028898
Jacobi stage: ITER: 36331 RES: 0.0027243
Jacobi stage: ITER: 36719 RES: 0.0025667
Jacobi stage: ITER: 37103 RES: 0.0024197
Jacobi stage: ITER: 37489 RES: 0.0022797
Jacobi stage: ITER: 37875 RES: 0.0021491
Jacobi stage: ITER: 38259 RES: 0.002026
Jacobi stage: ITER: 38645 RES: 0.0019088
Jacobi stage: ITER: 39031 RES: 0.0017995
Jacobi stage: ITER: 39415 RES: 0.0016964
Jacobi stage: ITER: 39801 RES: 0.0015983
Jacobi stage: ITER: 40187 RES: 0.0015067
Jacobi stage: ITER: 40571 RES: 0.0014204
Jacobi stage: ITER: 40959 RES: 0.0013382
Jacobi stage: ITER: 41343 RES: 0.0012616
Jacobi stage: ITER: 41729 RES: 0.0011886
Jacobi stage: ITER: 42115 RES: 0.0011205
Jacobi stage: ITER: 42499 RES: 0.0010563
Jacobi stage: ITER: 42885 RES: 0.00099522
Jacobi stage: ITER: 43271 RES: 0.00093822
Jacobi stage: ITER: 43655 RES: 0.00088448
Jacobi stage: ITER: 44041 RES: 0.00083331
Jacobi stage: ITER: 44427 RES: 0.00078558
Jacobi stage: ITER: 44811 RES: 0.00074058
Jacobi stage: ITER: 45197 RES: 0.00069773
Jacobi stage: ITER: 45583 RES: 0.00065777
Jacobi stage: ITER: 45967 RES: 0.00062009
Jacobi stage: ITER: 46353 RES: 0.00058422
Jacobi stage: ITER: 46739 RES: 0.00055076
Jacobi stage: ITER: 47123 RES: 0.00051921
Jacobi stage: ITER: 47509 RES: 0.00048917
Jacobi stage: ITER: 47895 RES: 0.00046115
Jacobi stage: ITER: 48279 RES: 0.00043474
Jacobi stage: ITER: 48665 RES: 0.00040959
Jacobi stage: ITER: 49051 RES: 0.00038613
Jacobi stage: ITER: 49435 RES: 0.00036401
Jacobi stage: ITER: 49821 RES: 0.00034295
Jacobi stage: ITER: 50207 RES: 0.00032331
Jacobi stage: ITER: 50591 RES: 0.00030479
Jacobi stage: ITER: 50977 RES: 0.00028715
Jacobi stage: ITER: 51363 RES: 0.00027071
Jacobi stage: ITER: 51747 RES: 0.0002552
Jacobi stage: ITER: 52133 RES: 0.00024044
Jacobi stage: ITER: 52519 RES: 0.00022666
Jacobi stage: ITER: 52903 RES: 0.00021368
Jacobi stage: ITER: 53289 RES: 0.00020132
Jacobi stage: ITER: 53675 RES: 0.00018979
Jacobi stage: ITER: 54059 RES: 0.00017892
Jacobi stage: ITER: 54445 RES: 0.00016857
Jacobi stage: ITER: 54831 RES: 0.00015891
Jacobi stage: ITER: 55215 RES: 0.00014981
Jacobi stage: ITER: 55601 RES: 0.00014114
Jacobi stage: ITER: 55987 RES: 0.00013306
Jacobi stage: ITER: 56373 RES: 0.00012536
Jacobi stage: ITER: 56759 RES: 0.00011818
Jacobi stage: ITER: 57143 RES: 0.00011141
Jacobi stage: ITER: 57529 RES: 0.00010496
Jacobi stage: ITER: 57915 RES: 9.8952e-05
Jacobi stage: ITER: 58299 RES: 9.3285e-05
Jacobi stage: ITER: 58685 RES: 8.7888e-05
Jacobi stage: ITER: 59071 RES: 8.2854e-05
Jacobi stage: ITER: 59455 RES: 7.8108e-05
Jacobi stage: ITER: 59843 RES: 7.3589e-05
Jacobi stage: ITER: 60227 RES: 6.9374e-05
Jacobi stage: ITER: 60613 RES: 6.536e-05
Jacobi stage: ITER: 60999 RES: 6.1617e-05
Jacobi stage: ITER: 61383 RES: 5.8087e-05
Jacobi stage: ITER: 61769 RES: 5.4727e-05
Jacobi stage: ITER: 62155 RES: 5.1592e-05
Jacobi stage: ITER: 62539 RES: 4.8637e-05
Jacobi stage: ITER: 62927 RES: 4.5823e-05
Jacobi stage: ITER: 63311 RES: 4.3198e-05
Jacobi stage: ITER: 63697 RES: 4.0699e-05
Jacobi stage: ITER: 64083 RES: 3.8368e-05
Jacobi stage: ITER: 64467 RES: 3.617e-05
Jacobi stage: ITER: 64855 RES: 3.4078e-05
Jacobi stage: ITER: 65239 RES: 3.2126e-05
Jacobi stage: ITER: 65625 RES: 3.0267e-05
Jacobi stage: ITER: 66011 RES: 2.8534e-05
Jacobi stage: ITER: 66395 RES: 2.6899e-05
Jacobi stage: ITER: 66781 RES: 2.5343e-05
Jacobi stage: ITER: 67167 RES: 2.3891e-05
Jacobi stage: ITER: 67551 RES: 2.2523e-05
Jacobi stage: ITER: 67937 RES: 2.122e-05
Jacobi stage: ITER: 68323 RES: 2.0004e-05
Jacobi stage: ITER: 68709 RES: 1.8847e-05
Jacobi stage: ITER: 69095 RES: 1.7768e-05
Jacobi stage: ITER: 69479 RES: 1.675e-05
Jacobi stage: ITER: 69865 RES: 1.5781e-05
Jacobi stage: ITER: 70251 RES: 1.4877e-05
Jacobi stage: ITER: 70635 RES: 1.4025e-05
Jacobi stage: ITER: 71021 RES: 1.3213e-05
Jacobi stage: ITER: 71407 RES: 1.2457e-05
Jacobi stage: ITER: 71791 RES: 1.1743e-05
Jacobi stage: ITER: 72179 RES: 1.1064e-05
Jacobi stage: ITER: 72563 RES: 1.043e-05
Jacobi stage: ITER: 72949 RES: 9.8265e-06
Jacobi stage: ITER: 73335 RES: 9.2637e-06
Jacobi stage: ITER: 73719 RES: 8.7331e-06
Jacobi stage: ITER: 74105 RES: 8.2278e-06
Jacobi stage: ITER: 74491 RES: 7.7566e-06
Jacobi stage: ITER: 74877 RES: 7.3078e-06
Jacobi stage: ITER: 75263 RES: 6.8892e-06
Jacobi stage: ITER: 75647 RES: 6.4946e-06
Jacobi stage: ITER: 76033 RES: 6.1189e-06
Jacobi stage: ITER: 76419 RES: 5.7684e-06
Jacobi stage: ITER: 76803 RES: 5.438e-06
Jacobi stage: ITER: 77189 RES: 5.1234e-06
Jacobi stage: ITER: 77575 RES: 4.8299e-06
Jacobi stage: ITER: 77961 RES: 4.5505e-06
Jacobi stage: ITER: 78347 RES: 4.2898e-06
Jacobi stage: ITER: 78731 RES: 4.0441e-06
Jacobi stage: ITER: 79119 RES: 3.8102e-06
Jacobi stage: ITER: 79503 RES: 3.5919e-06
Jacobi stage: ITER: 79889 RES: 3.3841e-06
Jacobi stage: ITER: 80275 RES: 3.1903e-06
Jacobi stage: ITER: 80659 RES: 3.0075e-06
Jacobi stage: ITER: 81045 RES: 2.8335e-06
Jacobi stage: ITER: 81431 RES: 2.6712e-06
Jacobi stage: ITER: 81817 RES: 2.5167e-06
Jacobi stage: ITER: 82721 RES: 2.1904e-06
Jacobi stage: ITER: 83107 RES: 2.065e-06
Jacobi stage: ITER: 83491 RES: 1.9467e-06
Jacobi stage: ITER: 83877 RES: 1.834e-06
Jacobi stage: ITER: 84263 RES: 1.729e-06
Jacobi stage: ITER: 84647 RES: 1.63e-06
Jacobi stage: ITER: 85033 RES: 1.5357e-06
Jacobi stage: ITER: 85419 RES: 1.4477e-06
Jacobi stage: ITER: 85803 RES: 1.3648e-06
Jacobi stage: ITER: 86189 RES: 1.2858e-06
Jacobi stage: ITER: 86575 RES: 1.2122e-06
Jacobi stage: ITER: 86959 RES: 1.1427e-06
Jacobi stage: ITER: 87345 RES: 1.0766e-06
Jacobi stage: ITER: 87731 RES: 1.015e-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 <TNL/Timer.h>
4#include <TNL/Matrices/SparseMatrix.h>
5#include <TNL/Devices/Sequential.h>
6#include <TNL/Devices/Cuda.h>
7#include <TNL/Solvers/Linear/Jacobi.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 );
26 matrix_ptr->setDimensions( size, size );
27 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
28
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 );
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;
61
62 /***
63 * Setup solver of the linear system.
64 */
66 LinearSolver solver;
67 solver.setMatrix( matrix_ptr );
68 solver.setOmega( 0.0005 );
69
70 /***
71 * Setup monitor of the iterative solver.
72 */
74 IterativeSolverMonitorType monitor;
75 TNL::Solvers::SolverMonitorThread mmonitorThread( monitor );
76 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
77 monitor.setVerbose( 1 );
78 monitor.setStage( "Jacobi stage:" );
79 TNL::Timer timer;
80 monitor.setTimer( timer );
81 timer.start();
82 solver.setSolverMonitor( monitor );
83 solver.setConvergenceResidue( 1.0e-6 );
84 solver.solve( b, x );
85 monitor.stopMainLoop();
87}
88
89int
90main( int argc, char* argv[] )
91{
93 iterativeLinearSolverExample< TNL::Devices::Sequential >();
94
95#ifdef __CUDACC__
97 iterativeLinearSolverExample< TNL::Devices::Cuda >();
98#endif
99}
The only changes are around the lines where we create an instance of TNL::Timer, connect it with the monitor using TNL::Solvers::SolverMonitor::setTimer and start the timer with TNL::Timer::start.
The result looks as follows:
Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
ELA:1.8375e-
ELA: 0.1001 Jacobi stage: ITER: 74601 RES: 7.6243e-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:1.059e-0
ELA: 0.1001 Jacobi stage: ITER: 407 RES: 0.87565
ELA: 0.20017 Jacobi stage: ITER: 809 RES: 0.77751
ELA: 0.30026 Jacobi stage: ITER: 1219 RES: 0.69783
ELA: 0.40037 Jacobi stage: ITER: 1623 RES: 0.6331
ELA: 0.50046 Jacobi stage: ITER: 2011 RES: 0.58087
ELA: 0.60054 Jacobi stage: ITER: 2405 RES: 0.53519
ELA: 0.70064 Jacobi stage: ITER: 2819 RES: 0.49396
ELA: 0.80072 Jacobi stage: ITER: 3229 RES: 0.45774
ELA: 0.9008 Jacobi stage: ITER: 3639 RES: 0.42579
ELA: 1.0009 Jacobi stage: ITER: 4043 RES: 0.39727
ELA: 1.1009 Jacobi stage: ITER: 4449 RES: 0.37105
ELA: 1.201 Jacobi stage: ITER: 4851 RES: 0.34749
ELA: 1.3011 Jacobi stage: ITER: 5269 RES: 0.32471
ELA: 1.4013 Jacobi stage: ITER: 5679 RES: 0.30425
ELA: 1.5013 Jacobi stage: ITER: 6067 RES: 0.28614
ELA: 1.6014 Jacobi stage: ITER: 6433 RES: 0.27007
ELA: 1.7015 Jacobi stage: ITER: 6839 RES: 0.25355
ELA: 1.8016 Jacobi stage: ITER: 7245 RES: 0.23795
ELA: 1.9017 Jacobi stage: ITER: 7653 RES: 0.22336
ELA: 2.0018 Jacobi stage: ITER: 7847 RES: 0.21682
ELA: 2.1019 Jacobi stage: ITER: 8019 RES: 0.21113
ELA: 2.202 Jacobi stage: ITER: 8199 RES: 0.20533
ELA: 2.3021 Jacobi stage: ITER: 8379 RES: 0.1997
ELA: 2.4022 Jacobi stage: ITER: 8555 RES: 0.19435
ELA: 2.5023 Jacobi stage: ITER: 8731 RES: 0.18915
ELA: 2.6023 Jacobi stage: ITER: 8907 RES: 0.18408
ELA: 2.7024 Jacobi stage: ITER: 9087 RES: 0.17905
ELA: 2.8025 Jacobi stage: ITER: 9259 RES: 0.17436
ELA: 2.9026 Jacobi stage: ITER: 9435 RES: 0.1697
ELA: 3.0027 Jacobi stage: ITER: 9609 RES: 0.16516
ELA: 3.1028 Jacobi stage: ITER: 9789 RES: 0.16065
ELA: 3.2029 Jacobi stage: ITER: 10063 RES: 0.15407
ELA: 3.303 Jacobi stage: ITER: 10459 RES: 0.14496
ELA: 3.4031 Jacobi stage: ITER: 10871 RES: 0.13606
ELA: 3.5032 Jacobi stage: ITER: 11265 RES: 0.12803
ELA: 3.6032 Jacobi stage: ITER: 11679 RES: 0.12017
ELA: 3.7033 Jacobi stage: ITER: 12067 RES: 0.11322
ELA: 3.8034 Jacobi stage: ITER: 12451 RES: 0.10673
ELA: 3.9035 Jacobi stage: ITER: 12827 RES: 0.10074
ELA: 4.0036 Jacobi stage: ITER: 13211 RES: 0.094967
ELA: 4.1037 Jacobi stage: ITER: 13599 RES: 0.089472
ELA: 4.2038 Jacobi stage: ITER: 13983 RES: 0.084347
ELA: 4.3038 Jacobi stage: ITER: 14369 RES: 0.079466
ELA: 4.4039 Jacobi stage: ITER: 14755 RES: 0.074914
ELA: 4.504 Jacobi stage: ITER: 15139 RES: 0.070623
ELA: 4.604 Jacobi stage: ITER: 15525 RES: 0.066537
ELA: 4.7041 Jacobi stage: ITER: 15911 RES: 0.062726
ELA: 4.8041 Jacobi stage: ITER: 16295 RES: 0.059133
ELA: 4.9042 Jacobi stage: ITER: 16681 RES: 0.055712
ELA: 5.0042 Jacobi stage: ITER: 17067 RES: 0.052521
ELA: 5.1043 Jacobi stage: ITER: 17451 RES: 0.049512
ELA: 5.2043 Jacobi stage: ITER: 17835 RES: 0.046676
ELA: 5.3043 Jacobi stage: ITER: 18221 RES: 0.043976
ELA: 5.4045 Jacobi stage: ITER: 18609 RES: 0.041432
ELA: 5.5045 Jacobi stage: ITER: 18995 RES: 0.039058
ELA: 5.6045 Jacobi stage: ITER: 19379 RES: 0.036821
ELA: 5.7046 Jacobi stage: ITER: 19763 RES: 0.034712
ELA: 5.8046 Jacobi stage: ITER: 20149 RES: 0.032704
ELA: 5.9046 Jacobi stage: ITER: 20535 RES: 0.030831
ELA: 6.0047 Jacobi stage: ITER: 20919 RES: 0.029065
ELA: 6.1048 Jacobi stage: ITER: 21305 RES: 0.027383
ELA: 6.2048 Jacobi stage: ITER: 21691 RES: 0.025815
ELA: 6.3049 Jacobi stage: ITER: 22075 RES: 0.024336
ELA: 6.4049 Jacobi stage: ITER: 22461 RES: 0.022928
ELA: 6.5049 Jacobi stage: ITER: 22847 RES: 0.021615
ELA: 6.605 Jacobi stage: ITER: 23231 RES: 0.020377
ELA: 6.705 Jacobi stage: ITER: 23617 RES: 0.019198
ELA: 6.805 Jacobi stage: ITER: 24001 RES: 0.018098
ELA: 6.9051 Jacobi stage: ITER: 24387 RES: 0.017062
ELA: 7.0051 Jacobi stage: ITER: 24771 RES: 0.016084
ELA: 7.1052 Jacobi stage: ITER: 25157 RES: 0.015154
ELA: 7.2052 Jacobi stage: ITER: 25543 RES: 0.014286
ELA: 7.3052 Jacobi stage: ITER: 25927 RES: 0.013468
ELA: 7.4053 Jacobi stage: ITER: 26313 RES: 0.012688
ELA: 7.5053 Jacobi stage: ITER: 26699 RES: 0.011962
ELA: 7.6054 Jacobi stage: ITER: 27083 RES: 0.011277
ELA: 7.7054 Jacobi stage: ITER: 27469 RES: 0.010624
ELA: 7.8054 Jacobi stage: ITER: 27855 RES: 0.010016
ELA: 7.9055 Jacobi stage: ITER: 28239 RES: 0.009442
ELA: 8.0055 Jacobi stage: ITER: 28625 RES: 0.0088957
ELA: 8.1055 Jacobi stage: ITER: 29011 RES: 0.0083862
ELA: 8.2056 Jacobi stage: ITER: 29395 RES: 0.0079058
ELA: 8.3056 Jacobi stage: ITER: 29781 RES: 0.0074484
ELA: 8.4056 Jacobi stage: ITER: 30167 RES: 0.0070218
ELA: 8.5057 Jacobi stage: ITER: 30551 RES: 0.0066196
ELA: 8.6057 Jacobi stage: ITER: 30935 RES: 0.0062405
ELA: 8.7058 Jacobi stage: ITER: 31321 RES: 0.0058794
ELA: 8.8058 Jacobi stage: ITER: 31707 RES: 0.0055427
ELA: 8.9058 Jacobi stage: ITER: 32091 RES: 0.0052252
ELA: 9.0059 Jacobi stage: ITER: 32477 RES: 0.0049229
ELA: 9.1059 Jacobi stage: ITER: 32863 RES: 0.0046409
ELA: 9.206 Jacobi stage: ITER: 33247 RES: 0.0043751
ELA: 9.306 Jacobi stage: ITER: 33633 RES: 0.004122
ELA: 9.406 Jacobi stage: ITER: 34019 RES: 0.0038859
ELA: 9.5061 Jacobi stage: ITER: 34403 RES: 0.0036633
ELA: 9.6061 Jacobi stage: ITER: 34789 RES: 0.0034514
ELA: 9.7061 Jacobi stage: ITER: 35175 RES: 0.0032537
ELA: 9.8062 Jacobi stage: ITER: 35559 RES: 0.0030673
ELA: 9.9062 Jacobi stage: ITER: 35945 RES: 0.0028898
ELA: 10.006 Jacobi stage: ITER: 36331 RES: 0.0027243
ELA: 10.106 Jacobi stage: ITER: 36715 RES: 0.0025683
ELA: 10.206 Jacobi stage: ITER: 37101 RES: 0.0024197
ELA: 10.306 Jacobi stage: ITER: 37487 RES: 0.0022811
ELA: 10.406 Jacobi stage: ITER: 37873 RES: 0.0021491
ELA: 10.506 Jacobi stage: ITER: 38259 RES: 0.002026
ELA: 10.606 Jacobi stage: ITER: 38643 RES: 0.00191
ELA: 10.707 Jacobi stage: ITER: 39029 RES: 0.0017995
ELA: 10.807 Jacobi stage: ITER: 39415 RES: 0.0016964
ELA: 10.907 Jacobi stage: ITER: 39799 RES: 0.0015992
ELA: 11.007 Jacobi stage: ITER: 40185 RES: 0.0015067
ELA: 11.107 Jacobi stage: ITER: 40571 RES: 0.0014204
ELA: 11.207 Jacobi stage: ITER: 40955 RES: 0.0013391
ELA: 11.307 Jacobi stage: ITER: 41341 RES: 0.0012616
ELA: 11.407 Jacobi stage: ITER: 41727 RES: 0.0011893
ELA: 11.507 Jacobi stage: ITER: 42111 RES: 0.0011212
ELA: 11.607 Jacobi stage: ITER: 42497 RES: 0.0010563
ELA: 11.707 Jacobi stage: ITER: 42883 RES: 0.00099583
ELA: 11.807 Jacobi stage: ITER: 43269 RES: 0.00093822
ELA: 11.907 Jacobi stage: ITER: 43655 RES: 0.00088448
ELA: 12.007 Jacobi stage: ITER: 44039 RES: 0.00083382
ELA: 12.107 Jacobi stage: ITER: 44425 RES: 0.00078558
ELA: 12.207 Jacobi stage: ITER: 44811 RES: 0.00074058
ELA: 12.307 Jacobi stage: ITER: 45195 RES: 0.00069816
ELA: 12.407 Jacobi stage: ITER: 45581 RES: 0.00065777
ELA: 12.507 Jacobi stage: ITER: 45967 RES: 0.00062009
ELA: 12.607 Jacobi stage: ITER: 46351 RES: 0.00058458
ELA: 12.707 Jacobi stage: ITER: 46737 RES: 0.00055076
ELA: 12.807 Jacobi stage: ITER: 47123 RES: 0.00051921
ELA: 12.907 Jacobi stage: ITER: 47507 RES: 0.00048947
ELA: 13.007 Jacobi stage: ITER: 47893 RES: 0.00046115
ELA: 13.108 Jacobi stage: ITER: 48279 RES: 0.00043474
ELA: 13.208 Jacobi stage: ITER: 48663 RES: 0.00040984
ELA: 13.308 Jacobi stage: ITER: 49049 RES: 0.00038613
ELA: 13.408 Jacobi stage: ITER: 49435 RES: 0.00036401
ELA: 13.508 Jacobi stage: ITER: 49819 RES: 0.00034316
ELA: 13.608 Jacobi stage: ITER: 50205 RES: 0.00032331
ELA: 13.708 Jacobi stage: ITER: 50591 RES: 0.00030479
ELA: 13.808 Jacobi stage: ITER: 50975 RES: 0.00028733
ELA: 13.908 Jacobi stage: ITER: 51361 RES: 0.00027071
ELA: 14.008 Jacobi stage: ITER: 51747 RES: 0.0002552
ELA: 14.108 Jacobi stage: ITER: 52131 RES: 0.00024058
ELA: 14.208 Jacobi stage: ITER: 52517 RES: 0.00022666
ELA: 14.308 Jacobi stage: ITER: 52903 RES: 0.00021368
ELA: 14.408 Jacobi stage: ITER: 53287 RES: 0.00020144
ELA: 14.508 Jacobi stage: ITER: 53675 RES: 0.00018979
ELA: 14.608 Jacobi stage: ITER: 54059 RES: 0.00017892
ELA: 14.708 Jacobi stage: ITER: 54445 RES: 0.00016857
ELA: 14.808 Jacobi stage: ITER: 54831 RES: 0.00015891
ELA: 14.908 Jacobi stage: ITER: 55215 RES: 0.00014981
ELA: 15.008 Jacobi stage: ITER: 55601 RES: 0.00014114
ELA: 15.108 Jacobi stage: ITER: 55987 RES: 0.00013306
ELA: 15.208 Jacobi stage: ITER: 56371 RES: 0.00012544
ELA: 15.308 Jacobi stage: ITER: 56757 RES: 0.00011818
ELA: 15.408 Jacobi stage: ITER: 57143 RES: 0.00011141
ELA: 15.509 Jacobi stage: ITER: 57527 RES: 0.00010503
ELA: 15.609 Jacobi stage: ITER: 57913 RES: 9.8952e-05
ELA: 15.709 Jacobi stage: ITER: 58299 RES: 9.3285e-05
ELA: 15.809 Jacobi stage: ITER: 58683 RES: 8.7942e-05
ELA: 15.909 Jacobi stage: ITER: 59069 RES: 8.2854e-05
ELA: 16.009 Jacobi stage: ITER: 59455 RES: 7.8108e-05
ELA: 16.243 Jacobi stage: ITER: 60357 RES: 6.7982e-05
ELA: 16.343 Jacobi stage: ITER: 60743 RES: 6.4088e-05
ELA: 16.443 Jacobi stage: ITER: 61127 RES: 6.0417e-05
ELA: 16.543 Jacobi stage: ITER: 61513 RES: 5.6921e-05
ELA: 16.643 Jacobi stage: ITER: 61899 RES: 5.3661e-05
ELA: 16.743 Jacobi stage: ITER: 62283 RES: 5.0588e-05
ELA: 16.843 Jacobi stage: ITER: 62669 RES: 4.7661e-05
ELA: 16.943 Jacobi stage: ITER: 63055 RES: 4.4931e-05
ELA: 17.043 Jacobi stage: ITER: 63441 RES: 4.2331e-05
ELA: 17.143 Jacobi stage: ITER: 63827 RES: 3.9907e-05
ELA: 17.243 Jacobi stage: ITER: 64211 RES: 3.7621e-05
ELA: 17.344 Jacobi stage: ITER: 64597 RES: 3.5444e-05
ELA: 17.444 Jacobi stage: ITER: 64983 RES: 3.3414e-05
ELA: 17.544 Jacobi stage: ITER: 65367 RES: 3.15e-05
ELA: 17.644 Jacobi stage: ITER: 65895 RES: 2.9046e-05
ELA: 17.744 Jacobi stage: ITER: 66491 RES: 2.6505e-05
ELA: 17.844 Jacobi stage: ITER: 67071 RES: 2.4246e-05
ELA: 17.944 Jacobi stage: ITER: 67651 RES: 2.218e-05
ELA: 18.044 Jacobi stage: ITER: 68231 RES: 2.0289e-05
ELA: 18.144 Jacobi stage: ITER: 68811 RES: 1.856e-05
ELA: 18.244 Jacobi stage: ITER: 69391 RES: 1.6978e-05
ELA: 18.344 Jacobi stage: ITER: 69971 RES: 1.5531e-05
ELA: 18.444 Jacobi stage: ITER: 70549 RES: 1.4207e-05
ELA: 18.544 Jacobi stage: ITER: 71129 RES: 1.2996e-05
ELA: 18.644 Jacobi stage: ITER: 71709 RES: 1.1888e-05
ELA: 18.744 Jacobi stage: ITER: 72289 RES: 1.0875e-05
ELA: 18.844 Jacobi stage: ITER: 72867 RES: 9.9541e-06
ELA: 18.944 Jacobi stage: ITER: 73447 RES: 9.1057e-06
ELA: 19.044 Jacobi stage: ITER: 74027 RES: 8.3295e-06
ELA: 19.144 Jacobi stage: ITER: 74607 RES: 7.6196e-06
ELA: 19.244 Jacobi stage: ITER: 75199 RES: 6.9573e-06
ELA: 19.344 Jacobi stage: ITER: 75777 RES: 6.3643e-06
ELA: 19.444 Jacobi stage: ITER: 76357 RES: 5.8218e-06
ELA: 19.544 Jacobi stage: ITER: 76937 RES: 5.3256e-06
ELA: 19.644 Jacobi stage: ITER: 77517 RES: 4.8716e-06
ELA: 19.744 Jacobi stage: ITER: 78095 RES: 4.4591e-06
ELA: 19.844 Jacobi stage: ITER: 78675 RES: 4.0791e-06
ELA: 19.944 Jacobi stage: ITER: 79255 RES: 3.7314e-06
ELA: 20.044 Jacobi stage: ITER: 79835 RES: 3.4133e-06
ELA: 20.145 Jacobi stage: ITER: 80415 RES: 3.1224e-06
ELA: 20.245 Jacobi stage: ITER: 80995 RES: 2.8563e-06
ELA: 20.345 Jacobi stage: ITER: 81575 RES: 2.6128e-06
ELA: 20.445 Jacobi stage: ITER: 82155 RES: 2.3901e-06
ELA: 20.545 Jacobi stage: ITER: 82735 RES: 2.1864e-06
ELA: 20.645 Jacobi stage: ITER: 83315 RES: 2e-06
ELA: 20.745 Jacobi stage: ITER: 83893 RES: 1.8295e-06
ELA: 20.845 Jacobi stage: ITER: 84473 RES: 1.6736e-06
ELA: 20.945 Jacobi stage: ITER: 85053 RES: 1.531e-06
ELA: 21.045 Jacobi stage: ITER: 85643 RES: 1.3987e-06
ELA: 21.145 Jacobi stage: ITER: 86223 RES: 1.2795e-06
ELA: 21.245 Jacobi stage: ITER: 86803 RES: 1.1705e-06
ELA: 21.345 Jacobi stage: ITER: 87383 RES: 1.0707e-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 );
26 matrix_ptr->setDimensions( size, size );
27 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
28
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 );
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;
61
62 /***
63 * Solve the linear system using diagonal (Jacobi) preconditioner.
64 */
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 );
75}
76
77int
78main( int argc, char* argv[] )
79{
81 iterativeLinearSolverExample< TNL::Devices::Sequential >();
82
83#ifdef __CUDACC__
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 in the solver setup. Similarly to the matrix of the linear system, the preconditioner needs to be passed to the solver as a std::shared_ptr. When the preconditioner object is created, we have to initialize it using the update method, which has to be called everytime the matrix of the linear system changes. This is important, for example, when solving time-dependent PDEs, but it does not happen in this example. Finally, we need to connect the solver with the preconditioner using the setPreconditioner method.
The result looks as follows:
Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Choosing the solver and preconditioner type at runtime
When developing a numerical solver, one often has to search for a combination of various methods and algorithms that fit given requirements the best. To make this easier, TNL provides the functions TNL::Solvers::getLinearSolver and TNL::Solvers::getPreconditioner for selecting the linear solver and preconditioner at runtime. The following example shows how to use these functions:
1#include <iostream>
2#include <memory>
3#include <TNL/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Host.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/LinearSolverTypeResolver.h>
7
8template< typename Device >
9void
10iterativeLinearSolverExample()
11{
12 /***
13 * Set the following matrix (dots represent zero matrix elements):
14 *
15 * / 2.5 -1 . . . \
16 * | -1 2.5 -1 . . |
17 * | . -1 2.5 -1. . |
18 * | . . -1 2.5 -1 |
19 * \ . . . -1 2.5 /
20 */
23 const int size( 5 );
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
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 );
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;
60
61 /***
62 * Solve the linear system using diagonal (Jacobi) preconditioner.
63 */
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 );
72}
73
74int
75main( int argc, char* argv[] )
76{
78 iterativeLinearSolverExample< TNL::Devices::Sequential >();
79
80#ifdef __CUDACC__
82 iterativeLinearSolverExample< TNL::Devices::Cuda >();
83#endif
84}
std::shared_ptr< Linear::Preconditioners::Preconditioner< MatrixType > > getPreconditioner(const std::string &name)
Function returning shared pointer with linear preconditioner given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:144
std::shared_ptr< Linear::LinearSolver< MatrixType > > getLinearSolver(const std::string &name)
Function returning shared pointer with linear solver given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:87
We still stay with the same problem and the only changes are in the solver setup. We first use TNL::Solvers::getLinearSolver to get a shared pointer holding the solver and then TNL::Solvers::getPreconditioner to get a shared pointer holding the preconditioner. The rest of the code is the same as in the previous examples with the only difference that we work with the pointer solver_ptr instead of the direct instance solver of the solver type.
The result looks as follows:
Solving linear system on host:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Solving linear system on CUDA device:
Row: 0 -> 0:2.5 1:-1
Row: 1 -> 0:-1 1:2.5 2:-1
Row: 2 -> 1:-1 2:2.5 3:-1
Row: 3 -> 2:-1 3:2.5 4:-1
Row: 4 -> 3:-1 4:2.5
Vector b = [ 1.5, 0.5, 0.5, 0.5, 1.5 ]
Vector x = [ 1, 1, 1, 1, 1 ]
Generated on for Template Numerical Library by
1.16.1