From 00103113c71bc68bc0e53c014f95dfb0c02d4300 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Thu, 30 Dec 2021 13:09:22 +0100
Subject: [PATCH] Added example on getLinearSolver and getPreconditioner.

---
 .../Examples/Solvers/Linear/CMakeLists.txt    |  1 +
 ...iveLinearSolverWithRuntimeTypesExample.cpp | 84 +++++++++++++++++++
 ...tiveLinearSolverWithRuntimeTypesExample.cu |  1 +
 3 files changed, 86 insertions(+)
 create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp
 create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu

diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt
index 0d2354617c..97c77e7317 100644
--- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt
+++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt
@@ -3,6 +3,7 @@ set( COMMON_EXAMPLES
    IterativeLinearSolverWithMonitorExample
    IterativeLinearSolverWithTimerExample
    IterativeLinearSolverWithPreconditionerExample
+   IterativeLinearSolverWithRuntimeTypesExample
 )
 
 if( BUILD_CUDA )
diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp
new file mode 100644
index 0000000000..cbe7707cfe
--- /dev/null
+++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp
@@ -0,0 +1,84 @@
+#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/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->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/IterativeLinearSolverWithRuntimeTypesExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu
new file mode 120000
index 0000000000..ab40ffd370
--- /dev/null
+++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu
@@ -0,0 +1 @@
+IterativeLinearSolverWithRuntimeTypesExample.cpp
\ No newline at end of file
-- 
GitLab