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

Introduction

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

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

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

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

Iterative solvers of linear systems

Basic setup

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

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

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: 21769 RES: 0.0255
Jacobi stage: ITER: 44555 RES: 0.00077028
Jacobi stage: ITER: 68441 RES: 1.9639e-05
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: 759 RES: 0.78894
Jacobi stage: ITER: 1819 RES: 0.60566
Jacobi stage: ITER: 2863 RES: 0.48987
Jacobi stage: ITER: 3791 RES: 0.41473
Jacobi stage: ITER: 4511 RES: 0.36739
Jacobi stage: ITER: 5187 RES: 0.32912
Jacobi stage: ITER: 5871 RES: 0.29513
Jacobi stage: ITER: 6395 RES: 0.27177
Jacobi stage: ITER: 7035 RES: 0.24593
Jacobi stage: ITER: 7471 RES: 0.22982
Jacobi stage: ITER: 8171 RES: 0.20622
Jacobi stage: ITER: 9307 RES: 0.17308
Jacobi stage: ITER: 10053 RES: 0.15426
Jacobi stage: ITER: 10657 RES: 0.14057
Jacobi stage: ITER: 11285 RES: 0.12763
Jacobi stage: ITER: 11927 RES: 0.11568
Jacobi stage: ITER: 12361 RES: 0.10818
Jacobi stage: ITER: 12727 RES: 0.1023
Jacobi stage: ITER: 13139 RES: 0.096023
Jacobi stage: ITER: 13771 RES: 0.087139
Jacobi stage: ITER: 14281 RES: 0.080548
Jacobi stage: ITER: 14609 RES: 0.07659
Jacobi stage: ITER: 15187 RES: 0.070104
Jacobi stage: ITER: 15567 RES: 0.066129
Jacobi stage: ITER: 16179 RES: 0.060196
Jacobi stage: ITER: 16857 RES: 0.054226
Jacobi stage: ITER: 17591 RES: 0.048459
Jacobi stage: ITER: 18271 RES: 0.043653
Jacobi stage: ITER: 18711 RES: 0.0408
Jacobi stage: ITER: 19421 RES: 0.036573
Jacobi stage: ITER: 19975 RES: 0.0336
Jacobi stage: ITER: 20671 RES: 0.030193
Jacobi stage: ITER: 21265 RES: 0.027552
Jacobi stage: ITER: 21865 RES: 0.025126
Jacobi stage: ITER: 22257 RES: 0.023658
Jacobi stage: ITER: 22731 RES: 0.022004
Jacobi stage: ITER: 23359 RES: 0.01998
Jacobi stage: ITER: 23811 RES: 0.01864
Jacobi stage: ITER: 24439 RES: 0.016926
Jacobi stage: ITER: 24803 RES: 0.016006
Jacobi stage: ITER: 25211 RES: 0.015033
Jacobi stage: ITER: 26167 RES: 0.01298
Jacobi stage: ITER: 26807 RES: 0.011765
Jacobi stage: ITER: 27539 RES: 0.010514
Jacobi stage: ITER: 28083 RES: 0.009671
Jacobi stage: ITER: 28711 RES: 0.0087817
Jacobi stage: ITER: 29367 RES: 0.0079399
Jacobi stage: ITER: 30013 RES: 0.0071877
Jacobi stage: ITER: 30753 RES: 0.0064154
Jacobi stage: ITER: 31503 RES: 0.0057191
Jacobi stage: ITER: 31803 RES: 0.0054615
Jacobi stage: ITER: 32279 RES: 0.0050765
Jacobi stage: ITER: 32927 RES: 0.0045955
Jacobi stage: ITER: 33547 RES: 0.0041781
Jacobi stage: ITER: 34165 RES: 0.0037985
Jacobi stage: ITER: 34795 RES: 0.0034492
Jacobi stage: ITER: 35427 RES: 0.0031301
Jacobi stage: ITER: 36059 RES: 0.0028406
Jacobi stage: ITER: 36687 RES: 0.0025794
Jacobi stage: ITER: 37319 RES: 0.0023407
Jacobi stage: ITER: 37951 RES: 0.0021242
Jacobi stage: ITER: 38581 RES: 0.0019277
Jacobi stage: ITER: 39213 RES: 0.0017493
Jacobi stage: ITER: 39845 RES: 0.0015875
Jacobi stage: ITER: 40477 RES: 0.0014406
Jacobi stage: ITER: 41107 RES: 0.0013082
Jacobi stage: ITER: 41739 RES: 0.0011871
Jacobi stage: ITER: 42371 RES: 0.0010773
Jacobi stage: ITER: 42999 RES: 0.00097825
Jacobi stage: ITER: 43631 RES: 0.00088774
Jacobi stage: ITER: 44271 RES: 0.00080463
Jacobi stage: ITER: 44901 RES: 0.00073019
Jacobi stage: ITER: 45531 RES: 0.00066304
Jacobi stage: ITER: 46171 RES: 0.00060096
Jacobi stage: ITER: 46801 RES: 0.00054537
Jacobi stage: ITER: 47433 RES: 0.00049491
Jacobi stage: ITER: 48063 RES: 0.0004494
Jacobi stage: ITER: 48695 RES: 0.00040783
Jacobi stage: ITER: 49325 RES: 0.0003701
Jacobi stage: ITER: 49975 RES: 0.00033504
Jacobi stage: ITER: 50607 RES: 0.00030404
Jacobi stage: ITER: 51237 RES: 0.00027591
Jacobi stage: ITER: 51867 RES: 0.00025054
Jacobi stage: ITER: 52499 RES: 0.00022736
Jacobi stage: ITER: 53131 RES: 0.00020633
Jacobi stage: ITER: 53763 RES: 0.00018724
Jacobi stage: ITER: 54395 RES: 0.00016992
Jacobi stage: ITER: 55023 RES: 0.00015429
Jacobi stage: ITER: 55653 RES: 0.00014002
Jacobi stage: ITER: 56297 RES: 0.00012683
Jacobi stage: ITER: 56927 RES: 0.00011517
Jacobi stage: ITER: 57559 RES: 0.00010451
Jacobi stage: ITER: 58189 RES: 9.4845e-05
Jacobi stage: ITER: 58819 RES: 8.6124e-05
Jacobi stage: ITER: 59451 RES: 7.8156e-05
Jacobi stage: ITER: 60081 RES: 7.0926e-05
Jacobi stage: ITER: 60713 RES: 6.4364e-05
Jacobi stage: ITER: 61345 RES: 5.8409e-05
Jacobi stage: ITER: 61979 RES: 5.3006e-05
Jacobi stage: ITER: 62613 RES: 4.8073e-05
Jacobi stage: ITER: 63251 RES: 4.3598e-05
Jacobi stage: ITER: 63887 RES: 3.9541e-05
Jacobi stage: ITER: 64519 RES: 3.5883e-05
Jacobi stage: ITER: 65155 RES: 3.2543e-05
Jacobi stage: ITER: 65787 RES: 2.9532e-05
Jacobi stage: ITER: 66417 RES: 2.68e-05
Jacobi stage: ITER: 67051 RES: 2.4321e-05
Jacobi stage: ITER: 67683 RES: 2.2071e-05
Jacobi stage: ITER: 68315 RES: 2.0029e-05
Jacobi stage: ITER: 68949 RES: 1.8165e-05
Jacobi stage: ITER: 69581 RES: 1.6484e-05
Jacobi stage: ITER: 70215 RES: 1.4959e-05
Jacobi stage: ITER: 70849 RES: 1.3567e-05
Jacobi stage: ITER: 71483 RES: 1.2312e-05
Jacobi stage: ITER: 72115 RES: 1.1173e-05
Jacobi stage: ITER: 72747 RES: 1.0139e-05
Jacobi stage: ITER: 73379 RES: 9.2013e-06
Jacobi stage: ITER: 74025 RES: 8.3295e-06
Jacobi stage: ITER: 74659 RES: 7.559e-06
Jacobi stage: ITER: 75291 RES: 6.8597e-06
Jacobi stage: ITER: 75923 RES: 6.225e-06
Jacobi stage: ITER: 76557 RES: 5.6457e-06
Jacobi stage: ITER: 77191 RES: 5.1234e-06
Jacobi stage: ITER: 77823 RES: 4.6494e-06
Jacobi stage: ITER: 78455 RES: 4.2193e-06
Jacobi stage: ITER: 79085 RES: 3.8289e-06
Jacobi stage: ITER: 79719 RES: 3.4747e-06
Jacobi stage: ITER: 80353 RES: 3.1513e-06
Jacobi stage: ITER: 80987 RES: 2.8598e-06
Jacobi stage: ITER: 81623 RES: 2.5936e-06
Jacobi stage: ITER: 82259 RES: 2.3522e-06
Jacobi stage: ITER: 82891 RES: 2.1346e-06
Jacobi stage: ITER: 83523 RES: 1.9371e-06
Jacobi stage: ITER: 84153 RES: 1.7579e-06
Jacobi stage: ITER: 84785 RES: 1.5953e-06
Jacobi stage: ITER: 85415 RES: 1.4486e-06
Jacobi stage: ITER: 86047 RES: 1.3146e-06
Jacobi stage: ITER: 86679 RES: 1.193e-06
Jacobi stage: ITER: 87315 RES: 1.0819e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]

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

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

The only changes 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.905e-0
ELA: 0.10191 Jacobi stage: ITER: 24555 RES: 0.016617
ELA: 0.20207 Jacobi stage: ITER: 47193 RES: 0.0005135
ELA: 0.30253 Jacobi stage: ITER: 69835 RES: 1.5859e-05
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:0.002498
ELA: 0.10504 Jacobi stage: ITER: 1231 RES: 0.69572
ELA: 0.20515 Jacobi stage: ITER: 2175 RES: 0.56113
ELA: 0.30531 Jacobi stage: ITER: 3079 RES: 0.4706
ELA: 0.40548 Jacobi stage: ITER: 4245 RES: 0.38386
ELA: 0.50917 Jacobi stage: ITER: 5031 RES: 0.3375
ELA: 0.60933 Jacobi stage: ITER: 5697 RES: 0.30328
ELA: 0.71272 Jacobi stage: ITER: 6293 RES: 0.27607
ELA: 0.81284 Jacobi stage: ITER: 6859 RES: 0.25276
ELA: 0.91915 Jacobi stage: ITER: 7631 RES: 0.22419
ELA: 1.0193 Jacobi stage: ITER: 8367 RES: 0.20007
ELA: 1.1194 Jacobi stage: ITER: 9019 RES: 0.18093
ELA: 1.2195 Jacobi stage: ITER: 9647 RES: 0.16425
ELA: 1.3197 Jacobi stage: ITER: 10451 RES: 0.14514
ELA: 1.4198 Jacobi stage: ITER: 11087 RES: 0.13162
ELA: 1.5292 Jacobi stage: ITER: 11773 RES: 0.11841
ELA: 1.6325 Jacobi stage: ITER: 12433 RES: 0.10699
ELA: 1.7326 Jacobi stage: ITER: 13141 RES: 0.095964
ELA: 1.8327 Jacobi stage: ITER: 13911 RES: 0.085285
ELA: 1.9358 Jacobi stage: ITER: 14535 RES: 0.077489
ELA: 2.0425 Jacobi stage: ITER: 14975 RES: 0.072425
ELA: 2.1426 Jacobi stage: ITER: 15607 RES: 0.065724
ELA: 2.2428 Jacobi stage: ITER: 16375 RES: 0.058411
ELA: 2.3429 Jacobi stage: ITER: 17183 RES: 0.051593
ELA: 2.4492 Jacobi stage: ITER: 17859 RES: 0.046505
ELA: 2.5493 Jacobi stage: ITER: 18443 RES: 0.042515
ELA: 2.6495 Jacobi stage: ITER: 18735 RES: 0.04065
ELA: 2.7496 Jacobi stage: ITER: 19153 RES: 0.03811
ELA: 2.8497 Jacobi stage: ITER: 19671 RES: 0.035206
ELA: 2.9498 Jacobi stage: ITER: 20271 RES: 0.032107
ELA: 3.05 Jacobi stage: ITER: 20567 RES: 0.03068
ELA: 3.1501 Jacobi stage: ITER: 21091 RES: 0.028307
ELA: 3.2525 Jacobi stage: ITER: 21425 RES: 0.026883
ELA: 3.3526 Jacobi stage: ITER: 21889 RES: 0.025034
ELA: 3.4528 Jacobi stage: ITER: 22539 RES: 0.022662
ELA: 3.553 Jacobi stage: ITER: 23321 RES: 0.020091
ELA: 3.6531 Jacobi stage: ITER: 23979 RES: 0.018165
ELA: 3.7533 Jacobi stage: ITER: 24649 RES: 0.016384
ELA: 3.8558 Jacobi stage: ITER: 25249 RES: 0.014941
ELA: 3.956 Jacobi stage: ITER: 25763 RES: 0.013811
ELA: 4.0561 Jacobi stage: ITER: 26483 RES: 0.012365
ELA: 4.1562 Jacobi stage: ITER: 27165 RES: 0.011132
ELA: 4.2591 Jacobi stage: ITER: 27639 RES: 0.010354
ELA: 4.3593 Jacobi stage: ITER: 28161 RES: 0.0095528
ELA: 4.4692 Jacobi stage: ITER: 28721 RES: 0.0087655
ELA: 4.5725 Jacobi stage: ITER: 29403 RES: 0.0078961
ELA: 4.6726 Jacobi stage: ITER: 30025 RES: 0.0071744
ELA: 4.7728 Jacobi stage: ITER: 30521 RES: 0.0066482
ELA: 4.8729 Jacobi stage: ITER: 30953 RES: 0.0062213
ELA: 4.973 Jacobi stage: ITER: 31283 RES: 0.0059157
ELA: 5.0732 Jacobi stage: ITER: 31679 RES: 0.0055666
ELA: 5.1734 Jacobi stage: ITER: 32443 RES: 0.0049502
ELA: 5.2735 Jacobi stage: ITER: 33083 RES: 0.0044867
ELA: 5.3737 Jacobi stage: ITER: 33671 RES: 0.0040992
ELA: 5.4739 Jacobi stage: ITER: 34287 RES: 0.0037292
ELA: 5.574 Jacobi stage: ITER: 34911 RES: 0.0033883
ELA: 5.6741 Jacobi stage: ITER: 35543 RES: 0.0030749
ELA: 5.7743 Jacobi stage: ITER: 36173 RES: 0.0027904
ELA: 5.8744 Jacobi stage: ITER: 36803 RES: 0.0025338
ELA: 5.9746 Jacobi stage: ITER: 37433 RES: 0.0022994
ELA: 6.0747 Jacobi stage: ITER: 38063 RES: 0.002088
ELA: 6.1748 Jacobi stage: ITER: 38695 RES: 0.0018948
ELA: 6.275 Jacobi stage: ITER: 39327 RES: 0.0017195
ELA: 6.3751 Jacobi stage: ITER: 39959 RES: 0.0015604
ELA: 6.4752 Jacobi stage: ITER: 40591 RES: 0.0014161
ELA: 6.5754 Jacobi stage: ITER: 41219 RES: 0.0012858
ELA: 6.6755 Jacobi stage: ITER: 41849 RES: 0.0011669
ELA: 6.7756 Jacobi stage: ITER: 42481 RES: 0.0010589
ELA: 6.8759 Jacobi stage: ITER: 43111 RES: 0.00096156
ELA: 6.976 Jacobi stage: ITER: 43743 RES: 0.0008726
ELA: 7.0761 Jacobi stage: ITER: 44371 RES: 0.00079236
ELA: 7.1763 Jacobi stage: ITER: 45003 RES: 0.00071906
ELA: 7.2764 Jacobi stage: ITER: 45633 RES: 0.00065254
ELA: 7.3766 Jacobi stage: ITER: 46263 RES: 0.00059253
ELA: 7.4767 Jacobi stage: ITER: 46895 RES: 0.00053772
ELA: 7.5769 Jacobi stage: ITER: 47527 RES: 0.00048797
ELA: 7.677 Jacobi stage: ITER: 48157 RES: 0.00044283
ELA: 7.7771 Jacobi stage: ITER: 48787 RES: 0.00040211
ELA: 7.8773 Jacobi stage: ITER: 49419 RES: 0.00036491
ELA: 7.9774 Jacobi stage: ITER: 50053 RES: 0.00033094
ELA: 8.0776 Jacobi stage: ITER: 50683 RES: 0.00030051
ELA: 8.1777 Jacobi stage: ITER: 51315 RES: 0.00027271
ELA: 8.2779 Jacobi stage: ITER: 51947 RES: 0.00024748
ELA: 8.378 Jacobi stage: ITER: 52577 RES: 0.00022459
ELA: 8.4781 Jacobi stage: ITER: 53211 RES: 0.00020381
ELA: 8.5782 Jacobi stage: ITER: 53841 RES: 0.00018495
ELA: 8.6784 Jacobi stage: ITER: 54471 RES: 0.00016795
ELA: 8.7785 Jacobi stage: ITER: 55103 RES: 0.00015241
ELA: 8.8787 Jacobi stage: ITER: 55733 RES: 0.00013831
ELA: 8.9788 Jacobi stage: ITER: 56363 RES: 0.00012559
ELA: 9.0789 Jacobi stage: ITER: 56995 RES: 0.00011397
ELA: 9.1791 Jacobi stage: ITER: 57625 RES: 0.00010343
ELA: 9.2791 Jacobi stage: ITER: 58255 RES: 9.3917e-05
ELA: 9.3793 Jacobi stage: ITER: 58885 RES: 8.5229e-05
ELA: 9.4794 Jacobi stage: ITER: 59519 RES: 7.7344e-05
ELA: 9.5795 Jacobi stage: ITER: 60147 RES: 7.0232e-05
ELA: 9.6797 Jacobi stage: ITER: 60779 RES: 6.3734e-05
ELA: 9.7798 Jacobi stage: ITER: 61409 RES: 5.7838e-05
ELA: 9.8799 Jacobi stage: ITER: 62037 RES: 5.252e-05
ELA: 9.9825 Jacobi stage: ITER: 62673 RES: 4.7632e-05
ELA: 10.083 Jacobi stage: ITER: 63295 RES: 4.3305e-05
ELA: 10.183 Jacobi stage: ITER: 63919 RES: 3.9347e-05
ELA: 10.283 Jacobi stage: ITER: 64541 RES: 3.5751e-05
ELA: 10.383 Jacobi stage: ITER: 65163 RES: 3.2503e-05
ELA: 10.483 Jacobi stage: ITER: 65787 RES: 2.9532e-05
ELA: 10.583 Jacobi stage: ITER: 66415 RES: 2.6817e-05
ELA: 10.683 Jacobi stage: ITER: 67043 RES: 2.4351e-05
ELA: 10.784 Jacobi stage: ITER: 67667 RES: 2.2125e-05
ELA: 10.884 Jacobi stage: ITER: 68295 RES: 2.0091e-05
ELA: 10.984 Jacobi stage: ITER: 68919 RES: 1.8254e-05
ELA: 11.084 Jacobi stage: ITER: 69543 RES: 1.6586e-05
ELA: 11.184 Jacobi stage: ITER: 70167 RES: 1.507e-05
ELA: 11.284 Jacobi stage: ITER: 70791 RES: 1.3693e-05
ELA: 11.384 Jacobi stage: ITER: 71417 RES: 1.2434e-05
ELA: 11.485 Jacobi stage: ITER: 72047 RES: 1.129e-05
ELA: 11.585 Jacobi stage: ITER: 72675 RES: 1.0252e-05
ELA: 11.685 Jacobi stage: ITER: 73303 RES: 9.3093e-06
ELA: 11.785 Jacobi stage: ITER: 73931 RES: 8.4533e-06
ELA: 11.885 Jacobi stage: ITER: 74557 RES: 7.676e-06
ELA: 11.985 Jacobi stage: ITER: 75183 RES: 6.9744e-06
ELA: 12.085 Jacobi stage: ITER: 75807 RES: 6.337e-06
ELA: 12.186 Jacobi stage: ITER: 76435 RES: 5.7542e-06
ELA: 12.286 Jacobi stage: ITER: 77063 RES: 5.2251e-06
ELA: 12.386 Jacobi stage: ITER: 77691 RES: 4.7446e-06
ELA: 12.486 Jacobi stage: ITER: 78319 RES: 4.3083e-06
ELA: 12.586 Jacobi stage: ITER: 78949 RES: 3.9098e-06
ELA: 12.686 Jacobi stage: ITER: 79575 RES: 3.5524e-06
ELA: 12.786 Jacobi stage: ITER: 80199 RES: 3.2277e-06
ELA: 12.887 Jacobi stage: ITER: 80825 RES: 2.9309e-06
ELA: 12.987 Jacobi stage: ITER: 81447 RES: 2.6647e-06
ELA: 13.087 Jacobi stage: ITER: 82071 RES: 2.4211e-06
ELA: 13.187 Jacobi stage: ITER: 82697 RES: 2.1985e-06
ELA: 13.287 Jacobi stage: ITER: 83327 RES: 1.9963e-06
ELA: 13.387 Jacobi stage: ITER: 83957 RES: 1.8116e-06
ELA: 13.487 Jacobi stage: ITER: 84585 RES: 1.6451e-06
ELA: 13.587 Jacobi stage: ITER: 85215 RES: 1.4938e-06
ELA: 13.688 Jacobi stage: ITER: 85845 RES: 1.3556e-06
ELA: 13.788 Jacobi stage: ITER: 86473 RES: 1.2309e-06
ELA: 13.888 Jacobi stage: ITER: 87099 RES: 1.1184e-06
ELA: 13.988 Jacobi stage: ITER: 87727 RES: 1.0156e-06
Vector x = [ 0.999999, 0.999999, 0.999998, 0.999999, 0.999999 ]

Setup with preconditioner

Preconditioners of iterative solvers can significantly improve the performance of the solver. In the case of the linear systems, they are used mainly with the Krylov subspace methods. Preconditioners cannot be used with the starionary methods (TNL::Solvers::Linear::Jacobi and TNL::Solvers::Linear::SOR). The following example shows how to setup an iterative solver of linear systems with preconditioning.

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

In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen 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 );
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}
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
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 ]