Introduction

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

Stationary methods
1. Jacobi method (TNL::Solvers::Linear::Jacobi)
2. Successive-overrelaxation method, SOR (TNL::Solvers::Linear::SOR)
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:

Diagonal or Jacobi (TNL::Solvers::Linear::Preconditioners::Diagonal)
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:

#include <iostream>
#include <memory>
#include <TNL/Matrices/SparseMatrix.h>
#include <TNL/Devices/Host.h>
#include <TNL/Devices/Cuda.h>
#include <TNL/Solvers/Linear/TFQMR.h>
 
template< typename Device >
void
iterativeLinearSolverExample()
{
   /***
    * Set the following matrix (dots represent zero matrix elements):
    *
    *   /  2.5 -1    .    .    .   \
    *   | -1    2.5 -1    .    .   |
    *   |  .   -1    2.5 -1.   .   |
    *   |  .    .   -1    2.5 -1   |
    *   \  .    .    .   -1    2.5 /
    */
   using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
   using Vector = TNL::Containers::Vector< double, Device >;
   const int size( 5 );
   auto matrix_ptr = std::make_shared< MatrixType >();
   matrix_ptr->setDimensions( size, size );
   matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
 
   auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
   {
      const int rowIdx = row.getRowIndex();
      if( rowIdx == 0 ) {
         row.setElement( 0, rowIdx, 2.5 );     // diagonal element
         row.setElement( 1, rowIdx + 1, -1 );  // element above the diagonal
      }
      else if( rowIdx == size - 1 ) {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
      }
      else {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
         row.setElement( 2, rowIdx + 1, -1.0 );  // element above the diagonal
      }
   };
 
   /***
    * Set the matrix elements.
    */
   matrix_ptr->forAllRows( f );
   std::cout << *matrix_ptr << std::endl;
 
   /***
    * Set the right-hand side vector.
    */
   Vector x( size, 1.0 );
   Vector b( size );
   matrix_ptr->vectorProduct( x, b );
   x = 0.0;
   std::cout << "Vector b = " << b << std::endl;
 
   /***
    * Solve the linear system.
    */
   using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
   LinearSolver solver;
   solver.setMatrix( matrix_ptr );
   solver.setConvergenceResidue( 1.0e-6 );
   solver.solve( b, x );
   std::cout << "Vector x = " << x << std::endl;
}
 
int
main( int argc, char* argv[] )
{
   std::cout << "Solving linear system on host: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Sequential >();
 
#ifdef __CUDACC__
   std::cout << "Solving linear system on CUDA device: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Cuda >();
#endif
}

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.

#include <iostream>
#include <memory>
#include <chrono>
#include <thread>
#include <TNL/Matrices/SparseMatrix.h>
#include <TNL/Devices/Sequential.h>
#include <TNL/Devices/Cuda.h>
#include <TNL/Solvers/Linear/Jacobi.h>
 
template< typename Device >
void
iterativeLinearSolverExample()
{
   /***
    * Set the following matrix (dots represent zero matrix elements):
    *
    *   /  2.5 -1    .    .    .   \
    *   | -1    2.5 -1    .    .   |
    *   |  .   -1    2.5 -1.   .   |
    *   |  .    .   -1    2.5 -1   |
    *   \  .    .    .   -1    2.5 /
    */
   using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
   using Vector = TNL::Containers::Vector< double, Device >;
   const int size( 5 );
   auto matrix_ptr = std::make_shared< MatrixType >();
   matrix_ptr->setDimensions( size, size );
   matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
 
   auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
   {
      const int rowIdx = row.getRowIndex();
      if( rowIdx == 0 ) {
         row.setElement( 0, rowIdx, 2.5 );     // diagonal element
         row.setElement( 1, rowIdx + 1, -1 );  // element above the diagonal
      }
      else if( rowIdx == size - 1 ) {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
      }
      else {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
         row.setElement( 2, rowIdx + 1, -1.0 );  // element above the diagonal
      }
   };
 
   /***
    * Set the matrix elements.
    */
   matrix_ptr->forAllRows( f );
   std::cout << *matrix_ptr << std::endl;
 
   /***
    * Set the right-hand side vector.
    */
   Vector x( size, 1.0 );
   Vector b( size );
   matrix_ptr->vectorProduct( x, b );
   x = 0.0;
   std::cout << "Vector b = " << b << std::endl;
 
   /***
    * Setup solver of the linear system.
    */
   using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
   LinearSolver solver;
   solver.setMatrix( matrix_ptr );
   solver.setOmega( 0.0005 );
 
   /***
    * Setup monitor of the iterative solver.
    */
   using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
   IterativeSolverMonitorType monitor;
   TNL::Solvers::SolverMonitorThread monitorThread( monitor );
   monitor.setRefreshRate( 10 );  // refresh rate in milliseconds
   monitor.setVerbose( 1 );
   monitor.setStage( "Jacobi stage:" );
   solver.setSolverMonitor( monitor );
   solver.setConvergenceResidue( 1.0e-6 );
   solver.solve( b, x );
   monitor.stopMainLoop();
   std::cout << "Vector x = " << x << std::endl;
}
 
int
main( int argc, char* argv[] )
{
   std::cout << "Solving linear system on host: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Sequential >();
 
#ifdef __CUDACC__
   std::cout << "Solving linear system on CUDA device: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Cuda >();
#endif
}

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:   15987 RES:    0.061998
  Jacobi stage: ITER:   38391 RES:   0.0019854
  Jacobi stage: ITER:   54651 RES:  0.00016337
  Jacobi stage: ITER:   75287 RES:  6.8639e-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:    4481 RES:     0.36909
  Jacobi stage: ITER:    9113 RES:     0.17828
  Jacobi stage: ITER:   13487 RES:    0.091025
  Jacobi stage: ITER:   17889 RES:    0.046277
  Jacobi stage: ITER:   20843 RES:    0.029406
  Jacobi stage: ITER:   22867 RES:    0.021549
  Jacobi stage: ITER:   27019 RES:    0.011388
  Jacobi stage: ITER:   30363 RES:   0.0068136
  Jacobi stage: ITER:   32649 RES:   0.0047945
  Jacobi stage: ITER:   37107 RES:   0.0024182
  Jacobi stage: ITER:   41563 RES:   0.0012197
  Jacobi stage: ITER:   45969 RES:  0.00061971
  Jacobi stage: ITER:   50447 RES:   0.0003116
  Jacobi stage: ITER:   54871 RES:  0.00015794
  Jacobi stage: ITER:   55735 RES:  0.00013831
  Jacobi stage: ITER:   57135 RES:  0.00011155
  Jacobi stage: ITER:   59075 RES:  8.2803e-05
  Jacobi stage: ITER:   60795 RES:  6.3578e-05
  Jacobi stage: ITER:   61291 RES:  5.8914e-05
  Jacobi stage: ITER:   63599 RES:  4.1329e-05
  Jacobi stage: ITER:   66387 RES:  2.6932e-05
  Jacobi stage: ITER:   68423 RES:  1.9699e-05
  Jacobi stage: ITER:   69113 RES:  1.7713e-05
  Jacobi stage: ITER:   70263 RES:  1.4849e-05
  Jacobi stage: ITER:   71015 RES:   1.323e-05
  Jacobi stage: ITER:   71955 RES:  1.1451e-05
  Jacobi stage: ITER:   73947 RES:  8.4325e-06
  Jacobi stage: ITER:   74871 RES:  7.3168e-06
  Jacobi stage: ITER:   75747 RES:  6.3956e-06
  Jacobi stage: ITER:   76719 RES:  5.5086e-06
  Jacobi stage: ITER:   77715 RES:  4.7272e-06
  Jacobi stage: ITER:   78611 RES:  4.1194e-06
  Jacobi stage: ITER:   80419 RES:  3.1205e-06
  Jacobi stage: ITER:   82299 RES:  2.3378e-06
  Jacobi stage: ITER:   83003 RES:  2.0982e-06
  Jacobi stage: ITER:   84197 RES:  1.7461e-06
  Jacobi stage: ITER:   85135 RES:  1.5123e-06
  Jacobi stage: ITER:   85887 RES:  1.3473e-06
  Jacobi stage: ITER:   86951 RES:  1.1441e-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:

#include <iostream>
#include <memory>
#include <chrono>
#include <thread>
#include <TNL/Timer.h>
#include <TNL/Matrices/SparseMatrix.h>
#include <TNL/Devices/Sequential.h>
#include <TNL/Devices/Cuda.h>
#include <TNL/Solvers/Linear/Jacobi.h>
 
template< typename Device >
void
iterativeLinearSolverExample()
{
   /***
    * Set the following matrix (dots represent zero matrix elements):
    *
    *   /  2.5 -1    .    .    .   \
    *   | -1    2.5 -1    .    .   |
    *   |  .   -1    2.5 -1.   .   |
    *   |  .    .   -1    2.5 -1   |
    *   \  .    .    .   -1    2.5 /
    */
   using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
   using Vector = TNL::Containers::Vector< double, Device >;
   const int size( 5 );
   auto matrix_ptr = std::make_shared< MatrixType >();
   matrix_ptr->setDimensions( size, size );
   matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
 
   auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
   {
      const int rowIdx = row.getRowIndex();
      if( rowIdx == 0 ) {
         row.setElement( 0, rowIdx, 2.5 );     // diagonal element
         row.setElement( 1, rowIdx + 1, -1 );  // element above the diagonal
      }
      else if( rowIdx == size - 1 ) {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
      }
      else {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
         row.setElement( 2, rowIdx + 1, -1.0 );  // element above the diagonal
      }
   };
 
   /***
    * Set the matrix elements.
    */
   matrix_ptr->forAllRows( f );
   std::cout << *matrix_ptr << std::endl;
 
   /***
    * Set the right-hand side vector.
    */
   Vector x( size, 1.0 );
   Vector b( size );
   matrix_ptr->vectorProduct( x, b );
   x = 0.0;
   std::cout << "Vector b = " << b << std::endl;
 
   /***
    * Setup solver of the linear system.
    */
   using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >;
   LinearSolver solver;
   solver.setMatrix( matrix_ptr );
   solver.setOmega( 0.0005 );
 
   /***
    * Setup monitor of the iterative solver.
    */
   using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >;
   IterativeSolverMonitorType monitor;
   TNL::Solvers::SolverMonitorThread mmonitorThread( monitor );
   monitor.setRefreshRate( 10 );  // refresh rate in milliseconds
   monitor.setVerbose( 1 );
   monitor.setStage( "Jacobi stage:" );
   TNL::Timer timer;
   monitor.setTimer( timer );
   timer.start();
   solver.setSolverMonitor( monitor );
   solver.setConvergenceResidue( 1.0e-6 );
   solver.solve( b, x );
   monitor.stopMainLoop();
   std::cout << "Vector x = " << x << std::endl;
}
 
int
main( int argc, char* argv[] )
{
   std::cout << "Solving linear system on host: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Sequential >();
 
#ifdef __CUDACC__
   std::cout << "Solving linear system on CUDA device: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Cuda >();
#endif
}

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:0.001821
 ELA: 0.10194  Jacobi stage: ITER:   10105 RES:     0.15303
 ELA: 0.20204  Jacobi stage: ITER:   20607 RES:    0.030473
 ELA: 0.30215  Jacobi stage: ITER:   32975 RES:   0.0045618
 ELA: 0.40225  Jacobi stage: ITER:   53981 RES:  0.00018102
 ELA: 0.50235  Jacobi stage: ITER:   75075 RES:  7.0911e-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: 0.01128
 ELA:  0.1114  Jacobi stage: ITER:    2675 RES:     0.50774
 ELA:  0.2115  Jacobi stage: ITER:    4419 RES:     0.37302
 ELA: 0.31161  Jacobi stage: ITER:    6645 RES:     0.26126
 ELA: 0.41172  Jacobi stage: ITER:    7509 RES:      0.2284
 ELA: 0.51182  Jacobi stage: ITER:    8989 RES:     0.18171
 ELA: 0.61192  Jacobi stage: ITER:   10035 RES:     0.15473
 ELA: 0.71202  Jacobi stage: ITER:   10455 RES:     0.14505
 ELA: 0.81213  Jacobi stage: ITER:   11223 RES:      0.1289
 ELA: 0.91223  Jacobi stage: ITER:   12875 RES:    0.099998
 ELA:  1.0123  Jacobi stage: ITER:   14695 RES:    0.075608
 ELA:  1.1124  Jacobi stage: ITER:   15723 RES:    0.064564
 ELA:  1.2125  Jacobi stage: ITER:   16639 RES:    0.056089
 ELA:  1.3126  Jacobi stage: ITER:   17433 RES:    0.049634
 ELA:  1.4127  Jacobi stage: ITER:   18371 RES:    0.042987
 ELA:  1.5128  Jacobi stage: ITER:   19617 RES:    0.035489
 ELA:  1.6129  Jacobi stage: ITER:   20777 RES:    0.029697
 ELA:   1.713  Jacobi stage: ITER:   21827 RES:    0.025281
 ELA:  1.8131  Jacobi stage: ITER:   23019 RES:    0.021051
 ELA:  1.9132  Jacobi stage: ITER:   23851 RES:    0.018526
 ELA:  2.0133  Jacobi stage: ITER:   26261 RES:     0.01279
 ELA:  2.1134  Jacobi stage: ITER:   27637 RES:    0.010354
 ELA:  2.2135  Jacobi stage: ITER:   28595 RES:   0.0089395
 ELA:  2.3136  Jacobi stage: ITER:   30341 RES:   0.0068345
 ELA:  2.4137  Jacobi stage: ITER:   31239 RES:   0.0059558
 ELA:  2.5138  Jacobi stage: ITER:   32177 RES:    0.005155
 ELA:  2.6155  Jacobi stage: ITER:   35009 RES:   0.0033367
 ELA:  2.7156  Jacobi stage: ITER:   39851 RES:   0.0015865
 ELA:  2.8157  Jacobi stage: ITER:   45127 RES:  0.00070549
 ELA:  2.9158  Jacobi stage: ITER:   50559 RES:  0.00030629
 ELA:  3.0159  Jacobi stage: ITER:   54849 RES:  0.00015842
 ELA:   3.116  Jacobi stage: ITER:   59831 RES:  7.3725e-05
 ELA:  3.2161  Jacobi stage: ITER:   64083 RES:  3.8368e-05
 ELA:  3.3162  Jacobi stage: ITER:   68629 RES:   1.908e-05
 ELA:  3.4163  Jacobi stage: ITER:   72835 RES:  1.0003e-05
 ELA:  3.5164  Jacobi stage: ITER:   77689 RES:  4.7446e-06
 ELA:  3.6222  Jacobi stage: ITER:   82715 RES:  2.1931e-06
 ELA:Vector x = [   3.7223  Jacobi stage: ITER:0.999999        ,  87829 RES:  9.9949e-07
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.

#include <iostream>
#include <memory>
#include <TNL/Matrices/SparseMatrix.h>
#include <TNL/Devices/Host.h>
#include <TNL/Devices/Cuda.h>
#include <TNL/Solvers/Linear/TFQMR.h>
#include <TNL/Solvers/Linear/Preconditioners/Diagonal.h>
 
template< typename Device >
void
iterativeLinearSolverExample()
{
   /***
    * Set the following matrix (dots represent zero matrix elements):
    *
    *   /  2.5 -1    .    .    .   \
    *   | -1    2.5 -1    .    .   |
    *   |  .   -1    2.5 -1.   .   |
    *   |  .    .   -1    2.5 -1   |
    *   \  .    .    .   -1    2.5 /
    */
   using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
   using Vector = TNL::Containers::Vector< double, Device >;
   const int size( 5 );
   auto matrix_ptr = std::make_shared< MatrixType >();
   matrix_ptr->setDimensions( size, size );
   matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
 
   auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
   {
      const int rowIdx = row.getRowIndex();
      if( rowIdx == 0 ) {
         row.setElement( 0, rowIdx, 2.5 );     // diagonal element
         row.setElement( 1, rowIdx + 1, -1 );  // element above the diagonal
      }
      else if( rowIdx == size - 1 ) {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
      }
      else {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
         row.setElement( 2, rowIdx + 1, -1.0 );  // element above the diagonal
      }
   };
 
   /***
    * Set the matrix elements.
    */
   matrix_ptr->forAllRows( f );
   std::cout << *matrix_ptr << std::endl;
 
   /***
    * Set the right-hand side vector.
    */
   Vector x( size, 1.0 );
   Vector b( size );
   matrix_ptr->vectorProduct( x, b );
   x = 0.0;
   std::cout << "Vector b = " << b << std::endl;
 
   /***
    * Solve the linear system using diagonal (Jacobi) preconditioner.
    */
   using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >;
   using Preconditioner = TNL::Solvers::Linear::Preconditioners::Diagonal< MatrixType >;
   auto preconditioner_ptr = std::make_shared< Preconditioner >();
   preconditioner_ptr->update( matrix_ptr );
   LinearSolver solver;
   solver.setMatrix( matrix_ptr );
   solver.setPreconditioner( preconditioner_ptr );
   solver.setConvergenceResidue( 1.0e-6 );
   solver.solve( b, x );
   std::cout << "Vector x = " << x << std::endl;
}
 
int
main( int argc, char* argv[] )
{
   std::cout << "Solving linear system on host: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Sequential >();
 
#ifdef __CUDACC__
   std::cout << "Solving linear system on CUDA device: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Cuda >();
#endif
}

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:

#include <iostream>
#include <memory>
#include <TNL/Matrices/SparseMatrix.h>
#include <TNL/Devices/Host.h>
#include <TNL/Devices/Cuda.h>
#include <TNL/Solvers/LinearSolverTypeResolver.h>
 
template< typename Device >
void
iterativeLinearSolverExample()
{
   /***
    * Set the following matrix (dots represent zero matrix elements):
    *
    *   /  2.5 -1    .    .    .   \
    *   | -1    2.5 -1    .    .   |
    *   |  .   -1    2.5 -1.   .   |
    *   |  .    .   -1    2.5 -1   |
    *   \  .    .    .   -1    2.5 /
    */
   using MatrixType = TNL::Matrices::SparseMatrix< double, Device >;
   using Vector = TNL::Containers::Vector< double, Device >;
   const int size( 5 );
   auto matrix_ptr = std::make_shared< MatrixType >();
   matrix_ptr->setDimensions( size, size );
   matrix_ptr->setRowCapacities( Vector( { 2, 3, 3, 3, 2 } ) );
 
   auto f = [ = ] __cuda_callable__( typename MatrixType::RowView & row ) mutable
   {
      const int rowIdx = row.getRowIndex();
      if( rowIdx == 0 ) {
         row.setElement( 0, rowIdx, 2.5 );     // diagonal element
         row.setElement( 1, rowIdx + 1, -1 );  // element above the diagonal
      }
      else if( rowIdx == size - 1 ) {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
      }
      else {
         row.setElement( 0, rowIdx - 1, -1.0 );  // element below the diagonal
         row.setElement( 1, rowIdx, 2.5 );       // diagonal element
         row.setElement( 2, rowIdx + 1, -1.0 );  // element above the diagonal
      }
   };
 
   /***
    * Set the matrix elements.
    */
   matrix_ptr->forAllRows( f );
   std::cout << *matrix_ptr << std::endl;
 
   /***
    * Set the right-hand side vector.
    */
   Vector x( size, 1.0 );
   Vector b( size );
   matrix_ptr->vectorProduct( x, b );
   x = 0.0;
   std::cout << "Vector b = " << b << std::endl;
 
   /***
    * Solve the linear system using diagonal (Jacobi) preconditioner.
    */
   auto solver_ptr = TNL::Solvers::getLinearSolver< MatrixType >( "tfqmr" );
   auto preconditioner_ptr = TNL::Solvers::getPreconditioner< MatrixType >( "diagonal" );
   preconditioner_ptr->update( matrix_ptr );
   solver_ptr->setMatrix( matrix_ptr );
   solver_ptr->setPreconditioner( preconditioner_ptr );
   solver_ptr->setConvergenceResidue( 1.0e-6 );
   solver_ptr->solve( b, x );
   std::cout << "Vector x = " << x << std::endl;
}
 
int
main( int argc, char* argv[] )
{
   std::cout << "Solving linear system on host: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Sequential >();
 
#ifdef __CUDACC__
   std::cout << "Solving linear system on CUDA device: " << std::endl;
   iterativeLinearSolverExample< TNL::Devices::Cuda >();
#endif
}

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 ]

Table of Contents

Introduction

Iterative solvers of linear systems

Basic setup

Setup with a solver monitor

Setup with preconditioner

Choosing the solver and preconditioner type at runtime