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: 46727 RES: 0.00055143
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: 4961 RES: 0.34123
Jacobi stage: ITER: 6747 RES: 0.25722
Jacobi stage: ITER: 7203 RES: 0.23959
Jacobi stage: ITER: 7667 RES: 0.22295
Jacobi stage: ITER: 8131 RES: 0.2075
Jacobi stage: ITER: 8575 RES: 0.19375
Jacobi stage: ITER: 9031 RES: 0.1806
Jacobi stage: ITER: 9497 RES: 0.16804
Jacobi stage: ITER: 9945 RES: 0.15684
Jacobi stage: ITER: 10405 RES: 0.14613
Jacobi stage: ITER: 10871 RES: 0.13606
Jacobi stage: ITER: 11331 RES: 0.12677
Jacobi stage: ITER: 11781 RES: 0.11827
Jacobi stage: ITER: 12219 RES: 0.1106
Jacobi stage: ITER: 12689 RES: 0.10287
Jacobi stage: ITER: 13155 RES: 0.095788
Jacobi stage: ITER: 13615 RES: 0.089252
Jacobi stage: ITER: 14069 RES: 0.083214
Jacobi stage: ITER: 14531 RES: 0.077537
Jacobi stage: ITER: 14909 RES: 0.07314
Jacobi stage: ITER: 15183 RES: 0.070147
Jacobi stage: ITER: 15455 RES: 0.067277
Jacobi stage: ITER: 15737 RES: 0.064405
Jacobi stage: ITER: 16027 RES: 0.061618
Jacobi stage: ITER: 16307 RES: 0.059024
Jacobi stage: ITER: 16591 RES: 0.056505
Jacobi stage: ITER: 16871 RES: 0.054126
Jacobi stage: ITER: 17245 RES: 0.051088
Jacobi stage: ITER: 17695 RES: 0.047691
Jacobi stage: ITER: 18147 RES: 0.044492
Jacobi stage: ITER: 18603 RES: 0.041483
Jacobi stage: ITER: 19069 RES: 0.038605
Jacobi stage: ITER: 19535 RES: 0.035949
Jacobi stage: ITER: 19999 RES: 0.033476
Jacobi stage: ITER: 20443 RES: 0.03127
Jacobi stage: ITER: 20883 RES: 0.029226
Jacobi stage: ITER: 21351 RES: 0.027199
Jacobi stage: ITER: 21817 RES: 0.025312
Jacobi stage: ITER: 22283 RES: 0.023571
Jacobi stage: ITER: 22751 RES: 0.021936
Jacobi stage: ITER: 23227 RES: 0.020389
Jacobi stage: ITER: 23639 RES: 0.019139
Jacobi stage: ITER: 24099 RES: 0.017833
Jacobi stage: ITER: 24565 RES: 0.016596
Jacobi stage: ITER: 25027 RES: 0.015464
Jacobi stage: ITER: 25453 RES: 0.01448
Jacobi stage: ITER: 25893 RES: 0.013534
Jacobi stage: ITER: 26333 RES: 0.01265
Jacobi stage: ITER: 26773 RES: 0.011823
Jacobi stage: ITER: 27213 RES: 0.01105
Jacobi stage: ITER: 27653 RES: 0.010328
Jacobi stage: ITER: 28093 RES: 0.0096531
Jacobi stage: ITER: 28533 RES: 0.0090223
Jacobi stage: ITER: 28973 RES: 0.0084327
Jacobi stage: ITER: 29413 RES: 0.0078816
Jacobi stage: ITER: 29853 RES: 0.0073665
Jacobi stage: ITER: 30293 RES: 0.0068851
Jacobi stage: ITER: 30733 RES: 0.0064352
Jacobi stage: ITER: 31173 RES: 0.0060146
Jacobi stage: ITER: 31613 RES: 0.0056215
Jacobi stage: ITER: 32053 RES: 0.0052542
Jacobi stage: ITER: 32491 RES: 0.0049138
Jacobi stage: ITER: 32931 RES: 0.0045927
Jacobi stage: ITER: 33371 RES: 0.0042926
Jacobi stage: ITER: 33811 RES: 0.004012
Jacobi stage: ITER: 34249 RES: 0.0037498
Jacobi stage: ITER: 34689 RES: 0.0035048
Jacobi stage: ITER: 35129 RES: 0.0032757
Jacobi stage: ITER: 35567 RES: 0.0030635
Jacobi stage: ITER: 36007 RES: 0.0028633
Jacobi stage: ITER: 36447 RES: 0.0026762
Jacobi stage: ITER: 36887 RES: 0.0025013
Jacobi stage: ITER: 37327 RES: 0.0023379
Jacobi stage: ITER: 37765 RES: 0.0021851
Jacobi stage: ITER: 38205 RES: 0.0020423
Jacobi stage: ITER: 38647 RES: 0.0019088
Jacobi stage: ITER: 39087 RES: 0.0017841
Jacobi stage: ITER: 39527 RES: 0.0016675
Jacobi stage: ITER: 39967 RES: 0.0015585
Jacobi stage: ITER: 40409 RES: 0.0014558
Jacobi stage: ITER: 40851 RES: 0.0013606
Jacobi stage: ITER: 41291 RES: 0.0012717
Jacobi stage: ITER: 41731 RES: 0.0011886
Jacobi stage: ITER: 42171 RES: 0.0011109
Jacobi stage: ITER: 42611 RES: 0.0010383
Jacobi stage: ITER: 43051 RES: 0.00097046
Jacobi stage: ITER: 43491 RES: 0.00090704
Jacobi stage: ITER: 43931 RES: 0.00084777
Jacobi stage: ITER: 44373 RES: 0.00079188
Jacobi stage: ITER: 44811 RES: 0.00074058
Jacobi stage: ITER: 45253 RES: 0.00069176
Jacobi stage: ITER: 45693 RES: 0.00064655
Jacobi stage: ITER: 46135 RES: 0.0006043
Jacobi stage: ITER: 46573 RES: 0.00056481
Jacobi stage: ITER: 47013 RES: 0.00052789
Jacobi stage: ITER: 47455 RES: 0.0004934
Jacobi stage: ITER: 47895 RES: 0.00046115
Jacobi stage: ITER: 48331 RES: 0.00043128
Jacobi stage: ITER: 48771 RES: 0.00040309
Jacobi stage: ITER: 49211 RES: 0.00037675
Jacobi stage: ITER: 49651 RES: 0.00035213
Jacobi stage: ITER: 50089 RES: 0.00032912
Jacobi stage: ITER: 50527 RES: 0.0003078
Jacobi stage: ITER: 50967 RES: 0.00028768
Jacobi stage: ITER: 51407 RES: 0.00026888
Jacobi stage: ITER: 51845 RES: 0.00025131
Jacobi stage: ITER: 52285 RES: 0.00023489
Jacobi stage: ITER: 52723 RES: 0.00021967
Jacobi stage: ITER: 53163 RES: 0.00020532
Jacobi stage: ITER: 53601 RES: 0.0001919
Jacobi stage: ITER: 54039 RES: 0.00017947
Jacobi stage: ITER: 54479 RES: 0.00016774
Jacobi stage: ITER: 54915 RES: 0.00015687
Jacobi stage: ITER: 55355 RES: 0.00014662
Jacobi stage: ITER: 55795 RES: 0.00013704
Jacobi stage: ITER: 56233 RES: 0.00012808
Jacobi stage: ITER: 56671 RES: 0.00011979
Jacobi stage: ITER: 57109 RES: 0.00011196
Jacobi stage: ITER: 57547 RES: 0.00010471
Jacobi stage: ITER: 57987 RES: 9.7864e-05
Jacobi stage: ITER: 58425 RES: 9.1469e-05
Jacobi stage: ITER: 58863 RES: 8.5543e-05
Jacobi stage: ITER: 59303 RES: 7.9953e-05
Jacobi stage: ITER: 59741 RES: 7.4728e-05
Jacobi stage: ITER: 60179 RES: 6.9887e-05
Jacobi stage: ITER: 60619 RES: 6.532e-05
Jacobi stage: ITER: 61057 RES: 6.1051e-05
Jacobi stage: ITER: 61495 RES: 5.7097e-05
Jacobi stage: ITER: 61935 RES: 5.3365e-05
Jacobi stage: ITER: 62373 RES: 4.9878e-05
Jacobi stage: ITER: 62811 RES: 4.6647e-05
Jacobi stage: ITER: 63251 RES: 4.3598e-05
Jacobi stage: ITER: 63689 RES: 4.0749e-05
Jacobi stage: ITER: 64127 RES: 3.811e-05
Jacobi stage: ITER: 64567 RES: 3.5619e-05
Jacobi stage: ITER: 65003 RES: 3.3312e-05
Jacobi stage: ITER: 65443 RES: 3.1135e-05
Jacobi stage: ITER: 65881 RES: 2.91e-05
Jacobi stage: ITER: 66319 RES: 2.7215e-05
Jacobi stage: ITER: 66759 RES: 2.5437e-05
Jacobi stage: ITER: 67195 RES: 2.3789e-05
Jacobi stage: ITER: 67635 RES: 2.2234e-05
Jacobi stage: ITER: 68073 RES: 2.0781e-05
Jacobi stage: ITER: 68503 RES: 1.9459e-05
Jacobi stage: ITER: 68943 RES: 1.8187e-05
Jacobi stage: ITER: 69381 RES: 1.6999e-05
Jacobi stage: ITER: 69819 RES: 1.5898e-05
Jacobi stage: ITER: 70259 RES: 1.4859e-05
Jacobi stage: ITER: 70697 RES: 1.3888e-05
Jacobi stage: ITER: 71135 RES: 1.2988e-05
Jacobi stage: ITER: 71575 RES: 1.2139e-05
Jacobi stage: ITER: 72013 RES: 1.1346e-05
Jacobi stage: ITER: 72451 RES: 1.0611e-05
Jacobi stage: ITER: 72891 RES: 9.9175e-06
Jacobi stage: ITER: 73329 RES: 9.2694e-06
Jacobi stage: ITER: 73767 RES: 8.6689e-06
Jacobi stage: ITER: 74207 RES: 8.1024e-06
Jacobi stage: ITER: 74645 RES: 7.5729e-06
Jacobi stage: ITER: 75083 RES: 7.0824e-06
Jacobi stage: ITER: 75523 RES: 6.6195e-06
Jacobi stage: ITER: 75961 RES: 6.1869e-06
Jacobi stage: ITER: 76399 RES: 5.7861e-06
Jacobi stage: ITER: 76839 RES: 5.408e-06
Jacobi stage: ITER: 77275 RES: 5.0577e-06
Jacobi stage: ITER: 77715 RES: 4.7272e-06
Jacobi stage: ITER: 78153 RES: 4.4182e-06
Jacobi stage: ITER: 78591 RES: 4.132e-06
Jacobi stage: ITER: 79029 RES: 3.862e-06
Jacobi stage: ITER: 79467 RES: 3.6118e-06
Jacobi stage: ITER: 79907 RES: 3.3758e-06
Jacobi stage: ITER: 80345 RES: 3.1552e-06
Jacobi stage: ITER: 80783 RES: 2.9508e-06
Jacobi stage: ITER: 81223 RES: 2.758e-06
Jacobi stage: ITER: 81661 RES: 2.5777e-06
Jacobi stage: ITER: 82099 RES: 2.4107e-06
Jacobi stage: ITER: 82537 RES: 2.2532e-06
Jacobi stage: ITER: 82975 RES: 2.1072e-06
Jacobi stage: ITER: 83415 RES: 1.9695e-06
Jacobi stage: ITER: 83853 RES: 1.8408e-06
Jacobi stage: ITER: 84291 RES: 1.7216e-06
Jacobi stage: ITER: 84729 RES: 1.6091e-06
Jacobi stage: ITER: 85167 RES: 1.5048e-06
Jacobi stage: ITER: 85605 RES: 1.4065e-06
Jacobi stage: ITER: 86043 RES: 1.3154e-06
Jacobi stage: ITER: 86483 RES: 1.2294e-06
Jacobi stage: ITER: 86921 RES: 1.1491e-06
Jacobi stage: ITER: 87359 RES: 1.0746e-06
Jacobi stage: ITER: 87797 RES: 1.0044e-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.447e-0
ELA: 0.1001 Jacobi stage: ITER: 48341 RES: 0.00043049
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.177e-0
ELA: 0.10011 Jacobi stage: ITER: 459 RES: 0.86172
ELA: 0.2002 Jacobi stage: ITER: 923 RES: 0.75394
ELA: 0.3003 Jacobi stage: ITER: 1383 RES: 0.67008
ELA: 0.40041 Jacobi stage: ITER: 1831 RES: 0.60406
ELA: 0.50059 Jacobi stage: ITER: 2279 RES: 0.54923
ELA: 0.60072 Jacobi stage: ITER: 2729 RES: 0.5023
ELA: 0.70081 Jacobi stage: ITER: 3185 RES: 0.46141
ELA: 0.8009 Jacobi stage: ITER: 3629 RES: 0.42639
ELA: 0.90099 Jacobi stage: ITER: 4091 RES: 0.39405
ELA: 1.0011 Jacobi stage: ITER: 4553 RES: 0.36474
ELA: 1.1012 Jacobi stage: ITER: 5001 RES: 0.33903
ELA: 1.2013 Jacobi stage: ITER: 5435 RES: 0.3163
ELA: 1.3014 Jacobi stage: ITER: 5895 RES: 0.29402
ELA: 1.4014 Jacobi stage: ITER: 6353 RES: 0.27348
ELA: 1.5015 Jacobi stage: ITER: 6813 RES: 0.2545
ELA: 1.6016 Jacobi stage: ITER: 7271 RES: 0.23707
ELA: 1.7017 Jacobi stage: ITER: 7727 RES: 0.22088
ELA: 1.8018 Jacobi stage: ITER: 8175 RES: 0.2061
ELA: 1.9019 Jacobi stage: ITER: 8447 RES: 0.19762
ELA: 2.002 Jacobi stage: ITER: 8707 RES: 0.18985
ELA: 2.1021 Jacobi stage: ITER: 8979 RES: 0.18205
ELA: 2.2022 Jacobi stage: ITER: 9263 RES: 0.17426
ELA: 2.3023 Jacobi stage: ITER: 9543 RES: 0.1669
ELA: 2.4024 Jacobi stage: ITER: 9823 RES: 0.15986
ELA: 2.5025 Jacobi stage: ITER: 10103 RES: 0.15312
ELA: 2.6025 Jacobi stage: ITER: 10411 RES: 0.14604
ELA: 2.7026 Jacobi stage: ITER: 10861 RES: 0.13623
ELA: 2.8027 Jacobi stage: ITER: 11311 RES: 0.12716
ELA: 2.9028 Jacobi stage: ITER: 11765 RES: 0.11856
ELA: 3.0029 Jacobi stage: ITER: 12219 RES: 0.1106
ELA: 3.1029 Jacobi stage: ITER: 12683 RES: 0.10299
ELA: 3.203 Jacobi stage: ITER: 13147 RES: 0.095906
ELA: 3.3031 Jacobi stage: ITER: 13607 RES: 0.089362
ELA: 3.4032 Jacobi stage: ITER: 14023 RES: 0.08383
ELA: 3.5033 Jacobi stage: ITER: 14473 RES: 0.078207
ELA: 3.6034 Jacobi stage: ITER: 14935 RES: 0.072871
ELA: 3.7035 Jacobi stage: ITER: 15405 RES: 0.067775
ELA: 3.8035 Jacobi stage: ITER: 15867 RES: 0.063151
ELA: 3.9036 Jacobi stage: ITER: 16339 RES: 0.058735
ELA: 4.0037 Jacobi stage: ITER: 16767 RES: 0.054997
ELA: 4.1038 Jacobi stage: ITER: 17187 RES: 0.051561
ELA: 4.2038 Jacobi stage: ITER: 17653 RES: 0.047985
ELA: 4.3039 Jacobi stage: ITER: 18119 RES: 0.044684
ELA: 4.404 Jacobi stage: ITER: 18539 RES: 0.041892
ELA: 4.5041 Jacobi stage: ITER: 18987 RES: 0.039106
ELA: 4.6042 Jacobi stage: ITER: 19427 RES: 0.036551
ELA: 4.7043 Jacobi stage: ITER: 19867 RES: 0.034162
ELA: 4.8044 Jacobi stage: ITER: 20307 RES: 0.03193
ELA: 4.9045 Jacobi stage: ITER: 20747 RES: 0.029843
ELA: 5.0046 Jacobi stage: ITER: 21187 RES: 0.027893
ELA: 5.1046 Jacobi stage: ITER: 21627 RES: 0.02607
ELA: 5.2047 Jacobi stage: ITER: 22067 RES: 0.024366
ELA: 5.3048 Jacobi stage: ITER: 22507 RES: 0.022774
ELA: 5.4049 Jacobi stage: ITER: 22947 RES: 0.021285
ELA: 5.505 Jacobi stage: ITER: 23387 RES: 0.019894
ELA: 5.6051 Jacobi stage: ITER: 23827 RES: 0.018594
ELA: 5.7052 Jacobi stage: ITER: 24267 RES: 0.017379
ELA: 5.8052 Jacobi stage: ITER: 24707 RES: 0.016243
ELA: 5.9053 Jacobi stage: ITER: 25147 RES: 0.015182
ELA: 6.0054 Jacobi stage: ITER: 25587 RES: 0.01419
ELA: 6.1055 Jacobi stage: ITER: 26027 RES: 0.013262
ELA: 6.2056 Jacobi stage: ITER: 26465 RES: 0.012396
ELA: 6.3057 Jacobi stage: ITER: 26905 RES: 0.011586
ELA: 6.4058 Jacobi stage: ITER: 27343 RES: 0.010835
ELA: 6.5059 Jacobi stage: ITER: 27783 RES: 0.010127
ELA: 6.6059 Jacobi stage: ITER: 28223 RES: 0.0094652
ELA: 6.706 Jacobi stage: ITER: 28663 RES: 0.0088466
ELA: 6.8061 Jacobi stage: ITER: 29103 RES: 0.0082685
ELA: 6.9062 Jacobi stage: ITER: 29541 RES: 0.0077281
ELA: 7.0063 Jacobi stage: ITER: 29981 RES: 0.0072231
ELA: 7.1064 Jacobi stage: ITER: 30421 RES: 0.0067511
ELA: 7.2069 Jacobi stage: ITER: 30863 RES: 0.0063099
ELA: 7.3073 Jacobi stage: ITER: 31303 RES: 0.0058975
ELA: 7.4074 Jacobi stage: ITER: 31743 RES: 0.0055121
ELA: 7.5079 Jacobi stage: ITER: 32183 RES: 0.0051519
ELA: 7.6079 Jacobi stage: ITER: 32623 RES: 0.0048152
ELA: 7.7084 Jacobi stage: ITER: 33065 RES: 0.0044977
ELA: 7.8089 Jacobi stage: ITER: 33507 RES: 0.0042038
ELA: 7.9093 Jacobi stage: ITER: 33947 RES: 0.0039291
ELA: 8.0098 Jacobi stage: ITER: 34387 RES: 0.0036723
ELA: 8.1103 Jacobi stage: ITER: 34829 RES: 0.0034302
ELA: 8.2104 Jacobi stage: ITER: 35269 RES: 0.003206
ELA: 8.3109 Jacobi stage: ITER: 35711 RES: 0.0029965
ELA: 8.411 Jacobi stage: ITER: 36149 RES: 0.0028007
ELA: 8.511 Jacobi stage: ITER: 36589 RES: 0.0026177
ELA: 8.6114 Jacobi stage: ITER: 37031 RES: 0.0024466
ELA: 8.7115 Jacobi stage: ITER: 37469 RES: 0.0022867
ELA: 8.812 Jacobi stage: ITER: 37911 RES: 0.0021373
ELA: 8.9124 Jacobi stage: ITER: 38351 RES: 0.0019976
ELA: 9.0129 Jacobi stage: ITER: 38793 RES: 0.0018659
ELA: 9.1133 Jacobi stage: ITER: 39233 RES: 0.001744
ELA: 9.2138 Jacobi stage: ITER: 39675 RES: 0.00163
ELA: 9.3139 Jacobi stage: ITER: 40115 RES: 0.0015235
ELA: 9.414 Jacobi stage: ITER: 40555 RES: 0.0014239
ELA: 9.5144 Jacobi stage: ITER: 40995 RES: 0.0013309
ELA: 9.6145 Jacobi stage: ITER: 41433 RES: 0.0012439
ELA: 9.7146 Jacobi stage: ITER: 41871 RES: 0.0011633
ELA: 9.815 Jacobi stage: ITER: 42311 RES: 0.0010873
ELA: 9.9154 Jacobi stage: ITER: 42751 RES: 0.0010162
ELA: 10.016 Jacobi stage: ITER: 43191 RES: 0.00094982
ELA: 10.116 Jacobi stage: ITER: 43631 RES: 0.00088774
ELA: 10.217 Jacobi stage: ITER: 44071 RES: 0.00082973
ELA: 10.317 Jacobi stage: ITER: 44511 RES: 0.00077551
ELA: 10.417 Jacobi stage: ITER: 44949 RES: 0.00072483
ELA: 10.518 Jacobi stage: ITER: 45389 RES: 0.00067746
ELA: 10.618 Jacobi stage: ITER: 45827 RES: 0.00063357
ELA: 10.718 Jacobi stage: ITER: 46267 RES: 0.00059217
ELA: 10.819 Jacobi stage: ITER: 46707 RES: 0.00055347
ELA: 10.919 Jacobi stage: ITER: 47147 RES: 0.0005173
ELA: 11.019 Jacobi stage: ITER: 47585 RES: 0.00048349
ELA: 11.12 Jacobi stage: ITER: 48025 RES: 0.0004519
ELA: 11.22 Jacobi stage: ITER: 48463 RES: 0.00042262
ELA: 11.32 Jacobi stage: ITER: 48903 RES: 0.000395
ELA: 11.421 Jacobi stage: ITER: 49343 RES: 0.00036919
ELA: 11.521 Jacobi stage: ITER: 49783 RES: 0.00034506
ELA: 11.621 Jacobi stage: ITER: 50219 RES: 0.00032271
ELA: 11.722 Jacobi stage: ITER: 50659 RES: 0.00030162
ELA: 11.822 Jacobi stage: ITER: 51099 RES: 0.00028191
ELA: 11.922 Jacobi stage: ITER: 51541 RES: 0.00026332
ELA: 12.023 Jacobi stage: ITER: 51979 RES: 0.00024627
ELA: 12.123 Jacobi stage: ITER: 52419 RES: 0.00023017
ELA: 12.223 Jacobi stage: ITER: 52859 RES: 0.00021513
ELA: 12.324 Jacobi stage: ITER: 53301 RES: 0.00020095
ELA: 12.424 Jacobi stage: ITER: 53741 RES: 0.00018782
ELA: 12.525 Jacobi stage: ITER: 54181 RES: 0.00017554
ELA: 12.625 Jacobi stage: ITER: 54621 RES: 0.00016407
ELA: 12.726 Jacobi stage: ITER: 55061 RES: 0.00015335
ELA: 12.826 Jacobi stage: ITER: 55503 RES: 0.00014333
ELA: 12.927 Jacobi stage: ITER: 55941 RES: 0.00013396
ELA: 13.027 Jacobi stage: ITER: 56379 RES: 0.00012528
ELA: 13.127 Jacobi stage: ITER: 56819 RES: 0.0001171
ELA: 13.227 Jacobi stage: ITER: 57255 RES: 0.00010951
ELA: 13.327 Jacobi stage: ITER: 57695 RES: 0.00010235
ELA: 13.427 Jacobi stage: ITER: 58135 RES: 9.5664e-05
ELA: 13.528 Jacobi stage: ITER: 58575 RES: 8.9413e-05
ELA: 13.628 Jacobi stage: ITER: 59015 RES: 8.3569e-05
ELA: 13.729 Jacobi stage: ITER: 59455 RES: 7.8108e-05
ELA: 13.829 Jacobi stage: ITER: 59893 RES: 7.3004e-05
ELA: 13.929 Jacobi stage: ITER: 60333 RES: 6.8233e-05
ELA: 14.029 Jacobi stage: ITER: 60771 RES: 6.3813e-05
ELA: 14.13 Jacobi stage: ITER: 61211 RES: 5.9642e-05
ELA: 14.23 Jacobi stage: ITER: 61655 RES: 5.571e-05
ELA: 14.33 Jacobi stage: ITER: 62093 RES: 5.207e-05
ELA: 14.43 Jacobi stage: ITER: 62531 RES: 4.8697e-05
ELA: 14.53 Jacobi stage: ITER: 62971 RES: 4.5514e-05
ELA: 14.631 Jacobi stage: ITER: 63411 RES: 4.254e-05
ELA: 14.731 Jacobi stage: ITER: 63851 RES: 3.976e-05
ELA: 14.832 Jacobi stage: ITER: 64293 RES: 3.7139e-05
ELA: 14.932 Jacobi stage: ITER: 64733 RES: 3.4712e-05
ELA: 15.033 Jacobi stage: ITER: 65173 RES: 3.2443e-05
ELA: 15.133 Jacobi stage: ITER: 65615 RES: 3.0323e-05
ELA: 15.234 Jacobi stage: ITER: 66055 RES: 2.8341e-05
ELA: 15.334 Jacobi stage: ITER: 66495 RES: 2.6489e-05
ELA: 15.435 Jacobi stage: ITER: 66935 RES: 2.4758e-05
ELA: 15.535 Jacobi stage: ITER: 67375 RES: 2.314e-05
ELA: 15.635 Jacobi stage: ITER: 67815 RES: 2.1628e-05
ELA: 15.736 Jacobi stage: ITER: 68255 RES: 2.0214e-05
ELA: 15.836 Jacobi stage: ITER: 68695 RES: 1.8893e-05
ELA: 15.937 Jacobi stage: ITER: 69135 RES: 1.7659e-05
ELA: 16.037 Jacobi stage: ITER: 69573 RES: 1.6505e-05
ELA: 16.137 Jacobi stage: ITER: 70011 RES: 1.5436e-05
ELA: 16.237 Jacobi stage: ITER: 70451 RES: 1.4427e-05
ELA: 16.338 Jacobi stage: ITER: 70891 RES: 1.3484e-05
ELA: 16.438 Jacobi stage: ITER: 71331 RES: 1.2603e-05
ELA: 16.539 Jacobi stage: ITER: 71771 RES: 1.1779e-05
ELA: 16.639 Jacobi stage: ITER: 72211 RES: 1.1009e-05
ELA: 16.739 Jacobi stage: ITER: 72651 RES: 1.029e-05
ELA: 16.84 Jacobi stage: ITER: 73091 RES: 9.6175e-06
ELA: 16.94 Jacobi stage: ITER: 73531 RES: 8.9889e-06
ELA: 17.041 Jacobi stage: ITER: 73971 RES: 8.4015e-06
ELA: 17.141 Jacobi stage: ITER: 74411 RES: 7.8525e-06
ELA: 17.241 Jacobi stage: ITER: 74849 RES: 7.3393e-06
ELA: 17.341 Jacobi stage: ITER: 75287 RES: 6.8639e-06
ELA: 17.441 Jacobi stage: ITER: 75727 RES: 6.4153e-06
ELA: 17.542 Jacobi stage: ITER: 76165 RES: 5.9961e-06
ELA: 17.642 Jacobi stage: ITER: 76607 RES: 5.6042e-06
ELA: 17.743 Jacobi stage: ITER: 77047 RES: 5.238e-06
ELA: 17.843 Jacobi stage: ITER: 77483 RES: 4.8987e-06
ELA: 17.943 Jacobi stage: ITER: 77923 RES: 4.5785e-06
ELA: 18.043 Jacobi stage: ITER: 78361 RES: 4.2793e-06
ELA: 18.143 Jacobi stage: ITER: 78799 RES: 4.0021e-06
ELA: 18.243 Jacobi stage: ITER: 79237 RES: 3.7406e-06
ELA: 18.343 Jacobi stage: ITER: 79675 RES: 3.4983e-06
ELA: 18.443 Jacobi stage: ITER: 80115 RES: 3.2696e-06
ELA: 18.543 Jacobi stage: ITER: 80551 RES: 3.0579e-06
ELA: 18.644 Jacobi stage: ITER: 80991 RES: 2.858e-06
ELA: 18.744 Jacobi stage: ITER: 81597 RES: 2.6032e-06
ELA: 18.844 Jacobi stage: ITER: 82275 RES: 2.3464e-06
ELA: 18.944 Jacobi stage: ITER: 82939 RES: 2.1189e-06
ELA: 19.044 Jacobi stage: ITER: 83601 RES: 1.9135e-06
ELA: 19.144 Jacobi stage: ITER: 84261 RES: 1.729e-06
ELA: 19.244 Jacobi stage: ITER: 84919 RES: 1.5633e-06
ELA: 19.344 Jacobi stage: ITER: 85591 RES: 1.41e-06
ELA: 19.444 Jacobi stage: ITER: 86251 RES: 1.274e-06
ELA: 19.545 Jacobi stage: ITER: 86911 RES: 1.1512e-06
ELA: 19.645 Jacobi stage: ITER: 87569 RES: 1.0402e-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