From 8b1748bdff1b6317a76ab4b71acac862e0229b6c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Wed, 29 Dec 2021 21:08:38 +0100
Subject: [PATCH] Writting documentation on linear preconditioners.
---
.../Examples/Solvers/Linear/CMakeLists.txt | 1 +
...eLinearSolverWithPreconditionerExample.cpp | 87 +++++++++++++++++++
...veLinearSolverWithPreconditionerExample.cu | 1 +
.../Solvers/Linear/tutorial_Linear_solvers.md | 12 ++-
src/TNL/Solvers/IterativeSolverMonitor.h | 18 ++++
src/TNL/Solvers/Linear/LinearSolver.h | 10 +++
src/TNL/Solvers/Linear/Preconditioners/ILU0.h | 1 +
.../Linear/Preconditioners/Preconditioner.h | 9 ++
8 files changed, 138 insertions(+), 1 deletion(-)
create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp
create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu
diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt
index a07ba232a4..0d2354617c 100644
--- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt
+++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt
@@ -2,6 +2,7 @@ set( COMMON_EXAMPLES
IterativeLinearSolverExample
IterativeLinearSolverWithMonitorExample
IterativeLinearSolverWithTimerExample
+ IterativeLinearSolverWithPreconditionerExample
)
if( BUILD_CUDA )
diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp
new file mode 100644
index 0000000000..3a21410744
--- /dev/null
+++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp
@@ -0,0 +1,87 @@
+#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
+}
diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu
new file mode 120000
index 0000000000..dc6a12a588
--- /dev/null
+++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu
@@ -0,0 +1 @@
+IterativeLinearSolverWithPreconditionerExample.cpp
\ No newline at end of file
diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md
index c58b28273f..25e709ae66 100644
--- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md
+++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md
@@ -79,4 +79,14 @@ The result looks as follows:
\include IterativeLinearSolverWithTimerExample.out
-## Setup with preconditioner
\ No newline at end of file
+## 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 (\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
diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h
index 651f38cded..8f42a6b28d 100644
--- a/src/TNL/Solvers/IterativeSolverMonitor.h
+++ b/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 >
diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h
index 858f081890..f2de9a8f4a 100644
--- a/src/TNL/Solvers/Linear/LinearSolver.h
+++ b/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
diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h
index 4b35864484..6ce3fef515 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h
+++ b/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.
*/
diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h
index 43b429da74..e0dc721e36 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h
+++ b/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
--
GitLab