Writting documentation on linear preconditioners.

Documentation/Examples/Solvers/Linear/CMakeLists.txt

@@ -2,6 +2,7 @@ set( COMMON_EXAMPLES
    IterativeLinearSolverExample
    IterativeLinearSolverWithMonitorExample
    IterativeLinearSolverWithTimerExample
+   IterativeLinearSolverWithPreconditionerExample
 )
 
 if( BUILD_CUDA )

Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp

+#include <iostream>
+#include <memory>
+#include <TNL/Algorithms/ParallelFor.h>
+#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.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 HAVE_CUDA
+   std::cout << "Solving linear system on CUDA device: " << std::endl;
+   iterativeLinearSolverExample< TNL::Devices::Cuda >();
+#endif
+}

Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu

+IterativeLinearSolverWithPreconditionerExample.cpp
\ No newline at end of file

Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md

@@ -79,4 +79,14 @@ The result looks as follows:
 
 \include IterativeLinearSolverWithTimerExample.out
 
-## Setup with preconditioner
+## Setup with preconditioner
\ No newline at end of file
+
+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 (\ref TNL::Solvers::Linear::Jacobi and \ref TNL::Solvers::Linear::SOR). The following example shows how to setup an iterative solver of linear systems with preconditioning.
+
+\includelineno Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp
+
+In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen on the lines (68-72). Similar to the matrix of the linear system, the preconditioner is passed to the solver by the means of  smart shared pointer (\ref std::shared_ptr). The instance is created on the lines 68 and 69. Next we just need to connect the solver with the preconditioner (line 72, \ref TNL::Solvers::Linear::LinearSolver).
+
+The result looks as follows:
+
+\include IterativeLinearSolverWithPreconditionerExample.out

src/TNL/Solvers/IterativeSolverMonitor.h

@@ -20,6 +20,24 @@ namespace TNL {
  *
  * \tparam Real is a type of the floating-point arithmetics.
  * \tparam Index is an indexing type.
+ *
+ * The following example shows how to use the iterative solver monitor for monitoring
+ * convergence of linear iterative solver:
+ *
+ * \includelineno Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp
+ *
+ * The result looks as follows:
+ *
+ * \include IterativeLinearSolverWithMonitorExample.out
+ *
+ * The following example shows how to employ timer (\ref TNL::Timer) to the monitor
+ * of iterative solvers:
+ *
+ * \includelineno Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp
+ *
+ * The result looks as follows:
+ *
+ * \include IterativeLinearSolverWithTimerExample.out
  */
 template< typename Real = double,
           typename Index = int >

src/TNL/Solvers/Linear/LinearSolver.h

@@ -34,6 +34,16 @@ namespace TNL {
  * \ref LinearSolver::setPreconditioner.
  *
  * \tparam Matrix is type of matrix representing the linear system.
+ *
+ * The following example demonstrates the use iterative linear solvers:
+ *
+ * \includelineno Solvers/Linear/IterativeLinearSolverExample.cpp
+ *
+ * The result looks as follows:
+ *
+ * \include IterativeLinearSolverExample.out
+ *
+ * See also \ref TNL::Solvers::IterativeSolverMonitor for monitoring of iterative solvers.
  */
 template< typename Matrix >
 class LinearSolver

src/TNL/Solvers/Linear/Preconditioners/ILU0.h

@@ -76,6 +76,7 @@ class ILU0_impl< Matrix, Real, Devices::Host, Index >
 : public Preconditioner< Matrix >
 {
    public:
+
       /**
        * \brief Floating point type used for computations.
        */

src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h

@@ -27,6 +27,15 @@ namespace TNL {
  * \brief Base class for preconditioners of of iterative solvers of linear systems.
  *
  * \tparam Matrix is type of matrix describing the linear system.
+ *
+ * The following example shows how to setup an iterative solver of linear systems with
+ * preconditioning:
+ *
+ * \includelineno Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp
+ *
+ * The result looks as follows:
+ *
+ * \include IterativeLinearSolverWithPreconditionerExample.out
  */
 template< typename Matrix >
 class Preconditioner