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:57
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: 48723 RES: 0.00040608
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: 1187 RES: 0.7035
Jacobi stage: ITER: 1795 RES: 0.6089
Jacobi stage: ITER: 2411 RES: 0.53476
Jacobi stage: ITER: 3011 RES: 0.47653
Jacobi stage: ITER: 3535 RES: 0.43361
Jacobi stage: ITER: 3943 RES: 0.40407
Jacobi stage: ITER: 4349 RES: 0.37726
Jacobi stage: ITER: 4753 RES: 0.35297
Jacobi stage: ITER: 5151 RES: 0.33103
Jacobi stage: ITER: 5559 RES: 0.31011
Jacobi stage: ITER: 5967 RES: 0.29069
Jacobi stage: ITER: 6371 RES: 0.2728
Jacobi stage: ITER: 6779 RES: 0.25594
Jacobi stage: ITER: 7171 RES: 0.24078
Jacobi stage: ITER: 7571 RES: 0.22629
Jacobi stage: ITER: 7973 RES: 0.21257
Jacobi stage: ITER: 8303 RES: 0.20206
Jacobi stage: ITER: 8485 RES: 0.1964
Jacobi stage: ITER: 8667 RES: 0.19102
Jacobi stage: ITER: 8839 RES: 0.18602
Jacobi stage: ITER: 9009 RES: 0.18116
Jacobi stage: ITER: 9183 RES: 0.17642
Jacobi stage: ITER: 9363 RES: 0.17159
Jacobi stage: ITER: 9541 RES: 0.1669
Jacobi stage: ITER: 9719 RES: 0.16244
Jacobi stage: ITER: 9893 RES: 0.1581
Jacobi stage: ITER: 10063 RES: 0.15407
Jacobi stage: ITER: 10243 RES: 0.14986
Jacobi stage: ITER: 10419 RES: 0.14586
Jacobi stage: ITER: 10779 RES: 0.138
Jacobi stage: ITER: 11191 RES: 0.12953
Jacobi stage: ITER: 11603 RES: 0.12158
Jacobi stage: ITER: 12015 RES: 0.11412
Jacobi stage: ITER: 12427 RES: 0.10712
Jacobi stage: ITER: 12843 RES: 0.10049
Jacobi stage: ITER: 13259 RES: 0.094269
Jacobi stage: ITER: 13681 RES: 0.088325
Jacobi stage: ITER: 14075 RES: 0.083163
Jacobi stage: ITER: 14467 RES: 0.078303
Jacobi stage: ITER: 14879 RES: 0.073501
Jacobi stage: ITER: 15303 RES: 0.068866
Jacobi stage: ITER: 15717 RES: 0.064603
Jacobi stage: ITER: 16139 RES: 0.060567
Jacobi stage: ITER: 16543 RES: 0.056923
Jacobi stage: ITER: 16947 RES: 0.053498
Jacobi stage: ITER: 17347 RES: 0.05031
Jacobi stage: ITER: 17759 RES: 0.047224
Jacobi stage: ITER: 18171 RES: 0.044329
Jacobi stage: ITER: 18575 RES: 0.041661
Jacobi stage: ITER: 18967 RES: 0.039227
Jacobi stage: ITER: 19359 RES: 0.036935
Jacobi stage: ITER: 19755 RES: 0.034755
Jacobi stage: ITER: 20153 RES: 0.032684
Jacobi stage: ITER: 20557 RES: 0.030717
Jacobi stage: ITER: 20967 RES: 0.028851
Jacobi stage: ITER: 21375 RES: 0.027099
Jacobi stage: ITER: 21787 RES: 0.025437
Jacobi stage: ITER: 22195 RES: 0.023892
Jacobi stage: ITER: 22595 RES: 0.022468
Jacobi stage: ITER: 22995 RES: 0.021129
Jacobi stage: ITER: 23407 RES: 0.019833
Jacobi stage: ITER: 23827 RES: 0.018594
Jacobi stage: ITER: 24231 RES: 0.017476
Jacobi stage: ITER: 24631 RES: 0.016434
Jacobi stage: ITER: 25019 RES: 0.015483
Jacobi stage: ITER: 25427 RES: 0.014543
Jacobi stage: ITER: 25827 RES: 0.013676
Jacobi stage: ITER: 26243 RES: 0.01283
Jacobi stage: ITER: 26659 RES: 0.012035
Jacobi stage: ITER: 27065 RES: 0.011304
Jacobi stage: ITER: 27465 RES: 0.010631
Jacobi stage: ITER: 27867 RES: 0.0099972
Jacobi stage: ITER: 28271 RES: 0.0093957
Jacobi stage: ITER: 28667 RES: 0.0088412
Jacobi stage: ITER: 29075 RES: 0.0083041
Jacobi stage: ITER: 29475 RES: 0.0078093
Jacobi stage: ITER: 29863 RES: 0.0073575
Jacobi stage: ITER: 30247 RES: 0.0069361
Jacobi stage: ITER: 30633 RES: 0.0065348
Jacobi stage: ITER: 31019 RES: 0.0061605
Jacobi stage: ITER: 31403 RES: 0.0058076
Jacobi stage: ITER: 31787 RES: 0.005475
Jacobi stage: ITER: 32173 RES: 0.0051582
Jacobi stage: ITER: 32559 RES: 0.0048628
Jacobi stage: ITER: 32943 RES: 0.0045842
Jacobi stage: ITER: 33327 RES: 0.0043217
Jacobi stage: ITER: 33713 RES: 0.0040716
Jacobi stage: ITER: 34099 RES: 0.0038384
Jacobi stage: ITER: 34483 RES: 0.0036186
Jacobi stage: ITER: 34865 RES: 0.0034113
Jacobi stage: ITER: 35251 RES: 0.0032159
Jacobi stage: ITER: 35635 RES: 0.0030317
Jacobi stage: ITER: 36019 RES: 0.0028581
Jacobi stage: ITER: 36405 RES: 0.0026927
Jacobi stage: ITER: 36791 RES: 0.0025385
Jacobi stage: ITER: 37175 RES: 0.0023931
Jacobi stage: ITER: 37559 RES: 0.002256
Jacobi stage: ITER: 37945 RES: 0.0021255
Jacobi stage: ITER: 38329 RES: 0.0020037
Jacobi stage: ITER: 38715 RES: 0.001889
Jacobi stage: ITER: 39099 RES: 0.0017808
Jacobi stage: ITER: 39483 RES: 0.0016788
Jacobi stage: ITER: 39869 RES: 0.0015817
Jacobi stage: ITER: 40255 RES: 0.0014911
Jacobi stage: ITER: 40639 RES: 0.0014057
Jacobi stage: ITER: 41023 RES: 0.0013251
Jacobi stage: ITER: 41409 RES: 0.0012485
Jacobi stage: ITER: 41795 RES: 0.001177
Jacobi stage: ITER: 42179 RES: 0.0011096
Jacobi stage: ITER: 42563 RES: 0.001046
Jacobi stage: ITER: 42949 RES: 0.00098548
Jacobi stage: ITER: 43333 RES: 0.00092904
Jacobi stage: ITER: 43719 RES: 0.00087583
Jacobi stage: ITER: 44103 RES: 0.00082566
Jacobi stage: ITER: 44489 RES: 0.00077789
Jacobi stage: ITER: 44873 RES: 0.00073334
Jacobi stage: ITER: 45259 RES: 0.00069133
Jacobi stage: ITER: 45643 RES: 0.00065173
Jacobi stage: ITER: 46027 RES: 0.00061441
Jacobi stage: ITER: 46417 RES: 0.0005785
Jacobi stage: ITER: 46805 RES: 0.00054503
Jacobi stage: ITER: 47191 RES: 0.00051381
Jacobi stage: ITER: 47581 RES: 0.00048379
Jacobi stage: ITER: 47969 RES: 0.0004558
Jacobi stage: ITER: 48357 RES: 0.00042943
Jacobi stage: ITER: 48743 RES: 0.00040483
Jacobi stage: ITER: 49133 RES: 0.00038118
Jacobi stage: ITER: 49519 RES: 0.00035934
Jacobi stage: ITER: 49907 RES: 0.00033855
Jacobi stage: ITER: 50293 RES: 0.00031897
Jacobi stage: ITER: 50679 RES: 0.0003007
Jacobi stage: ITER: 51063 RES: 0.00028347
Jacobi stage: ITER: 51451 RES: 0.00026707
Jacobi stage: ITER: 51837 RES: 0.00025162
Jacobi stage: ITER: 52223 RES: 0.00023721
Jacobi stage: ITER: 52609 RES: 0.00022348
Jacobi stage: ITER: 52999 RES: 0.00021055
Jacobi stage: ITER: 53387 RES: 0.00019837
Jacobi stage: ITER: 53771 RES: 0.00018701
Jacobi stage: ITER: 54155 RES: 0.0001763
Jacobi stage: ITER: 54541 RES: 0.0001661
Jacobi stage: ITER: 54927 RES: 0.00015659
Jacobi stage: ITER: 55311 RES: 0.00014762
Jacobi stage: ITER: 55697 RES: 0.00013908
Jacobi stage: ITER: 56083 RES: 0.00013111
Jacobi stage: ITER: 56467 RES: 0.0001236
Jacobi stage: ITER: 56883 RES: 0.00011595
Jacobi stage: ITER: 57267 RES: 0.00010931
Jacobi stage: ITER: 57653 RES: 0.00010298
Jacobi stage: ITER: 58039 RES: 9.7086e-05
Jacobi stage: ITER: 58423 RES: 9.1525e-05
Jacobi stage: ITER: 58809 RES: 8.6229e-05
Jacobi stage: ITER: 59193 RES: 8.129e-05
Jacobi stage: ITER: 59579 RES: 7.6634e-05
Jacobi stage: ITER: 59963 RES: 7.2245e-05
Jacobi stage: ITER: 60349 RES: 6.8065e-05
Jacobi stage: ITER: 60735 RES: 6.4167e-05
Jacobi stage: ITER: 61119 RES: 6.0491e-05
Jacobi stage: ITER: 61503 RES: 5.7026e-05
Jacobi stage: ITER: 61889 RES: 5.3727e-05
Jacobi stage: ITER: 62275 RES: 5.065e-05
Jacobi stage: ITER: 62659 RES: 4.7749e-05
Jacobi stage: ITER: 63045 RES: 4.4986e-05
Jacobi stage: ITER: 63431 RES: 4.2409e-05
Jacobi stage: ITER: 63815 RES: 3.998e-05
Jacobi stage: ITER: 64199 RES: 3.769e-05
Jacobi stage: ITER: 64585 RES: 3.551e-05
Jacobi stage: ITER: 64971 RES: 3.3476e-05
Jacobi stage: ITER: 65355 RES: 3.1558e-05
Jacobi stage: ITER: 65741 RES: 2.9733e-05
Jacobi stage: ITER: 66127 RES: 2.803e-05
Jacobi stage: ITER: 66511 RES: 2.6424e-05
Jacobi stage: ITER: 66895 RES: 2.4911e-05
Jacobi stage: ITER: 67281 RES: 2.3469e-05
Jacobi stage: ITER: 67667 RES: 2.2125e-05
Jacobi stage: ITER: 68051 RES: 2.0858e-05
Jacobi stage: ITER: 68437 RES: 1.9651e-05
Jacobi stage: ITER: 68821 RES: 1.8526e-05
Jacobi stage: ITER: 69207 RES: 1.7464e-05
Jacobi stage: ITER: 69591 RES: 1.6464e-05
Jacobi stage: ITER: 69977 RES: 1.5512e-05
Jacobi stage: ITER: 70363 RES: 1.4623e-05
Jacobi stage: ITER: 70747 RES: 1.3786e-05
Jacobi stage: ITER: 71131 RES: 1.2996e-05
Jacobi stage: ITER: 71517 RES: 1.2244e-05
Jacobi stage: ITER: 71903 RES: 1.1543e-05
Jacobi stage: ITER: 72287 RES: 1.0882e-05
Jacobi stage: ITER: 72673 RES: 1.0252e-05
Jacobi stage: ITER: 73059 RES: 9.6648e-06
Jacobi stage: ITER: 73443 RES: 9.1113e-06
Jacobi stage: ITER: 73827 RES: 8.5894e-06
Jacobi stage: ITER: 74213 RES: 8.0925e-06
Jacobi stage: ITER: 74599 RES: 7.6289e-06
Jacobi stage: ITER: 74983 RES: 7.192e-06
Jacobi stage: ITER: 75367 RES: 6.78e-06
Jacobi stage: ITER: 75753 RES: 6.3878e-06
Jacobi stage: ITER: 76139 RES: 6.0219e-06
Jacobi stage: ITER: 76523 RES: 5.677e-06
Jacobi stage: ITER: 76909 RES: 5.3485e-06
Jacobi stage: ITER: 77293 RES: 5.0422e-06
Jacobi stage: ITER: 77679 RES: 4.7534e-06
Jacobi stage: ITER: 78063 RES: 4.4811e-06
Jacobi stage: ITER: 78449 RES: 4.2219e-06
Jacobi stage: ITER: 78835 RES: 3.98e-06
Jacobi stage: ITER: 79219 RES: 3.7521e-06
Jacobi stage: ITER: 79603 RES: 3.5372e-06
Jacobi stage: ITER: 79989 RES: 3.3325e-06
Jacobi stage: ITER: 80375 RES: 3.1416e-06
Jacobi stage: ITER: 80759 RES: 2.9617e-06
Jacobi stage: ITER: 81145 RES: 2.7903e-06
Jacobi stage: ITER: 81529 RES: 2.6305e-06
Jacobi stage: ITER: 81915 RES: 2.4799e-06
Jacobi stage: ITER: 82299 RES: 2.3378e-06
Jacobi stage: ITER: 82685 RES: 2.2026e-06
Jacobi stage: ITER: 83071 RES: 2.0764e-06
Jacobi stage: ITER: 83455 RES: 1.9575e-06
Jacobi stage: ITER: 83839 RES: 1.8454e-06
Jacobi stage: ITER: 84225 RES: 1.7386e-06
Jacobi stage: ITER: 84611 RES: 1.639e-06
Jacobi stage: ITER: 84995 RES: 1.5451e-06
Jacobi stage: ITER: 85379 RES: 1.4566e-06
Jacobi stage: ITER: 85765 RES: 1.3724e-06
Jacobi stage: ITER: 86151 RES: 1.2937e-06
Jacobi stage: ITER: 86535 RES: 1.2196e-06
Jacobi stage: ITER: 86921 RES: 1.1491e-06
Jacobi stage: ITER: 87305 RES: 1.0833e-06
Jacobi stage: ITER: 87691 RES: 1.0212e-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:2.1019e-
ELA: 0.10014 Jacobi stage: ITER: 51235 RES: 0.00027591
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.3746e-
ELA: 0.10012 Jacobi stage: ITER: 585 RES: 0.82916
ELA: 0.20025 Jacobi stage: ITER: 1199 RES: 0.70136
ELA: 0.30038 Jacobi stage: ITER: 1813 RES: 0.6062
ELA: 0.4005 Jacobi stage: ITER: 2415 RES: 0.53433
ELA: 0.50062 Jacobi stage: ITER: 2921 RES: 0.48439
ELA: 0.60077 Jacobi stage: ITER: 3331 RES: 0.44959
ELA: 0.7009 Jacobi stage: ITER: 3737 RES: 0.41847
ELA: 0.80126 Jacobi stage: ITER: 4143 RES: 0.39062
ELA: 0.90139 Jacobi stage: ITER: 4547 RES: 0.36522
ELA: 1.0015 Jacobi stage: ITER: 4955 RES: 0.34167
ELA: 1.1016 Jacobi stage: ITER: 5361 RES: 0.31996
ELA: 1.2017 Jacobi stage: ITER: 5771 RES: 0.29984
ELA: 1.3019 Jacobi stage: ITER: 6179 RES: 0.28114
ELA: 1.402 Jacobi stage: ITER: 6571 RES: 0.26439
ELA: 1.5021 Jacobi stage: ITER: 6971 RES: 0.24839
ELA: 1.6022 Jacobi stage: ITER: 7373 RES: 0.23327
ELA: 1.7022 Jacobi stage: ITER: 7687 RES: 0.22226
ELA: 1.8024 Jacobi stage: ITER: 7867 RES: 0.21615
ELA: 1.9025 Jacobi stage: ITER: 8047 RES: 0.21022
ELA: 2.0026 Jacobi stage: ITER: 8221 RES: 0.20457
ELA: 2.1028 Jacobi stage: ITER: 8395 RES: 0.19921
ELA: 2.2029 Jacobi stage: ITER: 8573 RES: 0.19375
ELA: 2.303 Jacobi stage: ITER: 8755 RES: 0.18845
ELA: 2.4031 Jacobi stage: ITER: 8939 RES: 0.18318
ELA: 2.5033 Jacobi stage: ITER: 9115 RES: 0.17817
ELA: 2.6043 Jacobi stage: ITER: 9295 RES: 0.1734
ELA: 2.7053 Jacobi stage: ITER: 9469 RES: 0.16876
ELA: 2.8054 Jacobi stage: ITER: 9651 RES: 0.16415
ELA: 2.9055 Jacobi stage: ITER: 9831 RES: 0.15967
ELA: 3.0056 Jacobi stage: ITER: 10217 RES: 0.15041
ELA: 3.1057 Jacobi stage: ITER: 10629 RES: 0.14118
ELA: 3.2058 Jacobi stage: ITER: 11043 RES: 0.13251
ELA: 3.3059 Jacobi stage: ITER: 11453 RES: 0.12438
ELA: 3.406 Jacobi stage: ITER: 11859 RES: 0.11689
ELA: 3.5061 Jacobi stage: ITER: 12275 RES: 0.10965
ELA: 3.6062 Jacobi stage: ITER: 12691 RES: 0.10287
ELA: 3.7064 Jacobi stage: ITER: 13107 RES: 0.096497
ELA: 3.8065 Jacobi stage: ITER: 13501 RES: 0.090801
ELA: 3.9066 Jacobi stage: ITER: 13897 RES: 0.085442
ELA: 4.0067 Jacobi stage: ITER: 14307 RES: 0.080251
ELA: 4.1069 Jacobi stage: ITER: 14731 RES: 0.075191
ELA: 4.207 Jacobi stage: ITER: 15147 RES: 0.070536
ELA: 4.3071 Jacobi stage: ITER: 15567 RES: 0.066129
ELA: 4.4073 Jacobi stage: ITER: 15975 RES: 0.062112
ELA: 4.5074 Jacobi stage: ITER: 16377 RES: 0.058375
ELA: 4.6081 Jacobi stage: ITER: 16783 RES: 0.054862
ELA: 4.7082 Jacobi stage: ITER: 17195 RES: 0.051498
ELA: 4.8084 Jacobi stage: ITER: 17609 RES: 0.04831
ELA: 4.9085 Jacobi stage: ITER: 18011 RES: 0.045431
ELA: 5.0086 Jacobi stage: ITER: 18399 RES: 0.042803
ELA: 5.1087 Jacobi stage: ITER: 18803 RES: 0.040228
ELA: 5.2089 Jacobi stage: ITER: 19191 RES: 0.0379
ELA: 5.309 Jacobi stage: ITER: 19597 RES: 0.035598
ELA: 5.4091 Jacobi stage: ITER: 20007 RES: 0.033435
ELA: 5.5093 Jacobi stage: ITER: 20415 RES: 0.031404
ELA: 5.6094 Jacobi stage: ITER: 20827 RES: 0.029478
ELA: 5.7095 Jacobi stage: ITER: 21237 RES: 0.027671
ELA: 5.8103 Jacobi stage: ITER: 21645 RES: 0.02599
ELA: 5.9104 Jacobi stage: ITER: 22051 RES: 0.024426
ELA: 6.0105 Jacobi stage: ITER: 22451 RES: 0.022971
ELA: 6.1106 Jacobi stage: ITER: 22865 RES: 0.021549
ELA: 6.2107 Jacobi stage: ITER: 23285 RES: 0.020202
ELA: 6.3108 Jacobi stage: ITER: 23691 RES: 0.018987
ELA: 6.4109 Jacobi stage: ITER: 24087 RES: 0.017866
ELA: 6.5111 Jacobi stage: ITER: 24483 RES: 0.016812
ELA: 6.6112 Jacobi stage: ITER: 24891 RES: 0.015791
ELA: 6.7113 Jacobi stage: ITER: 25299 RES: 0.014832
ELA: 6.8114 Jacobi stage: ITER: 25713 RES: 0.013913
ELA: 6.9115 Jacobi stage: ITER: 26127 RES: 0.01306
ELA: 7.0117 Jacobi stage: ITER: 26535 RES: 0.012267
ELA: 7.1118 Jacobi stage: ITER: 26939 RES: 0.011529
ELA: 7.2119 Jacobi stage: ITER: 27335 RES: 0.010848
ELA: 7.312 Jacobi stage: ITER: 27743 RES: 0.010189
ELA: 7.4121 Jacobi stage: ITER: 28139 RES: 0.0095881
ELA: 7.5123 Jacobi stage: ITER: 28555 RES: 0.0089946
ELA: 7.6124 Jacobi stage: ITER: 28957 RES: 0.0084534
ELA: 7.7125 Jacobi stage: ITER: 29343 RES: 0.0079692
ELA: 7.8126 Jacobi stage: ITER: 29729 RES: 0.0075082
ELA: 7.9128 Jacobi stage: ITER: 30113 RES: 0.0070781
ELA: 8.0129 Jacobi stage: ITER: 30499 RES: 0.0066727
ELA: 8.113 Jacobi stage: ITER: 30883 RES: 0.0062905
ELA: 8.2131 Jacobi stage: ITER: 31269 RES: 0.0059266
ELA: 8.3133 Jacobi stage: ITER: 31655 RES: 0.0055871
ELA: 8.4134 Jacobi stage: ITER: 32039 RES: 0.0052671
ELA: 8.5135 Jacobi stage: ITER: 32425 RES: 0.0049624
ELA: 8.6136 Jacobi stage: ITER: 32811 RES: 0.0046781
ELA: 8.7138 Jacobi stage: ITER: 33195 RES: 0.0044102
ELA: 8.8139 Jacobi stage: ITER: 33579 RES: 0.0041576
ELA: 8.9143 Jacobi stage: ITER: 33967 RES: 0.003917
ELA: 9.0144 Jacobi stage: ITER: 34337 RES: 0.0036995
ELA: 9.1145 Jacobi stage: ITER: 34723 RES: 0.0034876
ELA: 9.2146 Jacobi stage: ITER: 35107 RES: 0.0032878
ELA: 9.3147 Jacobi stage: ITER: 35493 RES: 0.0030976
ELA: 9.4148 Jacobi stage: ITER: 35877 RES: 0.0029202
ELA: 9.5149 Jacobi stage: ITER: 36263 RES: 0.0027529
ELA: 9.615 Jacobi stage: ITER: 36647 RES: 0.0025953
ELA: 9.7151 Jacobi stage: ITER: 37031 RES: 0.0024466
ELA: 9.8152 Jacobi stage: ITER: 37417 RES: 0.0023051
ELA: 9.9153 Jacobi stage: ITER: 37803 RES: 0.002173
ELA: 10.015 Jacobi stage: ITER: 38187 RES: 0.0020486
ELA: 10.115 Jacobi stage: ITER: 38573 RES: 0.00193
ELA: 10.216 Jacobi stage: ITER: 38957 RES: 0.0018195
ELA: 10.316 Jacobi stage: ITER: 39343 RES: 0.0017153
ELA: 10.416 Jacobi stage: ITER: 39727 RES: 0.001617
ELA: 10.516 Jacobi stage: ITER: 40113 RES: 0.0015235
ELA: 10.616 Jacobi stage: ITER: 40499 RES: 0.0014362
ELA: 10.716 Jacobi stage: ITER: 40883 RES: 0.001354
ELA: 10.816 Jacobi stage: ITER: 41267 RES: 0.0012764
ELA: 10.916 Jacobi stage: ITER: 41653 RES: 0.0012026
ELA: 11.016 Jacobi stage: ITER: 42039 RES: 0.0011337
ELA: 11.116 Jacobi stage: ITER: 42423 RES: 0.0010687
ELA: 11.217 Jacobi stage: ITER: 42807 RES: 0.0010075
ELA: 11.317 Jacobi stage: ITER: 43193 RES: 0.00094923
ELA: 11.417 Jacobi stage: ITER: 43579 RES: 0.00089486
ELA: 11.517 Jacobi stage: ITER: 43963 RES: 0.00084361
ELA: 11.617 Jacobi stage: ITER: 44347 RES: 0.00079529
ELA: 11.717 Jacobi stage: ITER: 44733 RES: 0.00074928
ELA: 11.817 Jacobi stage: ITER: 45119 RES: 0.00070636
ELA: 11.917 Jacobi stage: ITER: 45503 RES: 0.0006659
ELA: 12.017 Jacobi stage: ITER: 45893 RES: 0.00062699
ELA: 12.118 Jacobi stage: ITER: 46281 RES: 0.00059071
ELA: 12.218 Jacobi stage: ITER: 46667 RES: 0.00055688
ELA: 12.318 Jacobi stage: ITER: 47057 RES: 0.00052434
ELA: 12.418 Jacobi stage: ITER: 47443 RES: 0.00049431
ELA: 12.518 Jacobi stage: ITER: 47833 RES: 0.00046542
ELA: 12.618 Jacobi stage: ITER: 48217 RES: 0.00043876
ELA: 12.718 Jacobi stage: ITER: 48609 RES: 0.00041312
ELA: 12.818 Jacobi stage: ITER: 48993 RES: 0.00038946
ELA: 12.919 Jacobi stage: ITER: 49383 RES: 0.00036693
ELA: 13.019 Jacobi stage: ITER: 49769 RES: 0.0003457
ELA: 13.119 Jacobi stage: ITER: 50155 RES: 0.0003259
ELA: 13.219 Jacobi stage: ITER: 50539 RES: 0.00030723
ELA: 13.319 Jacobi stage: ITER: 50927 RES: 0.00028946
ELA: 13.419 Jacobi stage: ITER: 51315 RES: 0.00027271
ELA: 13.519 Jacobi stage: ITER: 51699 RES: 0.00025709
ELA: 13.619 Jacobi stage: ITER: 52087 RES: 0.00024222
ELA: 13.719 Jacobi stage: ITER: 52477 RES: 0.00022806
ELA: 13.82 Jacobi stage: ITER: 52863 RES: 0.000215
ELA: 13.92 Jacobi stage: ITER: 53247 RES: 0.00020268
ELA: 14.02 Jacobi stage: ITER: 53633 RES: 0.00019096
ELA: 14.12 Jacobi stage: ITER: 54017 RES: 0.00018002
ELA: 14.22 Jacobi stage: ITER: 54401 RES: 0.00016971
ELA: 14.32 Jacobi stage: ITER: 54787 RES: 0.00015999
ELA: 14.42 Jacobi stage: ITER: 55171 RES: 0.00015083
ELA: 14.52 Jacobi stage: ITER: 55557 RES: 0.0001421
ELA: 14.62 Jacobi stage: ITER: 55943 RES: 0.00013396
ELA: 14.72 Jacobi stage: ITER: 56327 RES: 0.00012629
ELA: 14.82 Jacobi stage: ITER: 56711 RES: 0.00011905
ELA: 14.92 Jacobi stage: ITER: 57097 RES: 0.00011217
ELA: 15.02 Jacobi stage: ITER: 57481 RES: 0.00010574
ELA: 15.121 Jacobi stage: ITER: 57867 RES: 9.9685e-05
ELA: 15.221 Jacobi stage: ITER: 58251 RES: 9.3975e-05
ELA: 15.321 Jacobi stage: ITER: 58635 RES: 8.8592e-05
ELA: 15.421 Jacobi stage: ITER: 59021 RES: 8.3467e-05
ELA: 15.521 Jacobi stage: ITER: 59407 RES: 7.8686e-05
ELA: 15.621 Jacobi stage: ITER: 59791 RES: 7.4179e-05
ELA: 15.721 Jacobi stage: ITER: 60175 RES: 6.993e-05
ELA: 15.821 Jacobi stage: ITER: 60561 RES: 6.5884e-05
ELA: 15.921 Jacobi stage: ITER: 60947 RES: 6.2111e-05
ELA: 16.021 Jacobi stage: ITER: 61331 RES: 5.8553e-05
ELA: 16.121 Jacobi stage: ITER: 61715 RES: 5.5199e-05
ELA: 16.221 Jacobi stage: ITER: 62101 RES: 5.2006e-05
ELA: 16.321 Jacobi stage: ITER: 62485 RES: 4.9027e-05
ELA: 16.421 Jacobi stage: ITER: 62871 RES: 4.6219e-05
ELA: 16.521 Jacobi stage: ITER: 63255 RES: 4.3572e-05
ELA: 16.621 Jacobi stage: ITER: 63639 RES: 4.1076e-05
ELA: 16.721 Jacobi stage: ITER: 64025 RES: 3.8699e-05
ELA: 16.821 Jacobi stage: ITER: 64411 RES: 3.6483e-05
ELA: 16.922 Jacobi stage: ITER: 64795 RES: 3.4393e-05
ELA: 17.022 Jacobi stage: ITER: 65181 RES: 3.2403e-05
ELA: 17.122 Jacobi stage: ITER: 65565 RES: 3.0547e-05
ELA: 17.222 Jacobi stage: ITER: 65951 RES: 2.8798e-05
ELA: 17.322 Jacobi stage: ITER: 66335 RES: 2.7148e-05
ELA: 17.422 Jacobi stage: ITER: 66719 RES: 2.5593e-05
ELA: 17.522 Jacobi stage: ITER: 67105 RES: 2.4113e-05
ELA: 17.622 Jacobi stage: ITER: 67491 RES: 2.2731e-05
ELA: 17.722 Jacobi stage: ITER: 67875 RES: 2.1429e-05
ELA: 17.822 Jacobi stage: ITER: 68259 RES: 2.0202e-05
ELA: 17.922 Jacobi stage: ITER: 68645 RES: 1.9033e-05
ELA: 18.022 Jacobi stage: ITER: 69031 RES: 1.7943e-05
ELA: 18.122 Jacobi stage: ITER: 69415 RES: 1.6915e-05
ELA: 18.222 Jacobi stage: ITER: 69801 RES: 1.5937e-05
ELA: 18.322 Jacobi stage: ITER: 70185 RES: 1.5024e-05
ELA: 18.423 Jacobi stage: ITER: 70571 RES: 1.4163e-05
ELA: 18.523 Jacobi stage: ITER: 70955 RES: 1.3352e-05
ELA: 18.623 Jacobi stage: ITER: 71341 RES: 1.258e-05
ELA: 18.723 Jacobi stage: ITER: 71725 RES: 1.1859e-05
ELA: 18.823 Jacobi stage: ITER: 72111 RES: 1.118e-05
ELA: 18.923 Jacobi stage: ITER: 72495 RES: 1.0539e-05
ELA: 19.023 Jacobi stage: ITER: 72881 RES: 9.9297e-06
ELA: 19.123 Jacobi stage: ITER: 73267 RES: 9.3609e-06
ELA: 19.223 Jacobi stage: ITER: 73651 RES: 8.8248e-06
ELA: 19.323 Jacobi stage: ITER: 74035 RES: 8.3193e-06
ELA: 19.424 Jacobi stage: ITER: 74421 RES: 7.838e-06
ELA: 19.524 Jacobi stage: ITER: 74807 RES: 7.3891e-06
ELA: 19.624 Jacobi stage: ITER: 75191 RES: 6.9658e-06
ELA: 19.724 Jacobi stage: ITER: 75577 RES: 6.5628e-06
ELA: 19.824 Jacobi stage: ITER: 75961 RES: 6.1869e-06
ELA: 19.924 Jacobi stage: ITER: 76347 RES: 5.8325e-06
ELA: 20.024 Jacobi stage: ITER: 76731 RES: 5.4985e-06
ELA: 20.124 Jacobi stage: ITER: 77117 RES: 5.1804e-06
ELA: 20.224 Jacobi stage: ITER: 77503 RES: 4.8836e-06
ELA: 20.324 Jacobi stage: ITER: 77887 RES: 4.6039e-06
ELA: 20.425 Jacobi stage: ITER: 78273 RES: 4.3375e-06
ELA: 20.525 Jacobi stage: ITER: 78659 RES: 4.0891e-06
ELA: 20.625 Jacobi stage: ITER: 79043 RES: 3.8549e-06
ELA: 20.725 Jacobi stage: ITER: 79427 RES: 3.6341e-06
ELA: 20.825 Jacobi stage: ITER: 79813 RES: 3.4238e-06
ELA: 20.925 Jacobi stage: ITER: 80199 RES: 3.2277e-06
ELA: 21.025 Jacobi stage: ITER: 80583 RES: 3.0429e-06
ELA: 21.125 Jacobi stage: ITER: 80969 RES: 2.8668e-06
ELA: 21.225 Jacobi stage: ITER: 81353 RES: 2.7026e-06
ELA: 21.325 Jacobi stage: ITER: 81739 RES: 2.5478e-06
ELA: 21.426 Jacobi stage: ITER: 82123 RES: 2.4019e-06
ELA: 21.526 Jacobi stage: ITER: 82509 RES: 2.2629e-06
ELA: 21.626 Jacobi stage: ITER: 82895 RES: 2.1333e-06
ELA: 21.726 Jacobi stage: ITER: 83279 RES: 2.0111e-06
ELA: 21.826 Jacobi stage: ITER: 83663 RES: 1.8959e-06
ELA: 21.926 Jacobi stage: ITER: 84049 RES: 1.7862e-06
ELA: 22.026 Jacobi stage: ITER: 84435 RES: 1.6839e-06
ELA: 22.126 Jacobi stage: ITER: 84819 RES: 1.5875e-06
ELA: 22.226 Jacobi stage: ITER: 85205 RES: 1.4956e-06
ELA: 22.327 Jacobi stage: ITER: 85591 RES: 1.41e-06
ELA: 22.427 Jacobi stage: ITER: 85975 RES: 1.3292e-06
ELA: 22.527 Jacobi stage: ITER: 86359 RES: 1.2531e-06
ELA: 22.627 Jacobi stage: ITER: 86745 RES: 1.1806e-06
ELA: 22.727 Jacobi stage: ITER: 87131 RES: 1.1129e-06
ELA: 22.827 Jacobi stage: ITER: 87637 RES: 1.0294e-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