Template Numerical Library version\ main:52827a2
Loading...
Searching...
No Matches
Linear solvers

Table of Contents

Introduction

Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the following iterative methods:

  1. Stationary methods
    1. Jacobi method (TNL::Solvers::Linear::Jacobi)
    2. Successive-overrelaxation method, SOR (TNL::Solvers::Linear::SOR)
  2. Krylov subspace methods
    1. Conjugate gradient method, CG (TNL::Solvers::Linear::CG)
    2. Biconjugate gradient stabilized method, BICGStab (TNL::Solvers::Linear::BICGStab)
    3. Biconjugate gradient stabilized method, BICGStab(l) (TNL::Solvers::Linear::BICGStabL)
    4. Transpose-free quasi-minimal residual method, TFQMR (TNL::Solvers::Linear::TFQMR)
    5. Generalized minimal residual method, GMRES (TNL::Solvers::Linear::GMRES) with various methods of orthogonalization:
      1. Classical Gramm-Schmidt, CGS
      2. Classical Gramm-Schmidt with reorthogonalization, CGSR
      3. Modified Gramm-Schmidt, MGS
      4. Modified Gramm-Schmidt with reorthogonalization, MGSR
      5. Compact WY form of the Householder reflections, CWY

The iterative solvers (not the stationary solvers like TNL::Solvers::Linear::Jacobi and TNL::Solvers::Linear::SOR) can be combined with the following preconditioners:

  1. Diagonal or Jacobi (TNL::Solvers::Linear::Preconditioners::Diagonal)
  2. ILU (Incomplete LU) - CPU only currently
    1. ILU(0) (TNL::Solvers::Linear::Preconditioners::ILU0)
    2. ILUT (ILU with thresholding) (TNL::Solvers::Linear::Preconditioners::ILUT)

Iterative solvers of linear systems

Basic setup

All iterative solvers for linear systems can be found in the namespace TNL::Solvers::Linear. The following example shows the use the iterative solvers:

1#include <iostream>
2#include <memory>
3#include <TNL/Matrices/SparseMatrix.h>
4#include <TNL/Devices/Host.h>
5#include <TNL/Devices/Cuda.h>
6#include <TNL/Solvers/Linear/TFQMR.h>
7
8template< typename Device >
9void
10iterativeLinearSolverExample()
11{
12 /***
13 * Set the following matrix (dots represent zero matrix elements):
14 *
15 * / 2.5 -1 . . . \
16 * | -1 2.5 -1 . . |
17 * | . -1 2.5 -1. . |
18 * | . . -1 2.5 -1 |
19 * \ . . . -1 2.5 /
20 */
21 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
22 using Vector = TNL::Containers::Vector< double, Device >;
23 const int size( 5 );
24 auto matrix_ptr = std::make_shared< MatrixType >();
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
28 auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
29 {
30 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 ) {
32 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
33 row.setElement( 1, rowIdx + 1, -1 ); // element above the diagonal
34 }
35 else if( rowIdx == size - 1 ) {
36 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
37 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
38 }
39 else {
40 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
41 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
42 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
43 }
44 };
45
46 /***
47 * Set the matrix elements.
48 */
49 matrix_ptr->forAllRows( f );
50 std::cout << *matrix_ptr << std::endl;
51
52 /***
53 * Set the right-hand side vector.
54 */
55 Vector x( size, 1.0 );
56 Vector b( size );
57 matrix_ptr->vectorProduct( x, b );
58 x = 0.0;
59 std::cout << "Vector b = " << b << std::endl;
60
61 /***
62 * Solve the linear system.
63 */
64 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
65 LinearSolver solver;
66 solver.setMatrix( matrix_ptr );
67 solver.setConvergenceResidue( 1.0e-6 );
68 solver.solve( b, x );
69 std::cout << "Vector x = " << x << std::endl;
70}
71
72int
73main( int argc, char* argv[] )
74{
75 std::cout << "Solving linear system on host: " << std::endl;
76 iterativeLinearSolverExample< TNL::Devices::Sequential >();
77
78#ifdef __CUDACC__
79 std::cout << "Solving linear system on CUDA device: " << std::endl;
80 iterativeLinearSolverExample< TNL::Devices::Cuda >();
81#endif
82}
__cuda_callable__
#define __cuda_callable__
Definition Macros.h:49
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition Vector.h:36
TNL::Matrices::SparseMatrix
Implementation of sparse matrix, i.e. matrix storing only non-zero elements.
Definition SparseMatrix.h:57
TNL::Solvers::Linear::LinearSolver::setMatrix
void setMatrix(const MatrixPointer &matrix)
Set the matrix of the linear system.
Definition LinearSolver.h:120
TNL::Solvers::Linear::TFQMR
Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method.
Definition TFQMR.h:21
std::endl
T endl(T... args)
std::make_shared
T make_shared(T... args)

In this example we solve a linear system \( A \vec x = \vec b \) where

\[ A = \left( \begin{array}{cccc} 2.5 & -1 & & & \\ -1 & 2.5 & -1 & & \\ & -1 & 2.5 & -1 & \\ & & -1 & 2.5 & -1 \\ & & & -1 & 2.5 \\ \end{array} \right) \]

The right-hand side vector \(\vec b \) is set to \(( 1.5, 0.5, 0.5, 0.5, 1.5 )^T \) so that the exact solution is \( \vec x = ( 1, 1, 1, 1, 1 )^T \). The 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 <chrono>
4#include <thread>
5#include <TNL/Matrices/SparseMatrix.h>
6#include <TNL/Devices/Sequential.h>
7#include <TNL/Devices/Cuda.h>
8#include <TNL/Solvers/Linear/Jacobi.h>
9
10template< typename Device >
11void
12iterativeLinearSolverExample()
13{
14 /***
15 * Set the following matrix (dots represent zero matrix elements):
16 *
17 * / 2.5 -1 . . . \
18 * | -1 2.5 -1 . . |
19 * | . -1 2.5 -1. . |
20 * | . . -1 2.5 -1 |
21 * \ . . . -1 2.5 /
22 */
23 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
24 using Vector = TNL::Containers::Vector< double, Device >;
25 const int size( 5 );
26 auto matrix_ptr = std::make_shared< MatrixType >();
27 matrix_ptr->setDimensions( size, size );
28 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
29
30 auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
31 {
32 const int rowIdx = row.getRowIndex();
33 if( rowIdx == 0 ) {
34 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
35 row.setElement( 1, rowIdx + 1, -1 ); // element above the diagonal
36 }
37 else if( rowIdx == size - 1 ) {
38 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
39 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
40 }
41 else {
42 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
43 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
44 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
45 }
46 };
47
48 /***
49 * Set the matrix elements.
50 */
51 matrix_ptr->forAllRows( f );
52 std::cout << *matrix_ptr << std::endl;
53
54 /***
55 * Set the right-hand side vector.
56 */
57 Vector x( size, 1.0 );
58 Vector b( size );
59 matrix_ptr->vectorProduct( x, b );
60 x = 0.0;
61 std::cout << "Vector b = " << b << std::endl;
62
63 /***
64 * Setup solver of the linear system.
65 */
66 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
67 LinearSolver solver;
68 solver.setMatrix( matrix_ptr );
69 solver.setOmega( 0.0005 );
70
71 /***
72 * Setup monitor of the iterative solver.
73 */
74 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
75 IterativeSolverMonitorType monitor;
76 TNL::Solvers::SolverMonitorThread monitorThread( monitor );
77 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
78 monitor.setVerbose( 1 );
79 monitor.setStage( "Jacobi stage:" );
80 solver.setSolverMonitor( monitor );
81 solver.setConvergenceResidue( 1.0e-6 );
82 solver.solve( b, x );
83 monitor.stopMainLoop();
84 std::cout << "Vector x = " << x << std::endl;
85}
86
87int
88main( int argc, char* argv[] )
89{
90 std::cout << "Solving linear system on host: " << std::endl;
91 iterativeLinearSolverExample< TNL::Devices::Sequential >();
92
93#ifdef __CUDACC__
94 std::cout << "Solving linear system on CUDA device: " << std::endl;
95 iterativeLinearSolverExample< TNL::Devices::Cuda >();
96#endif
97}
TNL::Solvers::Linear::Jacobi
Iterative solver of linear systems based on the Jacobi method.
Definition Jacobi.h:21
TNL::Solvers::SolverMonitorThread
A RAII wrapper for launching the SolverMonitor's main loop in a separate thread.
Definition SolverMonitor.h:133

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:

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: 40067 RES: 0.0015338
Jacobi stage: ITER: 81289 RES: 2.7293e-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: 4255 RES: 0.38335
Jacobi stage: ITER: 8953 RES: 0.18273
Jacobi stage: ITER: 10329 RES: 0.14784
Jacobi stage: ITER: 11123 RES: 0.13089
Jacobi stage: ITER: 11707 RES: 0.11966
Jacobi stage: ITER: 12317 RES: 0.10892
Jacobi stage: ITER: 12939 RES: 0.09902
Jacobi stage: ITER: 13555 RES: 0.090079
Jacobi stage: ITER: 14189 RES: 0.081694
Jacobi stage: ITER: 14813 RES: 0.074227
Jacobi stage: ITER: 15429 RES: 0.067525
Jacobi stage: ITER: 16055 RES: 0.061354
Jacobi stage: ITER: 16689 RES: 0.055643
Jacobi stage: ITER: 17303 RES: 0.050651
Jacobi stage: ITER: 17933 RES: 0.045965
Jacobi stage: ITER: 18565 RES: 0.041713
Jacobi stage: ITER: 19165 RES: 0.03804
Jacobi stage: ITER: 19809 RES: 0.034457
Jacobi stage: ITER: 20443 RES: 0.03127
Jacobi stage: ITER: 21075 RES: 0.028377
Jacobi stage: ITER: 21669 RES: 0.025894
Jacobi stage: ITER: 22273 RES: 0.0236
Jacobi stage: ITER: 22923 RES: 0.021364
Jacobi stage: ITER: 23575 RES: 0.019328
Jacobi stage: ITER: 24227 RES: 0.017486
Jacobi stage: ITER: 24879 RES: 0.01582
Jacobi stage: ITER: 25529 RES: 0.014312
Jacobi stage: ITER: 26179 RES: 0.012948
Jacobi stage: ITER: 26831 RES: 0.011722
Jacobi stage: ITER: 27483 RES: 0.010605
Jacobi stage: ITER: 28133 RES: 0.009594
Jacobi stage: ITER: 28785 RES: 0.0086797
Jacobi stage: ITER: 29435 RES: 0.0078574
Jacobi stage: ITER: 30087 RES: 0.0071086
Jacobi stage: ITER: 30739 RES: 0.0064312
Jacobi stage: ITER: 31391 RES: 0.0058183
Jacobi stage: ITER: 32041 RES: 0.0052639
Jacobi stage: ITER: 32691 RES: 0.0047652
Jacobi stage: ITER: 33337 RES: 0.0043137
Jacobi stage: ITER: 33983 RES: 0.0039074
Jacobi stage: ITER: 34629 RES: 0.0035372
Jacobi stage: ITER: 35275 RES: 0.0032041
Jacobi stage: ITER: 35921 RES: 0.0029005
Jacobi stage: ITER: 36567 RES: 0.0026273
Jacobi stage: ITER: 37215 RES: 0.0023784
Jacobi stage: ITER: 37861 RES: 0.0021531
Jacobi stage: ITER: 38507 RES: 0.0019503
Jacobi stage: ITER: 39155 RES: 0.0017655
Jacobi stage: ITER: 39801 RES: 0.0015983
Jacobi stage: ITER: 40447 RES: 0.0014477
Jacobi stage: ITER: 41095 RES: 0.0013106
Jacobi stage: ITER: 41741 RES: 0.0011864
Jacobi stage: ITER: 42387 RES: 0.0010747
Jacobi stage: ITER: 43035 RES: 0.00097285
Jacobi stage: ITER: 43681 RES: 0.00088068
Jacobi stage: ITER: 44327 RES: 0.00079774
Jacobi stage: ITER: 44973 RES: 0.00072216
Jacobi stage: ITER: 45619 RES: 0.00065414
Jacobi stage: ITER: 46265 RES: 0.00059217
Jacobi stage: ITER: 46911 RES: 0.0005364
Jacobi stage: ITER: 47561 RES: 0.00048528
Jacobi stage: ITER: 48213 RES: 0.00043903
Jacobi stage: ITER: 48865 RES: 0.00039719
Jacobi stage: ITER: 49515 RES: 0.00035956
Jacobi stage: ITER: 50167 RES: 0.0003253
Jacobi stage: ITER: 50819 RES: 0.0002943
Jacobi stage: ITER: 51471 RES: 0.00026625
Jacobi stage: ITER: 52123 RES: 0.00024088
Jacobi stage: ITER: 52775 RES: 0.00021793
Jacobi stage: ITER: 53425 RES: 0.00019716
Jacobi stage: ITER: 54079 RES: 0.00017837
Jacobi stage: ITER: 54729 RES: 0.00016137
Jacobi stage: ITER: 55381 RES: 0.00014599
Jacobi stage: ITER: 56031 RES: 0.00013216
Jacobi stage: ITER: 56683 RES: 0.00011957
Jacobi stage: ITER: 57335 RES: 0.00010817
Jacobi stage: ITER: 57987 RES: 9.7864e-05
Jacobi stage: ITER: 58639 RES: 8.8538e-05
Jacobi stage: ITER: 59293 RES: 8.0051e-05
Jacobi stage: ITER: 59945 RES: 7.2423e-05
Jacobi stage: ITER: 60597 RES: 6.5521e-05
Jacobi stage: ITER: 61249 RES: 5.9277e-05
Jacobi stage: ITER: 61899 RES: 5.3661e-05
Jacobi stage: ITER: 62551 RES: 4.8547e-05
Jacobi stage: ITER: 63203 RES: 4.3921e-05
Jacobi stage: ITER: 63855 RES: 3.9736e-05
Jacobi stage: ITER: 64509 RES: 3.5927e-05
Jacobi stage: ITER: 65161 RES: 3.2503e-05
Jacobi stage: ITER: 65813 RES: 2.9406e-05
Jacobi stage: ITER: 66465 RES: 2.6603e-05
Jacobi stage: ITER: 67117 RES: 2.4068e-05
Jacobi stage: ITER: 67769 RES: 2.1775e-05
Jacobi stage: ITER: 68421 RES: 1.9699e-05
Jacobi stage: ITER: 69075 RES: 1.7822e-05
Jacobi stage: ITER: 69727 RES: 1.6124e-05
Jacobi stage: ITER: 70373 RES: 1.4596e-05
Jacobi stage: ITER: 71019 RES: 1.3221e-05
Jacobi stage: ITER: 71667 RES: 1.1969e-05
Jacobi stage: ITER: 72313 RES: 1.0835e-05
Jacobi stage: ITER: 72959 RES: 9.8144e-06
Jacobi stage: ITER: 73605 RES: 8.8846e-06
Jacobi stage: ITER: 74251 RES: 8.0478e-06
Jacobi stage: ITER: 74899 RES: 7.2854e-06
Jacobi stage: ITER: 75547 RES: 6.5952e-06
Jacobi stage: ITER: 76193 RES: 5.9703e-06
Jacobi stage: ITER: 76839 RES: 5.408e-06
Jacobi stage: ITER: 77489 RES: 4.8926e-06
Jacobi stage: ITER: 78141 RES: 4.4264e-06
Jacobi stage: ITER: 78795 RES: 4.0046e-06
Jacobi stage: ITER: 79445 RES: 3.6229e-06
Jacobi stage: ITER: 80099 RES: 3.2777e-06
Jacobi stage: ITER: 80751 RES: 2.9653e-06
Jacobi stage: ITER: 81403 RES: 2.6828e-06
Jacobi stage: ITER: 82055 RES: 2.4271e-06
Jacobi stage: ITER: 82707 RES: 2.1958e-06
Jacobi stage: ITER: 83359 RES: 1.9865e-06
Jacobi stage: ITER: 84011 RES: 1.7972e-06
Jacobi stage: ITER: 84663 RES: 1.626e-06
Jacobi stage: ITER: 85315 RES: 1.471e-06
Jacobi stage: ITER: 85967 RES: 1.3308e-06
Jacobi stage: ITER: 86619 RES: 1.204e-06
Jacobi stage: ITER: 87271 RES: 1.0893e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]

The monitoring of the solver can be improved by time elapsed since the beginning of the computation as demonstrated in the following example:

1#include <iostream>
2#include <memory>
3#include <chrono>
4#include <thread>
5#include <TNL/Timer.h>
6#include <TNL/Matrices/SparseMatrix.h>
7#include <TNL/Devices/Sequential.h>
8#include <TNL/Devices/Cuda.h>
9#include <TNL/Solvers/Linear/Jacobi.h>
10
11template< typename Device >
12void
13iterativeLinearSolverExample()
14{
15 /***
16 * Set the following matrix (dots represent zero matrix elements):
17 *
18 * / 2.5 -1 . . . \
19 * | -1 2.5 -1 . . |
20 * | . -1 2.5 -1. . |
21 * | . . -1 2.5 -1 |
22 * \ . . . -1 2.5 /
23 */
24 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
25 using Vector = TNL::Containers::Vector< double, Device >;
26 const int size( 5 );
27 auto matrix_ptr = std::make_shared< MatrixType >();
28 matrix_ptr->setDimensions( size, size );
29 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
30
31 auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
32 {
33 const int rowIdx = row.getRowIndex();
34 if( rowIdx == 0 ) {
35 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
36 row.setElement( 1, rowIdx + 1, -1 ); // element above the diagonal
37 }
38 else if( rowIdx == size - 1 ) {
39 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
40 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
41 }
42 else {
43 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
44 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
45 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
46 }
47 };
48
49 /***
50 * Set the matrix elements.
51 */
52 matrix_ptr->forAllRows( f );
53 std::cout << *matrix_ptr << std::endl;
54
55 /***
56 * Set the right-hand side vector.
57 */
58 Vector x( size, 1.0 );
59 Vector b( size );
60 matrix_ptr->vectorProduct( x, b );
61 x = 0.0;
62 std::cout << "Vector b = " << b << std::endl;
63
64 /***
65 * Setup solver of the linear system.
66 */
67 using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
68 LinearSolver solver;
69 solver.setMatrix( matrix_ptr );
70 solver.setOmega( 0.0005 );
71
72 /***
73 * Setup monitor of the iterative solver.
74 */
75 using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
76 IterativeSolverMonitorType monitor;
77 TNL::Solvers::SolverMonitorThread mmonitorThread( monitor );
78 monitor.setRefreshRate( 10 ); // refresh rate in milliseconds
79 monitor.setVerbose( 1 );
80 monitor.setStage( "Jacobi stage:" );
81 TNL::Timer timer;
82 monitor.setTimer( timer );
83 timer.start();
84 solver.setSolverMonitor( monitor );
85 solver.setConvergenceResidue( 1.0e-6 );
86 solver.solve( b, x );
87 monitor.stopMainLoop();
88 std::cout << "Vector x = " << x << std::endl;
89}
90
91int
92main( int argc, char* argv[] )
93{
94 std::cout << "Solving linear system on host: " << std::endl;
95 iterativeLinearSolverExample< TNL::Devices::Sequential >();
96
97#ifdef __CUDACC__
98 std::cout << "Solving linear system on CUDA device: " << std::endl;
99 iterativeLinearSolverExample< TNL::Devices::Cuda >();
100#endif
101}
TNL::Timer
Class for real time, CPU time and CPU cycles measuring.
Definition Timer.h:25
TNL::Timer::start
void start()
Starts timer.
Definition Timer.hpp:42

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:4.651e-0
ELA: 0.1001 Jacobi stage: ITER: 41375 RES: 0.0012546
ELA: 0.20022 Jacobi stage: ITER: 81273 RES: 2.736e-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:2.8301e-
ELA: 0.10013 Jacobi stage: ITER: 735 RES: 0.79433
ELA: 0.20023 Jacobi stage: ITER: 1511 RES: 0.64986
ELA: 0.30034 Jacobi stage: ITER: 2097 RES: 0.57012
ELA: 0.40045 Jacobi stage: ITER: 2711 RES: 0.50423
ELA: 0.50056 Jacobi stage: ITER: 3345 RES: 0.4483
ELA: 0.60069 Jacobi stage: ITER: 3965 RES: 0.40242
ELA: 0.70081 Jacobi stage: ITER: 4603 RES: 0.36187
ELA: 0.801 Jacobi stage: ITER: 5223 RES: 0.32722
ELA: 0.90117 Jacobi stage: ITER: 5851 RES: 0.29607
ELA: 1.0013 Jacobi stage: ITER: 6487 RES: 0.26788
ELA: 1.1014 Jacobi stage: ITER: 7119 RES: 0.24274
ELA: 1.2016 Jacobi stage: ITER: 7741 RES: 0.22034
ELA: 1.3017 Jacobi stage: ITER: 8379 RES: 0.1997
ELA: 1.4019 Jacobi stage: ITER: 9013 RES: 0.18104
ELA: 1.502 Jacobi stage: ITER: 9623 RES: 0.16486
ELA: 1.6022 Jacobi stage: ITER: 10269 RES: 0.14922
ELA: 1.7023 Jacobi stage: ITER: 10907 RES: 0.13531
ELA: 1.8024 Jacobi stage: ITER: 11535 RES: 0.12286
ELA: 1.9025 Jacobi stage: ITER: 12119 RES: 0.11232
ELA: 2.0026 Jacobi stage: ITER: 12757 RES: 0.1018
ELA: 2.1027 Jacobi stage: ITER: 13409 RES: 0.092094
ELA: 2.2028 Jacobi stage: ITER: 14059 RES: 0.083368
ELA: 2.3029 Jacobi stage: ITER: 14711 RES: 0.075422
ELA: 2.403 Jacobi stage: ITER: 15363 RES: 0.068234
ELA: 2.5031 Jacobi stage: ITER: 16013 RES: 0.061732
ELA: 2.6032 Jacobi stage: ITER: 16665 RES: 0.055849
ELA: 2.7033 Jacobi stage: ITER: 17315 RES: 0.050558
ELA: 2.8034 Jacobi stage: ITER: 17967 RES: 0.04574
ELA: 2.9035 Jacobi stage: ITER: 18619 RES: 0.041381
ELA: 3.0036 Jacobi stage: ITER: 19271 RES: 0.037437
ELA: 3.1037 Jacobi stage: ITER: 19921 RES: 0.03387
ELA: 3.2038 Jacobi stage: ITER: 20571 RES: 0.030661
ELA: 3.3039 Jacobi stage: ITER: 21223 RES: 0.027739
ELA: 3.404 Jacobi stage: ITER: 21875 RES: 0.025095
ELA: 3.5041 Jacobi stage: ITER: 22525 RES: 0.022704
ELA: 3.6042 Jacobi stage: ITER: 23175 RES: 0.020553
ELA: 3.7043 Jacobi stage: ITER: 23821 RES: 0.018606
ELA: 3.8045 Jacobi stage: ITER: 24467 RES: 0.016853
ELA: 3.9045 Jacobi stage: ITER: 25111 RES: 0.015266
ELA: 4.0046 Jacobi stage: ITER: 25759 RES: 0.01382
ELA: 4.1047 Jacobi stage: ITER: 26405 RES: 0.01251
ELA: 4.2048 Jacobi stage: ITER: 27051 RES: 0.011332
ELA: 4.3049 Jacobi stage: ITER: 27697 RES: 0.010259
ELA: 4.405 Jacobi stage: ITER: 28343 RES: 0.0092923
ELA: 4.5051 Jacobi stage: ITER: 28991 RES: 0.008412
ELA: 4.6052 Jacobi stage: ITER: 29639 RES: 0.007615
ELA: 4.7053 Jacobi stage: ITER: 30283 RES: 0.0068978
ELA: 4.8054 Jacobi stage: ITER: 30931 RES: 0.0062443
ELA: 4.9055 Jacobi stage: ITER: 31579 RES: 0.0056527
ELA: 5.0056 Jacobi stage: ITER: 32225 RES: 0.0051172
ELA: 5.1057 Jacobi stage: ITER: 32871 RES: 0.0046352
ELA: 5.2058 Jacobi stage: ITER: 33519 RES: 0.0041961
ELA: 5.3059 Jacobi stage: ITER: 34165 RES: 0.0037985
ELA: 5.406 Jacobi stage: ITER: 34811 RES: 0.0034408
ELA: 5.5061 Jacobi stage: ITER: 35457 RES: 0.0031148
ELA: 5.6062 Jacobi stage: ITER: 36103 RES: 0.0028214
ELA: 5.7063 Jacobi stage: ITER: 36749 RES: 0.0025541
ELA: 5.8064 Jacobi stage: ITER: 37395 RES: 0.0023136
ELA: 5.9065 Jacobi stage: ITER: 38047 RES: 0.0020931
ELA: 6.0066 Jacobi stage: ITER: 38699 RES: 0.0018936
ELA: 6.1067 Jacobi stage: ITER: 39351 RES: 0.0017132
ELA: 6.2068 Jacobi stage: ITER: 40001 RES: 0.0015499
ELA: 6.3069 Jacobi stage: ITER: 40653 RES: 0.0014022
ELA: 6.407 Jacobi stage: ITER: 41305 RES: 0.0012686
ELA: 6.507 Jacobi stage: ITER: 41957 RES: 0.0011477
ELA: 6.6071 Jacobi stage: ITER: 42607 RES: 0.001039
ELA: 6.7072 Jacobi stage: ITER: 43259 RES: 0.00093995
ELA: 6.8073 Jacobi stage: ITER: 43911 RES: 0.00085037
ELA: 6.9074 Jacobi stage: ITER: 44563 RES: 0.00076934
ELA: 7.0075 Jacobi stage: ITER: 45215 RES: 0.00069602
ELA: 7.1076 Jacobi stage: ITER: 45867 RES: 0.00062969
ELA: 7.2077 Jacobi stage: ITER: 46519 RES: 0.00056968
ELA: 7.3078 Jacobi stage: ITER: 47171 RES: 0.0005154
ELA: 7.4079 Jacobi stage: ITER: 47823 RES: 0.00046628
ELA: 7.508 Jacobi stage: ITER: 48475 RES: 0.00042185
ELA: 7.6081 Jacobi stage: ITER: 49127 RES: 0.00038164
ELA: 7.7082 Jacobi stage: ITER: 49779 RES: 0.00034528
ELA: 7.8082 Jacobi stage: ITER: 50431 RES: 0.00031237
ELA: 7.9083 Jacobi stage: ITER: 51083 RES: 0.0002826
ELA: 8.0084 Jacobi stage: ITER: 51735 RES: 0.00025567
ELA: 8.1085 Jacobi stage: ITER: 52387 RES: 0.00023131
ELA: 8.2086 Jacobi stage: ITER: 53039 RES: 0.00020926
ELA: 8.3087 Jacobi stage: ITER: 53691 RES: 0.00018932
ELA: 8.4088 Jacobi stage: ITER: 54343 RES: 0.00017128
ELA: 8.5089 Jacobi stage: ITER: 54995 RES: 0.00015496
ELA: 8.609 Jacobi stage: ITER: 55647 RES: 0.00014019
ELA: 8.709 Jacobi stage: ITER: 56299 RES: 0.00012683
ELA: 8.8092 Jacobi stage: ITER: 56951 RES: 0.00011474
ELA: 8.9093 Jacobi stage: ITER: 57603 RES: 0.00010381
ELA: 9.0094 Jacobi stage: ITER: 58255 RES: 9.3917e-05
ELA: 9.1095 Jacobi stage: ITER: 58907 RES: 8.4967e-05
ELA: 9.2096 Jacobi stage: ITER: 59559 RES: 7.687e-05
ELA: 9.3097 Jacobi stage: ITER: 60211 RES: 6.9545e-05
ELA: 9.4098 Jacobi stage: ITER: 60859 RES: 6.2956e-05
ELA: 9.5099 Jacobi stage: ITER: 61505 RES: 5.6991e-05
ELA: 9.61 Jacobi stage: ITER: 62151 RES: 5.1624e-05
ELA: 9.7101 Jacobi stage: ITER: 62799 RES: 4.6733e-05
ELA: 9.8102 Jacobi stage: ITER: 63445 RES: 4.2305e-05
ELA: 9.9103 Jacobi stage: ITER: 64091 RES: 3.8321e-05
ELA: 10.01 Jacobi stage: ITER: 64735 RES: 3.4712e-05
ELA: 10.111 Jacobi stage: ITER: 65383 RES: 3.1423e-05
ELA: 10.211 Jacobi stage: ITER: 66031 RES: 2.8446e-05
ELA: 10.311 Jacobi stage: ITER: 66677 RES: 2.5751e-05
ELA: 10.411 Jacobi stage: ITER: 67323 RES: 2.3326e-05
ELA: 10.511 Jacobi stage: ITER: 67975 RES: 2.1103e-05
ELA: 10.611 Jacobi stage: ITER: 68627 RES: 1.9092e-05
ELA: 10.711 Jacobi stage: ITER: 69279 RES: 1.7272e-05
ELA: 10.811 Jacobi stage: ITER: 69931 RES: 1.5626e-05
ELA: 10.911 Jacobi stage: ITER: 70583 RES: 1.4137e-05
ELA: 11.011 Jacobi stage: ITER: 71235 RES: 1.279e-05
ELA: 11.111 Jacobi stage: ITER: 71887 RES: 1.1571e-05
ELA: 11.212 Jacobi stage: ITER: 72539 RES: 1.0468e-05
ELA: 11.312 Jacobi stage: ITER: 73191 RES: 9.4709e-06
ELA: 11.412 Jacobi stage: ITER: 73843 RES: 8.5683e-06
ELA: 11.512 Jacobi stage: ITER: 74495 RES: 7.7518e-06
ELA: 11.612 Jacobi stage: ITER: 75147 RES: 7.0131e-06
ELA: 11.712 Jacobi stage: ITER: 75799 RES: 6.3447e-06
ELA: 11.812 Jacobi stage: ITER: 76453 RES: 5.7366e-06
ELA: 11.912 Jacobi stage: ITER: 77103 RES: 5.1931e-06
ELA: 12.012 Jacobi stage: ITER: 77757 RES: 4.6953e-06
ELA: 12.112 Jacobi stage: ITER: 78505 RES: 4.1857e-06
ELA: 12.212 Jacobi stage: ITER: 79497 RES: 3.5941e-06
ELA: 12.313 Jacobi stage: ITER: 80491 RES: 3.0862e-06
ELA: 12.413 Jacobi stage: ITER: 81485 RES: 2.6484e-06
ELA: 12.513 Jacobi stage: ITER: 82479 RES: 2.2741e-06
ELA: 12.613 Jacobi stage: ITER: 83475 RES: 1.9515e-06
ELA: 12.713 Jacobi stage: ITER: 84471 RES: 1.6746e-06
ELA: 12.813 Jacobi stage: ITER: 85467 RES: 1.4371e-06
ELA: 12.913 Jacobi stage: ITER: 86461 RES: 1.2332e-06
ELA: 13.013 Jacobi stage: ITER: 87481 RES: 1.0544e-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 */
22 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
23 using Vector = TNL::Containers::Vector< double, Device >;
24 const int size( 5 );
25 auto matrix_ptr = std::make_shared< MatrixType >();
26 matrix_ptr->setDimensions( size, size );
27 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
28
29 auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
30 {
31 const int rowIdx = row.getRowIndex();
32 if( rowIdx == 0 ) {
33 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
34 row.setElement( 1, rowIdx + 1, -1 ); // element above the diagonal
35 }
36 else if( rowIdx == size - 1 ) {
37 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
38 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
39 }
40 else {
41 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
42 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
43 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
44 }
45 };
46
47 /***
48 * Set the matrix elements.
49 */
50 matrix_ptr->forAllRows( f );
51 std::cout << *matrix_ptr << std::endl;
52
53 /***
54 * Set the right-hand side vector.
55 */
56 Vector x( size, 1.0 );
57 Vector b( size );
58 matrix_ptr->vectorProduct( x, b );
59 x = 0.0;
60 std::cout << "Vector b = " << b << std::endl;
61
62 /***
63 * Solve the linear system using diagonal (Jacobi) preconditioner.
64 */
65 using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
66 using Preconditioner = TNL::Solvers::Linear::Preconditioners::Diagonal< MatrixType >;
67 auto preconditioner_ptr = std::make_shared< Preconditioner >();
68 preconditioner_ptr->update( matrix_ptr );
69 LinearSolver solver;
70 solver.setMatrix( matrix_ptr );
71 solver.setPreconditioner( preconditioner_ptr );
72 solver.setConvergenceResidue( 1.0e-6 );
73 solver.solve( b, x );
74 std::cout << "Vector x = " << x << std::endl;
75}
76
77int
78main( int argc, char* argv[] )
79{
80 std::cout << "Solving linear system on host: " << std::endl;
81 iterativeLinearSolverExample< TNL::Devices::Sequential >();
82
83#ifdef __CUDACC__
84 std::cout << "Solving linear system on CUDA device: " << std::endl;
85 iterativeLinearSolverExample< TNL::Devices::Cuda >();
86#endif
87}
TNL::Solvers::Linear::Preconditioners::Diagonal
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 */
21 using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
22 using Vector = TNL::Containers::Vector< double, Device >;
23 const int size( 5 );
24 auto matrix_ptr = std::make_shared< MatrixType >();
25 matrix_ptr->setDimensions( size, size );
26 matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
27
28 auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
29 {
30 const int rowIdx = row.getRowIndex();
31 if( rowIdx == 0 ) {
32 row.setElement( 0, rowIdx, 2.5 ); // diagonal element
33 row.setElement( 1, rowIdx + 1, -1 ); // element above the diagonal
34 }
35 else if( rowIdx == size - 1 ) {
36 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
37 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
38 }
39 else {
40 row.setElement( 0, rowIdx - 1, -1.0 ); // element below the diagonal
41 row.setElement( 1, rowIdx, 2.5 ); // diagonal element
42 row.setElement( 2, rowIdx + 1, -1.0 ); // element above the diagonal
43 }
44 };
45
46 /***
47 * Set the matrix elements.
48 */
49 matrix_ptr->forAllRows( f );
50 std::cout << *matrix_ptr << std::endl;
51
52 /***
53 * Set the right-hand side vector.
54 */
55 Vector x( size, 1.0 );
56 Vector b( size );
57 matrix_ptr->vectorProduct( x, b );
58 x = 0.0;
59 std::cout << "Vector b = " << b << std::endl;
60
61 /***
62 * Solve the linear system using diagonal (Jacobi) preconditioner.
63 */
64 auto solver_ptr = TNL::Solvers::getLinearSolver< MatrixType >( "tfqmr" );
65 auto preconditioner_ptr = TNL::Solvers::getPreconditioner< MatrixType >( "diagonal" );
66 preconditioner_ptr->update( matrix_ptr );
67 solver_ptr->setMatrix( matrix_ptr );
68 solver_ptr->setPreconditioner( preconditioner_ptr );
69 solver_ptr->setConvergenceResidue( 1.0e-6 );
70 solver_ptr->solve( b, x );
71 std::cout << "Vector x = " << x << std::endl;
72}
73
74int
75main( int argc, char* argv[] )
76{
77 std::cout << "Solving linear system on host: " << std::endl;
78 iterativeLinearSolverExample< TNL::Devices::Sequential >();
79
80#ifdef __CUDACC__
81 std::cout << "Solving linear system on CUDA device: " << std::endl;
82 iterativeLinearSolverExample< TNL::Devices::Cuda >();
83#endif
84}
TNL::Solvers::getLinearSolver
std::shared_ptr< Linear::LinearSolver< MatrixType > > getLinearSolver(std::string name)
Function returning shared pointer with linear solver given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:84
TNL::Solvers::getPreconditioner
std::shared_ptr< Linear::Preconditioners::Preconditioner< MatrixType > > getPreconditioner(std::string name)
Function returning shared pointer with linear preconditioner given by its name in a form of a string.
Definition LinearSolverTypeResolver.h:138

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 Fri Jul 19 2024 12:11:20 for Template Numerical Library by doxygen 1.11.0