diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 7aa7364299a853bc261ebd6fba0b225204218121..b27168ca0ac6630480f6882a4210280e25781506 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -2,6 +2,7 @@ ADD_SUBDIRECTORY( Algorithms ) ADD_SUBDIRECTORY( Containers ) ADD_SUBDIRECTORY( Pointers ) ADD_SUBDIRECTORY( Matrices ) +ADD_SUBDIRECTORY( Solvers ) set( COMMON_EXAMPLES ) diff --git a/Documentation/Examples/Solvers/CMakeLists.txt b/Documentation/Examples/Solvers/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..79b9712aa0e935b0a8d27ef31e91c9ec260885a7 --- /dev/null +++ b/Documentation/Examples/Solvers/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory( Linear ) diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..97c77e73175b0225d68f98b868e73b53e3e0c1d5 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -0,0 +1,28 @@ +set( COMMON_EXAMPLES + IterativeLinearSolverExample + IterativeLinearSolverWithMonitorExample + IterativeLinearSolverWithTimerExample + IterativeLinearSolverWithPreconditionerExample + IterativeLinearSolverWithRuntimeTypesExample +) + +if( BUILD_CUDA ) + foreach( target IN ITEMS ${COMMON_EXAMPLES} ) + cuda_add_executable( ${target}-cuda ${target}.cu OPTIONS ) + add_custom_command( COMMAND ${target}-cuda > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out OUTPUT ${target}.out ) + set( CUDA_OUTPUTS ${CUDA_OUTPUTS} ${target}.out ) + endforeach() +else() + foreach( target IN ITEMS ${COMMON_EXAMPLES} ) + add_executable( ${target} ${target}.cpp ) + add_custom_command( COMMAND ${target} > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out OUTPUT ${target}.out ) + set( HOST_OUTPUTS ${HOST_OUTPUTS} ${target}.out ) + endforeach() +endif() + +IF( BUILD_CUDA ) + ADD_CUSTOM_TARGET( RunLinearSolversExamples-cuda ALL DEPENDS ${CUDA_OUTPUTS} ) +ELSE() + ADD_CUSTOM_TARGET( RunLinearSolversExamples ALL DEPENDS ${HOST_OUTPUTS} ) +ENDIF() + diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp new file mode 100644 index 0000000000000000000000000000000000000000..222a8d86616cbffe00b36df7ea809f78a6f34906 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp @@ -0,0 +1,82 @@ +#include +#include +#include +#include +#include +#include +#include + +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.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/IterativeLinearSolverExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cu new file mode 120000 index 0000000000000000000000000000000000000000..17884b9012a6daec0d453ad0f99165e34583b079 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverExample.cpp \ No newline at end of file diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp new file mode 100644 index 0000000000000000000000000000000000000000..92cf9c840276d9384b0ce19c25ad5aca9cc06674 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp @@ -0,0 +1,97 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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.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::Host >(); + +#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/IterativeLinearSolverWithMonitorExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cu new file mode 120000 index 0000000000000000000000000000000000000000..c0474305b0e999447fb99a4205e2265263c8f8d7 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithMonitorExample.cpp \ No newline at end of file diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp new file mode 100644 index 0000000000000000000000000000000000000000..035dbb6ce2fb8014166c4f201c031f278d5dda56 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp @@ -0,0 +1,87 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +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 0000000000000000000000000000000000000000..dc6a12a588680a497eed939cf34de678a3f1f87b --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithPreconditionerExample.cpp \ No newline at end of file diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cbe7707cfe494d990d48528c7672b97e36dbaa63 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp @@ -0,0 +1,84 @@ +#include +#include +#include +#include +#include +#include +#include + +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 0000000000000000000000000000000000000000..ab40ffd37040b07ed10b5dc8a4d10c57a88796c3 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithRuntimeTypesExample.cpp \ No newline at end of file diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a9ae53c1d5dca47e46d348c9595c94af8780992d --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp @@ -0,0 +1,101 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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.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 HAVE_CUDA + std::cout << "Solving linear system on CUDA device: " << std::endl; + iterativeLinearSolverExample< TNL::Devices::Cuda >(); +#endif +} diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cu new file mode 120000 index 0000000000000000000000000000000000000000..4c497154d238d93a1ca7ff5c151bb51472e378e4 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithTimerExample.cpp \ No newline at end of file diff --git a/Documentation/Tutorials/CMakeLists.txt b/Documentation/Tutorials/CMakeLists.txt index 05ed1f33cc0ea1fd0c653439b9cd40f50a34510c..e30aab1e24035f047b4bff709089a46c0070c89c 100644 --- a/Documentation/Tutorials/CMakeLists.txt +++ b/Documentation/Tutorials/CMakeLists.txt @@ -5,3 +5,4 @@ add_subdirectory( ReductionAndScan ) add_subdirectory( Pointers ) add_subdirectory( Matrices ) add_subdirectory( Meshes ) +add_subdirectory( Solvers ) diff --git a/Documentation/Tutorials/Solvers/CMakeLists.txt b/Documentation/Tutorials/Solvers/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..79b9712aa0e935b0a8d27ef31e91c9ec260885a7 --- /dev/null +++ b/Documentation/Tutorials/Solvers/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory( Linear ) diff --git a/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt b/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md new file mode 100644 index 0000000000000000000000000000000000000000..b79dbbe1b766954f63ecadd354b8cfd4d0cb6b1d --- /dev/null +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -0,0 +1,105 @@ +\page tutorial_Linear_solvers Linear solvers tutorial + +[TOC] + +# Introduction + +Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the followiing iterative methods: + +1. Stationary methods + 1. [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method) (\ref TNL::Solvers::Linear::Jacobi) + 2. [Successive-overrelaxation method, SOR]([https://en.wikipedia.org/wiki/Successive_over-relaxation]) (\ref TNL::Solvers::Linear::SOR) +2. Krylov subspace methods + 1. [Conjugate gradient method, CG](https://en.wikipedia.org/wiki/Conjugate_gradient_method) (\ref TNL::Solvers::Linear::CG) + 2. [Biconjugate gradient stabilized method, BICGStab](https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method) (\ref TNL::Solvers::Linear::BICGStab) + 3. [Biconjugate gradient stabilized method, BICGStab(l)](https://dspace.library.uu.nl/bitstream/handle/1874/16827/sleijpen_93_bicgstab.pdf) (\ref TNL::Solvers::Linear::BICGStabL) + 4. [Transpose-free quasi-minimal residual method, TFQMR]([https://second.wiki/wiki/algoritmo_tfqmr]) (\ref TNL::Solvers::Linear::TFQMR) + 5. [Generalized minimal residual method, GMRES](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method) (\ref 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 \ref TNL::Solvers::Linear::Jacobi and \ref TNL::Solevrs::Linear::SOR) can be combined with the following preconditioners: + +1. [Diagonal or Jacobi](http://netlib.org/linalg/html_templates/node55.html)) - \ref TNL::Solvers::Linear::Preconditioners::Diagonal +2. ILU (Incomplete LU) - CPU only currently + 1. [ILU(0)](https://en.wikipedia.org/wiki/Incomplete_LU_factorization) \ref TNL::Solvers::Linear::Preconditioners::ILU0 + 2. [ILUT (ILU with thresholding)](https://www-users.cse.umn.edu/~saad/PDF/umsi-92-38.pdf) \ref TNL::Solvers::Linear::Preconditioners::ILUT + +# Iterative solvers of linear systems + +## Basic setup + +All iterative solvers for linear systems can be found in the namespace \ref TNL::Solvers::Linear. The following example shows the use the iterative solvers: + +\includelineno Solvers/Linear/IterativeLinearSolverExample.cpp + +In this example we solve a linear system \f$ A \vec x = \vec b \f$ where + +\f[ +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) +\f] + +The right-hand side vector \f$\vec b \f$ is set to \f$( 1.5, 0.5, 0.5, 0.5, 1.5 )^T \f$ so that the exact solution is \f$ \vec x = ( 1, 1, 1, 1, 1 )^T\f$. The matrix elements of \f$A $\f$ is set on the lines 12-51 by the means of the method \ref TNL::Matrices::SparseMatrix::forAllElements. In this example, we use the sparse matrix but any other matrix type can be used as well (see the namespace \ref TNL::Matrices). Next we set the solution vector \f$ \vec x = ( 1, 1, 1, 1, 1 )^T\f$ (line 57) and multiply it with matrix \f$ A \f$ to get the right-hand side vector \f$\vec b\f$ (lines 58-59). Finally, we reset the vector \f$\vec x \f$ to zero vector. + +To solve the linear system, we use TFQMR method (line 66), as an example. Other solvers can be used as well (see the namespace \ref TNL::Solvers::Linear). The solver needs only one template parameter which is the matrix type. Next we create an instance of the solver (line 67 ) and set the matrix of the linear system (line 68). Note, that matrix is passed to the solver as a shared smart pointer (\ref std::shared_ptr). This is why we created an instance of the smart pointer on the line 24 instead of the sparse matrix itself. The solver is executed on the line 69 by calling the method \ref TNL::Solvers::Linear::LinearSolver::solve. The method accepts the right-hand side vector \f$ \vec b\f$ and the solution vector \f$ \vec x\f$. + +The result looks as follows: + +\include IterativeLinearSolverExample.out + +## 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 of the solver status in general. For this purpose, TNL offers solver monitors. The solver monitor prints (or somehow visualizes) the number of iterations, the residue of the current solution approximation or some other metrics. Sometimes such information is printed after each iteration or after every ten iterations. The problem of this approach is the fact that one iteration of the solver may take only few milliseconds but also several minutes. In the former case, the monitor creates overwhelming amount of output which may even slowdown the solver. In the later case, the user waits long time for update of the solver status. The monitor in TNL rather runs in separate thread and it refreshes the status of the solver in preset time periods. The use of the iterative solver monitor is demonstrated in the following example. + +\includelineno Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp + +On the lines 1-70, we setup 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. On the line 71, we set the relaxation parameter \f$ \omega \f$ of the Jacobi solver to 0.0005 (\ref TNL::Solvers::Linear::Jacobi). The reason is to slowdown the convergence because we want to see some iterations in this example. Next we create an instance of the solver monitor (lines 76 and 77) and we create a special thread for the monitor (line 78, \ref TNL::Solvers::SolverMonitorThread ). We set the refresh rate of the monitor to 10 milliseconds (line 79, \ref TNL::Solvers::SolverMonitor::setRefreshRate). We set a verbosity of the monitor to 1 (line 80 \ref TNL::Solvers::IterativeSolverMonitor::setVerbosity ). Next we set a name of the solver stage (line 81, \ref TNL::Solvers::IterativeSolverMonitor::setStage). The monitor stages serve for distinguishing between different phases or stages of more complex solvers (for example when the linear system solver is embedded into a time dependent PDE solver). Next we connect the solver with the monitor (line 82, \ref TNL::Solvers::IterativeSolver::setSolverMonitor). Finally we start the solver (line 83, \ref TNL::Solvers::Linear::Jacobi::start) and when the solver finishes we have to stop the monitor (line 84, \ref TNL::Solvers::SolverMonitor::stopMainLoop). + +The result looks as follows: + +\include IterativeLinearSolverWithMonitorExample.out + +The monitoring of the solver can be improved by time elapsed since the beginning of the computation as demonstrated in the following example: + +\includelineno Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp + +The only changes happen on lines 83-85 where we create an instance of TNL timer (line 83, \ref TNL::Timer), connect it with the monitor (line 84, \ref TNL::Solvers::SolverMonitor::setTimer) and start the timer (line 85, \ref TNL::Timer::start). + +The result looks as follows: + +\include IterativeLinearSolverWithTimerExample.out + +## 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 have to initialize the preconditioner (line 70, \ref TNL::Solvers::Linear::Preconditioners::Preconditioner::update). The method `update` 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 (line 73, \ref TNL::Solvers::Linear::LinearSolver). + +The result looks as follows: + +\include IterativeLinearSolverWithPreconditionerExample.out + +## 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 offers choosing the type of both linear solver and preconditioner at runtime by means of functions \ref TNL::Solvers::getLinearSolver and \ref TNL::Solvers::getPreconditioner. The following example shows how to use these functions: + +\includelineno Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp + +We still stay with the same problem and the only changes can be seen on lines 66-70. We first create an instance of shared pointer holding the solver (line 66, \ref TNL::Solvers::getLinearSolver) and the same with the preconditioner (line 67, \ref TNL::Solvers::getPreconditioner). 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: + +\include IterativeLinearSolverWithRuntimeTypesExample.out \ No newline at end of file diff --git a/Documentation/Tutorials/index.md b/Documentation/Tutorials/index.md index 031de3faee9cd64c347950dec75983fd71a27a36..900641f1bf42d81181f31c72128265c3545075b2 100644 --- a/Documentation/Tutorials/index.md +++ b/Documentation/Tutorials/index.md @@ -10,5 +10,6 @@ 6. [Sorting](tutorial_Sorting.html) 7. [Cross-device pointers](tutorial_Pointers.html) 8. [Matrices](tutorial_Matrices.html) -9. [Segments aka sparse formats](tutorial_Segments.html) -10. [Unstructured meshes](tutorial_Meshes.html) +9. [Linear solvers](tutorial_Linear_solvers.html) +10. [Segments aka sparse formats](tutorial_Segments.html) +11. [Unstructured meshes](tutorial_Meshes.html) diff --git a/src/Benchmarks/ODESolvers/Euler.hpp b/src/Benchmarks/ODESolvers/Euler.hpp index 840c0a45aab49d7412417219b1ac97b39286fbaa..72fb4867a856ff21108c6990c690acffdff53598 100644 --- a/src/Benchmarks/ODESolvers/Euler.hpp +++ b/src/Benchmarks/ODESolvers/Euler.hpp @@ -81,8 +81,6 @@ bool Euler< Problem, SolverMonitor >::solve( DofVectorPointer& u ) this->resetIterations(); this->setResidue( this->getConvergenceResidue() + 1.0 ); - this -> refreshSolverMonitor(); - /**** * Start the main loop */ @@ -128,17 +126,12 @@ bool Euler< Problem, SolverMonitor >::solve( DofVectorPointer& u ) currentTau = this -> getStopTime() - time; //we don't want to keep such tau else this -> tau = currentTau; - this->refreshSolverMonitor(); - /**** * Check stop conditions. */ if( time >= this->getStopTime() || ( this -> getConvergenceResidue() != 0.0 && this->getResidue() < this -> getConvergenceResidue() ) ) - { - this -> refreshSolverMonitor(); return true; - } if( this -> cflCondition != 0.0 ) { diff --git a/src/Benchmarks/ODESolvers/Merson.hpp b/src/Benchmarks/ODESolvers/Merson.hpp index 3bb6f9c96f772062d1f9c18c6fd35a406322ce3f..2fa53ed92893f2c926e20e525dfb4b2616bd1289 100644 --- a/src/Benchmarks/ODESolvers/Merson.hpp +++ b/src/Benchmarks/ODESolvers/Merson.hpp @@ -157,8 +157,6 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& u ) this->resetIterations(); this->setResidue( this->getConvergenceResidue() + 1.0 ); - this->refreshSolverMonitor(); - /**** * Start the main loop */ @@ -196,7 +194,6 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& u ) if( ! this->nextIteration() ) return false; } - this->refreshSolverMonitor(); /**** * Compute the new time step. @@ -221,12 +218,8 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& u ) //cerr << "this->getConvergenceResidue() = " << this->getConvergenceResidue() << std::endl; if( time >= this->getStopTime() || ( this->getConvergenceResidue() != 0.0 && this->getResidue() < this->getConvergenceResidue() ) ) - { - this->refreshSolverMonitor( true ); return true; - } } - this->refreshSolverMonitor( true ); return this->checkConvergence(); }; diff --git a/src/Python/pytnl/tnl/SparseMatrix.h b/src/Python/pytnl/tnl/SparseMatrix.h index aa0ea33941f0c1179b3f482b1a2ee4e1a4537cf0..2c66c6a98f40ed2eff61f5c6cd0ccaac02ff9fa5 100644 --- a/src/Python/pytnl/tnl/SparseMatrix.h +++ b/src/Python/pytnl/tnl/SparseMatrix.h @@ -149,7 +149,6 @@ void export_Matrix( py::module & m, const char* name ) // TODO: these two don't work //.def("addMatrix", &Matrix::addMatrix) //.def("getTransposition", &Matrix::getTransposition) - .def("performSORIteration", &Matrix::template performSORIteration< VectorType, VectorType >) // TODO: export for more types .def("assign", []( Matrix& matrix, const Matrix& other ) -> Matrix& { return matrix = other; diff --git a/src/TNL/Algorithms/AtomicOperations.h b/src/TNL/Algorithms/AtomicOperations.h index dbc1d18a7d16a7e39fe475ea2880da34f16c25fe..4833e33ce81ab5931156a6421e2037e16bc0d251 100644 --- a/src/TNL/Algorithms/AtomicOperations.h +++ b/src/TNL/Algorithms/AtomicOperations.h @@ -41,6 +41,21 @@ struct AtomicOperations< Devices::Host > } }; +template<> +struct AtomicOperations< Devices::Sequential > +{ + // this is __cuda_callable__ only to silence nvcc warnings (all methods inside class + // template specializations must have the same execution space specifier, otherwise + // nvcc complains) + TNL_NVCC_HD_WARNING_DISABLE + template< typename Value > + __cuda_callable__ + static void add( Value& v, const Value& a ) + { + v += a; + } +}; + template<> struct AtomicOperations< Devices::Cuda > { diff --git a/src/TNL/Atomic.h b/src/TNL/Atomic.h index 0d78a660ff95bbaea1acc21f337e34bfaaca78c4..ff8d419a1be1a3b04b2c76580dd5d442e0afff6e 100644 --- a/src/TNL/Atomic.h +++ b/src/TNL/Atomic.h @@ -15,6 +15,7 @@ #include // std::atomic #include +#include #include // double-precision atomicAdd function for Maxwell and older GPUs @@ -96,6 +97,18 @@ public: } }; +template< typename T > +class Atomic< T, Devices::Sequential > : public Atomic< T, Devices::Host > +{ + using Base = Atomic< T, Devices::Host >; + public: + + using Base::Atomic; + using Base::operator=; + using Base::fetch_max; + using Base::fetch_min; +}; + template< typename T > class Atomic< T, Devices::Cuda > { diff --git a/src/TNL/File.h b/src/TNL/File.h index cef110e1633537ecb5b13770cc805dccb1e2786f..9048c022e58a8e590345cfdfccf7538696435ba3 100644 --- a/src/TNL/File.h +++ b/src/TNL/File.h @@ -58,9 +58,9 @@ class File * Opens file with given \e fileName in some \e mode from \ref std::ios_base::openmode. * Note that the file is always opened in binary mode, i.e. \ref std::ios_base::binary * is always added to \e mode. - * + * * Throws \ref std::ios_base::failure on failure. - * + * * \param fileName String which indicates file name. * \param mode Indicates in what mode the file will be opened - see \ref std::ios_base::openmode. */ @@ -69,7 +69,7 @@ class File /** * \brief Closes the file. - * + * * Throws \ref std::ios_base::failure on failure. */ void close(); @@ -133,7 +133,7 @@ class File * \tparam Index type of index by which the elements are indexed. * \param buffer buffer that is going to be saved to the file. * \param elements number of elements saved to the file. - * + * * See \ref File::load for examples. */ template< typename Type, typename TargetType = Type, typename Allocator = Allocators::Host< Type > > diff --git a/src/TNL/Matrices/DenseMatrix.h b/src/TNL/Matrices/DenseMatrix.h index 3d0e2d62da85d102c502a5770cbe3fc2be013408..d46b844e66d55812eed984d0d5c99fa082a3acd2 100644 --- a/src/TNL/Matrices/DenseMatrix.h +++ b/src/TNL/Matrices/DenseMatrix.h @@ -906,12 +906,6 @@ class DenseMatrix : public Matrix< Real, Device, Index, RealAllocator > void getTransposition( const Matrix& matrix, const RealType& matrixMultiplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment operator with exactly the same type of the dense matrix. * diff --git a/src/TNL/Matrices/DenseMatrix.hpp b/src/TNL/Matrices/DenseMatrix.hpp index a42421aa7a661fba0ec7fde144a335322af6f157..ae3f7ee6bb7a7621a06f55ef93eb9f0a833f7c6a 100644 --- a/src/TNL/Matrices/DenseMatrix.hpp +++ b/src/TNL/Matrices/DenseMatrix.hpp @@ -1037,28 +1037,6 @@ void DenseMatrix< Real, Device, Index, Organization, RealAllocator >::getTranspo } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization, - typename RealAllocator > - template< typename Vector1, typename Vector2 > -void DenseMatrix< Real, Device, Index, Organization, RealAllocator >::performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ), diagonalValue; - for( IndexType i = 0; i < this->getColumns(); i++ ) - { - if( i == row ) - diagonalValue = this->getElement( row, row ); - else - sum += this->getElement( row, i ) * x[ i ]; - } - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / diagonalValue * ( b[ row ] - sum ); -} - template< typename Real, typename Device, typename Index, diff --git a/src/TNL/Matrices/DenseMatrixView.h b/src/TNL/Matrices/DenseMatrixView.h index 89a2219b3bddcc9b9139d8984694010fa298b61b..712a565ae255a9997bf0e3fc113e566dc56b2941 100644 --- a/src/TNL/Matrices/DenseMatrixView.h +++ b/src/TNL/Matrices/DenseMatrixView.h @@ -856,12 +856,6 @@ class DenseMatrixView : public MatrixView< Real, Device, Index > void getTransposition( const Matrix& matrix, const RealType& matrixMultiplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment operator with DenseMatrix. * diff --git a/src/TNL/Matrices/DenseMatrixView.hpp b/src/TNL/Matrices/DenseMatrixView.hpp index 3a44269d1d78923c5b0e40091c9b20d9eb303483..336fcbb9aa8ac3a5ddfae3f9528058783a3ac769 100644 --- a/src/TNL/Matrices/DenseMatrixView.hpp +++ b/src/TNL/Matrices/DenseMatrixView.hpp @@ -1006,28 +1006,6 @@ void DenseMatrixView< Real, Device, Index, Organization >::getTransposition( con } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization > - template< typename Vector1, typename Vector2 > -void DenseMatrixView< Real, Device, Index, Organization >::performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ), diagonalValue; - for( IndexType i = 0; i < this->getColumns(); i++ ) - { - if( i == row ) - diagonalValue = this->getElement( row, row ); - else - sum += this->getElement( row, i ) * x[ i ]; - } - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / diagonalValue * ( b[ row ] - sum ); -} - - template< typename Real, typename Device, typename Index, diff --git a/src/TNL/Matrices/DistributedMatrix.h b/src/TNL/Matrices/DistributedMatrix.h index 250521f8e1e22d04353e3937f5ad18af6de4432a..cdc6387ccb3147a710a6e956b30b6a49da3abc42 100644 --- a/src/TNL/Matrices/DistributedMatrix.h +++ b/src/TNL/Matrices/DistributedMatrix.h @@ -118,13 +118,6 @@ public: vectorProduct( const InVector& inVector, OutVector& outVector ) const; - // FIXME: does not work for distributed matrices, here only due to common interface - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - protected: LocalRangeType localRowRange; IndexType rows = 0; // global rows count diff --git a/src/TNL/Matrices/DistributedMatrix_impl.h b/src/TNL/Matrices/DistributedMatrix_impl.h index 498868509513834bc52d910e0cc861d49a8e92e8..e66e86ad054a17a215bdc247210e471596f51bda 100644 --- a/src/TNL/Matrices/DistributedMatrix_impl.h +++ b/src/TNL/Matrices/DistributedMatrix_impl.h @@ -266,17 +266,5 @@ vectorProduct( const InVector& inVector, // outVector.startSynchronization(); } -template< typename Matrix > - template< typename Vector1, typename Vector2 > -void -DistributedMatrix< Matrix >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - getLocalMatrix().performSORIteration( b, row, x, omega ); -} - } // namespace Matrices } // namespace TNL diff --git a/src/TNL/Matrices/LambdaMatrix.h b/src/TNL/Matrices/LambdaMatrix.h index 2b788ed5270540a58582153b8c37e87dd346c2a4..0278a00ef40e83a651dcfa4f902db4545326ebde 100644 --- a/src/TNL/Matrices/LambdaMatrix.h +++ b/src/TNL/Matrices/LambdaMatrix.h @@ -490,13 +490,6 @@ class LambdaMatrix const IndexType begin = 0, IndexType end = 0 ) const; - - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Method for printing the matrix to output stream. * diff --git a/src/TNL/Matrices/LambdaMatrix.hpp b/src/TNL/Matrices/LambdaMatrix.hpp index 730c8b5b1d092d03ce1914b31722fdd5c88966dc..e348a6dba6229e6c23c38e81d87fe4c834794606 100644 --- a/src/TNL/Matrices/LambdaMatrix.hpp +++ b/src/TNL/Matrices/LambdaMatrix.hpp @@ -401,22 +401,6 @@ sequentialForAllRows( Function&& function ) const sequentialForRows( (IndexType) 0, this->getRows(), function ); } -template< typename MatrixElementsLambda, - typename CompressedRowLengthsLambda, - typename Real, - typename Device, - typename Index > - template< typename Vector1, typename Vector2 > -void -LambdaMatrix< MatrixElementsLambda, CompressedRowLengthsLambda, Real, Device, Index >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - throw Exceptions::NotImplementedError("performSORIteration is not implemented for SparseMatrixView"); -} - template< typename MatrixElementsLambda, typename CompressedRowLengthsLambda, typename Real, diff --git a/src/TNL/Matrices/MultidiagonalMatrix.h b/src/TNL/Matrices/MultidiagonalMatrix.h index 3a74b0b44ff0d0b5529e87416b36790113a02f6f..d56f1fd0619ae582f90e0557cbf30068cc5aeaaf 100644 --- a/src/TNL/Matrices/MultidiagonalMatrix.h +++ b/src/TNL/Matrices/MultidiagonalMatrix.h @@ -1053,13 +1053,6 @@ class MultidiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator > void getTransposition( const MultidiagonalMatrix< Real2, Device, Index2 >& matrix, const RealType& matrixMultiplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - __cuda_callable__ - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of exactly the same matrix type. * diff --git a/src/TNL/Matrices/MultidiagonalMatrix.hpp b/src/TNL/Matrices/MultidiagonalMatrix.hpp index 244188831cf05b7d2dc5218060ff35c589e0d016..6ece15b4ca76b060d2e6d027ebe8efba2c252d8a 100644 --- a/src/TNL/Matrices/MultidiagonalMatrix.hpp +++ b/src/TNL/Matrices/MultidiagonalMatrix.hpp @@ -808,29 +808,6 @@ getTransposition( const MultidiagonalMatrix< Real2, Device, Index2 >& matrix, } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization, - typename RealAllocator, - typename IndexAllocator > - template< typename Vector1, typename Vector2 > -__cuda_callable__ -void MultidiagonalMatrix< Real, Device, Index, Organization, RealAllocator, IndexAllocator >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ); - if( row > 0 ) - sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ]; - if( row < this->getColumns() - 1 ) - sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ]; - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum ); -} - - // copy assignment template< typename Real, typename Device, diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.h b/src/TNL/Matrices/MultidiagonalMatrixView.h index 869ddb9732e8541fecf9847c1ed648f189012335..a0e13874d03f2072eda90dd19dd65a38067d2367 100644 --- a/src/TNL/Matrices/MultidiagonalMatrixView.h +++ b/src/TNL/Matrices/MultidiagonalMatrixView.h @@ -811,13 +811,6 @@ class MultidiagonalMatrixView : public MatrixView< Real, Device, Index > void getTransposition( const MultidiagonalMatrixView< Real2, Device, Index2 >& matrix, const RealType& matrixMultiplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - __cuda_callable__ - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of exactly the same matrix type. * diff --git a/src/TNL/Matrices/MultidiagonalMatrixView.hpp b/src/TNL/Matrices/MultidiagonalMatrixView.hpp index 550158e4d458261e65d89fb5bd5d4018f4eb9a32..65293b809e671abc4a57ce73b99741d90aefefa7 100644 --- a/src/TNL/Matrices/MultidiagonalMatrixView.hpp +++ b/src/TNL/Matrices/MultidiagonalMatrixView.hpp @@ -774,28 +774,6 @@ getTransposition( const MultidiagonalMatrixView< Real2, Device, Index2 >& matrix } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization > - template< typename Vector1, typename Vector2 > -__cuda_callable__ -void -MultidiagonalMatrixView< Real, Device, Index, Organization >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ); - if( row > 0 ) - sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ]; - if( row < this->getColumns() - 1 ) - sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ]; - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum ); -} - - template< typename Real, typename Device, typename Index, diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h index fa3441b01874b7f5c78a528fda7236dbc6f03838..76d47b50dc1db36bee622ab09b19e7032fe3fb34 100644 --- a/src/TNL/Matrices/SparseMatrix.h +++ b/src/TNL/Matrices/SparseMatrix.h @@ -1043,12 +1043,6 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator > const RealType& matrixMultiplicator = 1.0 ); */ - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of exactly the same matrix type. * diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp index 93cd9b173f786fdf7403cb7ef06ac6c753d6afb7..b7df383b0aecaf1ee1867f3253de4c5dfc55ce87 100644 --- a/src/TNL/Matrices/SparseMatrix.hpp +++ b/src/TNL/Matrices/SparseMatrix.hpp @@ -815,25 +815,6 @@ getTransposition( const SparseMatrix< Real2, Device, Index2 >& matrix, }*/ -template< typename Real, - typename Device, - typename Index, - typename MatrixType, - template< typename, typename, typename > class Segments, - typename ComputeReal, - typename RealAllocator, - typename IndexAllocator > -template< typename Vector1, typename Vector2 > -void -SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAllocator, IndexAllocator >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - this->view.performSORIteration( b, row, x, omega ); -} - // copy assignment template< typename Real, typename Device, diff --git a/src/TNL/Matrices/SparseMatrixView.h b/src/TNL/Matrices/SparseMatrixView.h index eaf692ac66d58c1778caa9847f36c65094cd7a4d..704d8c7ed98ccf3c10a0de9ddd1adf5b1b8dad83 100644 --- a/src/TNL/Matrices/SparseMatrixView.h +++ b/src/TNL/Matrices/SparseMatrixView.h @@ -821,12 +821,6 @@ class SparseMatrixView : public MatrixView< Real, Device, Index > const IndexType begin = 0, IndexType end = 0 ) const; - template< typename Vector1, typename Vector2 > - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of any matrix type. * . diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp index 2af26efd5986fd268f04ada7e1074479e7178de4..3d0ce734274506c98e249358f81b3f36e018d4a8 100644 --- a/src/TNL/Matrices/SparseMatrixView.hpp +++ b/src/TNL/Matrices/SparseMatrixView.hpp @@ -601,7 +601,7 @@ forElements( IndexType begin, IndexType end, Function& function ) const if( localIdx < columns ) { if( isBinary() ) - function( rowIdx, localIdx, columns_view[ globalIdx ], 1 ); + function( rowIdx, localIdx, columns_view[ globalIdx ], ( RealType ) 1.0 ); else function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] ); } @@ -827,23 +827,6 @@ getTransposition( const SparseMatrixView< Real2, Device, Index2 >& matrix, }*/ -template< typename Real, - typename Device, - typename Index, - typename MatrixType, - template< typename, typename > class SegmentsView, - typename ComputeReal > -template< typename Vector1, typename Vector2 > -void -SparseMatrixView< Real, Device, Index, MatrixType, SegmentsView, ComputeReal >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - throw Exceptions::NotImplementedError("performSORIteration is not implemented for SparseMatrixView"); -} - template< typename Real, typename Device, typename Index, diff --git a/src/TNL/Matrices/TridiagonalMatrix.h b/src/TNL/Matrices/TridiagonalMatrix.h index 45e6d132fecd289df3ba1ccf608eebf39a1342d7..bb47f68126cef4172f43688583989a047f8685ac 100644 --- a/src/TNL/Matrices/TridiagonalMatrix.h +++ b/src/TNL/Matrices/TridiagonalMatrix.h @@ -945,13 +945,6 @@ class TridiagonalMatrix : public Matrix< Real, Device, Index, RealAllocator > void getTransposition( const TridiagonalMatrix< Real2, Device, Index2 >& matrix, const RealType& matrixTriplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - __cuda_callable__ - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of exactly the same matrix type. * diff --git a/src/TNL/Matrices/TridiagonalMatrix.hpp b/src/TNL/Matrices/TridiagonalMatrix.hpp index 93ebcd8b39e68ba0f70079497a8b3cb64d5e5e43..80236b7de7ca023ea4100f58d0bb7204add17b1c 100644 --- a/src/TNL/Matrices/TridiagonalMatrix.hpp +++ b/src/TNL/Matrices/TridiagonalMatrix.hpp @@ -658,27 +658,6 @@ void TridiagonalMatrix< Real, Device, Index, Organization, RealAllocator >::getT } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization, - typename RealAllocator > - template< typename Vector1, typename Vector2 > -__cuda_callable__ -void TridiagonalMatrix< Real, Device, Index, Organization, RealAllocator >::performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ); - if( row > 0 ) - sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ]; - if( row < this->getColumns() - 1 ) - sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ]; - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum ); -} - - // copy assignment template< typename Real, typename Device, diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h index c8e0ecdcad2e6f9b4fc45f2b56e09a0070866c5a..3420f68154b8192b4e33270535890d3d50dc3665 100644 --- a/src/TNL/Matrices/TridiagonalMatrixView.h +++ b/src/TNL/Matrices/TridiagonalMatrixView.h @@ -779,13 +779,6 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index > void getTransposition( const TridiagonalMatrixView< Real2, Device, Index2 >& matrix, const RealType& matrixMultiplicator = 1.0 ); - template< typename Vector1, typename Vector2 > - __cuda_callable__ - void performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega = 1.0 ) const; - /** * \brief Assignment of exactly the same matrix type. * diff --git a/src/TNL/Matrices/TridiagonalMatrixView.hpp b/src/TNL/Matrices/TridiagonalMatrixView.hpp index c67073381047ef8a282bae91fcb5e470f5a484d1..64da406244ea7bd44d71c4e025b7a10a356f9627 100644 --- a/src/TNL/Matrices/TridiagonalMatrixView.hpp +++ b/src/TNL/Matrices/TridiagonalMatrixView.hpp @@ -757,28 +757,6 @@ getTransposition( const TridiagonalMatrixView< Real2, Device, Index2 >& matrix, } } -template< typename Real, - typename Device, - typename Index, - ElementsOrganization Organization > - template< typename Vector1, typename Vector2 > -__cuda_callable__ -void -TridiagonalMatrixView< Real, Device, Index, Organization >:: -performSORIteration( const Vector1& b, - const IndexType row, - Vector2& x, - const RealType& omega ) const -{ - RealType sum( 0.0 ); - if( row > 0 ) - sum += this->getElementFast( row, row - 1 ) * x[ row - 1 ]; - if( row < this->getColumns() - 1 ) - sum += this->getElementFast( row, row + 1 ) * x[ row + 1 ]; - x[ row ] = ( 1.0 - omega ) * x[ row ] + omega / this->getElementFast( row, row ) * ( b[ row ] - sum ); -} - - template< typename Real, typename Device, typename Index, diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index 97d3f1c565b418ec93ae2a408e0656e4a5168a1f..a0e1ccd0de55520166922286bb3dd97d549ff001 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -17,84 +17,230 @@ #include namespace TNL { -namespace Solvers { - + namespace Solvers { + +/** + * \brief Base class for iterative solvers. + * + * \tparam Real is a floating point type used for computations. + * \tparam Index is an indexing type. + * \tparam IterativeSolverMonitor< Real, Index > is type of an object used for monitoring of the convergence. + */ template< typename Real, typename Index, typename SolverMonitor = IterativeSolverMonitor< Real, Index > > class IterativeSolver { -public: - using SolverMonitorType = SolverMonitor; - - IterativeSolver() = default; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ); - - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ); - - void setMaxIterations( const Index& maxIterations ); - - const Index& getMaxIterations() const; - - void setMinIterations( const Index& minIterations ); - - const Index& getMinIterations() const; - - const Index& getIterations() const; - - void setConvergenceResidue( const Real& convergenceResidue ); - - const Real& getConvergenceResidue() const; - - void setDivergenceResidue( const Real& divergenceResidue ); - - const Real& getDivergenceResidue() const; - - void setResidue( const Real& residue ); - - const Real& getResidue() const; - - void setRefreshRate( const Index& refreshRate ); - - void setSolverMonitor( SolverMonitorType& solverMonitor ); - - void resetIterations(); - - bool nextIteration(); - - bool checkNextIteration(); - - bool checkConvergence(); - - void refreshSolverMonitor( bool force = false ); - -protected: - Index maxIterations = 1000000000; - - Index minIterations = 0; - - Index currentIteration = 0; - - Real convergenceResidue = 1e-6; - - // If the current residue is greater than divergenceResidue, the solver is stopped. - Real divergenceResidue = std::numeric_limits< float >::max(); - - Real currentResidue = 0; - - SolverMonitor* solverMonitor = nullptr; - - Index refreshRate = 1; - - String residualHistoryFileName = ""; - - std::ofstream residualHistoryFile; + public: + + /** + * \brief Type of an object used for monitoring of the convergence. + */ + using SolverMonitorType = SolverMonitor; + + /** + * \brief Default constructor. + */ + IterativeSolver() = default; + + /** + * \brief This method defines configuration entries for setup of the iterative solver. + * + * The following entries are defined: + * + * \e max-iterations - maximal number of iterations the solver \b may perform. + * + * \e min-iterations - minimal number of iterations the solver \b must perform. + * + * \e convergence-residue - convergence occurs when the residue drops \b bellow this limit. + * + * \e divergence-residue - divergence occurs when the residue \b exceeds given limit. + * + * \e refresh-rate - number of milliseconds between solver monitor refreshes. + * + * \e residual-history-file - path to the file where the residual history will be saved. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ); + + /** + * \brief Method for setup of the iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ); + + /** + * \brief Sets the maximal number of iterations the solver is \b allowed to perform. + * + * If the number of iterations performed by the solver exceeds given limit, the divergence occurs. + * + * \param maxIterations maximal number of allowed iterations. + */ + void setMaxIterations( const Index& maxIterations ); + + /** + * \brief Gets the maximal number of iterations the solver is \b allowed to perform. + * + * See \ref IterativeSolver::setMaxIterations. + * + * \return maximal number of allowed iterations. + */ + const Index& getMaxIterations() const; + + /** + * \brief Sets the minimal number of iterations the solver is \b supposed to do. + * + * \param minIterations minimal number of iterations the solver is supposed to do. + */ + void setMinIterations( const Index& minIterations ); + + /** + * \brief Gets the minimal number of iterations the solver is \b supposed to do. + * + * \return minimal number of iterations the solver is supposed to do. + */ + const Index& getMinIterations() const; + + /** + * \brief Gets the number of iterations performed by the solver so far. + * + * \return number of iterations performed so far. + */ + const Index& getIterations() const; + + /** + * \brief Sets the threshold for the convergence. + * + * The convergence occurs when the residue drops \b bellow this limit. + * + * \param convergenceResidue is threshold for the convergence. + */ + void setConvergenceResidue( const Real& convergenceResidue ); + + /** + * \brief Gets the the convergence threshold. + * + * See \ref IterativeSolver::setConvergenceResidue. + * + * \return the convergence threshold. + */ + const Real& getConvergenceResidue() const; + + /** + * \brief Sets the residue limit for the divergence criterion. + * + * The divergence occurs when the residue \b exceeds the limit. + * + * \param divergenceResidue the residue limit of the divergence. + */ + void setDivergenceResidue( const Real& divergenceResidue ); + + /** + * \brief Gets the limit for the divergence criterion. + * + * See \ref IterativeSolver::setDivergenceResidue. + * + * \return the residue limit fo the divergence. + */ + const Real& getDivergenceResidue() const; + + /** + * \brief Sets the residue reached at the current iteration. + * + * \param residue reached at the current iteration. + */ + void setResidue( const Real& residue ); + + /** + * \brief Gets the residue reached at the current iteration. + * + * \return residue reached at the current iteration. + */ + const Real& getResidue() const; + + /** + * \brief Sets the refresh rate (in milliseconds) for the solver monitor. + * + * \param refreshRate of the solver monitor in milliseconds. + */ + void setRefreshRate( const Index& refreshRate ); + + /** + * \brief Sets the solver monitor object. + * + * The solver monitor is an object for monitoring the status of the iterative solver. + * Usually it prints the number of iterations, current residue or elapsed time. + * + * \param solverMonitor is an object for monitoring the iterative solver. + */ + void setSolverMonitor( SolverMonitorType& solverMonitor ); + + /** + * \brief Sets the the number of the current iterations to zero. + */ + void resetIterations(); + + /** + * \brief Proceeds to the next iteration. + * + * \return \e true if the solver is allowed to do the next iteration. + * \return \e false if the solver is \b not allowed to do the next iteration. This may + * happen because the divergence occurred. + */ + bool nextIteration(); + + /** + * \brief Checks if the solver is allowed to the next iteration. + * + * \return true \e true if the solver is allowed to do the next iteration. + * \return \e false if the solver is \b not allowed to do the next iteration. This may + * happen because the divergence occurred. + */ + bool checkNextIteration(); + + /** + * \brief Checks whether the convergence occurred already. + * + * \return \e true if the convergence already occured. + * \return \e false if the convergence did not occur yet. + */ + bool checkConvergence(); + + /** + * \brief Refreshes the solver monitor. + */ + //void refreshSolverMonitor(); + + protected: + Index maxIterations = 1000000000; + + Index minIterations = 0; + + Index currentIteration = 0; + + Real convergenceResidue = 1e-6; + + // If the current residue is greater than divergenceResidue, the solver is stopped. + Real divergenceResidue = std::numeric_limits< Real >::max(); + + Real currentResidue = 0; + + SolverMonitor* solverMonitor = nullptr; + + Index refreshRate = 1; + + String residualHistoryFileName = ""; + + std::ofstream residualHistoryFile; }; -} // namespace Solvers + } // namespace Solvers } // namespace TNL #include diff --git a/src/TNL/Solvers/IterativeSolver.hpp b/src/TNL/Solvers/IterativeSolver.hpp index f894f37119e247247e6a7c0e265e7a341acc6063..0a78dff89e982ce1bcc244cc5f7df33e111a1cf7 100644 --- a/src/TNL/Solvers/IterativeSolver.hpp +++ b/src/TNL/Solvers/IterativeSolver.hpp @@ -127,8 +127,6 @@ bool IterativeSolver< Real, Index, SolverMonitor >:: checkNextIteration() { - this->refreshSolverMonitor(); - if( std::isnan( this->getResidue() ) || this->getIterations() > this->getMaxIterations() || ( this->getResidue() > this->getDivergenceResidue() && this->getIterations() >= this->getMinIterations() ) || @@ -236,6 +234,7 @@ IterativeSolver< Real, Index, SolverMonitor >:: setRefreshRate( const Index& refreshRate ) { this->refreshRate = refreshRate; + this->solverMonitor->setRefreshRate( this->refreshRate ); } template< typename Real, typename Index, typename SolverMonitor > @@ -246,18 +245,18 @@ setSolverMonitor( SolverMonitorType& solverMonitor ) this->solverMonitor = &solverMonitor; } -template< typename Real, typename Index, typename SolverMonitor > +/*template< typename Real, typename Index, typename SolverMonitor > void IterativeSolver< Real, Index, SolverMonitor >:: -refreshSolverMonitor( bool force ) +refreshSolverMonitor() { if( this->solverMonitor ) { this->solverMonitor->setIterations( this->getIterations() ); this->solverMonitor->setResidue( this->getResidue() ); - this->solverMonitor->setRefreshRate( this-> refreshRate ); + this->solverMonitor->setRefreshRate( this->refreshRate ); } -} +}*/ } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h index 07dee8de99a0a151e40bffbbae082b81469069e7..8f42a6b28d7f8ac2a4bdfcf0e925e94ad09c60ac 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.h +++ b/src/TNL/Solvers/IterativeSolverMonitor.h @@ -13,51 +13,133 @@ #include namespace TNL { -namespace Solvers { - -template< typename Real, typename Index> + namespace Solvers { + +/** + * \brief Object for monitoring convergence of iterative solvers. + * + * \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 > class IterativeSolverMonitor : public SolverMonitor { -public: - typedef Index IndexType; - typedef Real RealType; - - IterativeSolverMonitor(); - - void setStage( const std::string& stage ); - - void setTime( const RealType& time ); - - void setTimeStep( const RealType& timeStep ); - - void setIterations( const IndexType& iterations ); - - void setResidue( const RealType& residue ); - - void setVerbose( const IndexType& verbose ); - - void setNodesPerIteration( const IndexType& nodes ); - - virtual void refresh(); - -protected: - int getLineWidth(); - - std::string stage, saved_stage; - - std::atomic_bool saved, attributes_changed; - - RealType time, saved_time, timeStep, saved_timeStep, residue, saved_residue, elapsed_time_before_refresh, last_mlups; - //TODO: Move MLUPS to LBM solver only i.e create solver monitor for LBM - - IndexType iterations, saved_iterations, iterations_before_refresh; - - IndexType verbose; - - IndexType nodesPerIteration; + public: + + /** + * \brief A type of the floating-point arithmetics. + */ + using RealType = Real; + + /** + * \brief A type for indexing. + */ + using IndexType = Index; + + /** + * \brief Construct with no parameters. + */ + IterativeSolverMonitor(); + + /** + * \brief This method can be used for naming a stage of the monitored solver. + * + * The stage name can be used to differ between various stages of iterative solvers. + * + * \param stage is name of the solver stage. + */ + void setStage( const std::string& stage ); + + /** + * \brief Set the time of the simulated evolution if it is time dependent. + * + * This can be used for example when solving parabolic or hyperbolic PDEs. + * + * \param time time of the simulated evolution. + */ + void setTime( const RealType& time ); + + /** + * \brief Set the time step for time dependent iterative solvers. + * + * \param timeStep time step of the time dependent iterative solver. + */ + void setTimeStep( const RealType& timeStep ); + + /** + * \brief Set number of the current iteration. + * + * \param iterations is number of the current iteration. + */ + void setIterations( const IndexType& iterations ); + + /** + * \brief Set residue of the current approximation of the solution. + * + * \param residue is a residue of the current approximation of the solution. + */ + void setResidue( const RealType& residue ); + + /** + * \brief Set up the verbosity of the monitor. + * + * \param verbose is the new value of the verbosity of the monitor. + */ + void setVerbose( const IndexType& verbose ); + + /** + * \brief Set the number of nodes of the numerical mesh or lattice. + * + * This can be used to compute the number of nodes processed per one second. + * + * \param nodes is number of nodes of the numerical mesh or lattice. + */ + void setNodesPerIteration( const IndexType& nodes ); + + /** + * \brief Causes that the monitor prints out the status of the solver. + */ + virtual void refresh(); + + protected: + + int getLineWidth(); + + std::string stage, saved_stage; + + std::atomic_bool saved, attributes_changed; + + RealType time, saved_time, timeStep, saved_timeStep, residue, saved_residue, elapsed_time_before_refresh, last_mlups; + //TODO: Move MLUPS to LBM solver only i.e create solver monitor for LBM + + IndexType iterations, saved_iterations, iterations_before_refresh; + + // TODO: move verbose to SolverMonitor + IndexType verbose; + + IndexType nodesPerIteration; }; -} // namespace Solvers + } // namespace Solvers } // namespace TNL #include diff --git a/src/TNL/Solvers/IterativeSolverMonitor.hpp b/src/TNL/Solvers/IterativeSolverMonitor.hpp index a1fa8c45434ee43828ec74ccc7fff823796ba79e..996f47f959d93274aad6d1f6806b0d50f947bf7c 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.hpp +++ b/src/TNL/Solvers/IterativeSolverMonitor.hpp @@ -152,13 +152,19 @@ void IterativeSolverMonitor< Real, Index > :: refresh() // verbose == 1, attributes were not updated since the last refresh return; - print_item( " ELA:" ); - print_item( real_to_string( getElapsedTime(), 5 ), 8 ); - print_item( " T:" ); - print_item( real_to_string( (saved) ? saved_time : time, 5 ), 8 ); - if( (saved) ? saved_timeStep : timeStep > 0 ) { - print_item( " TAU:" ); - print_item( real_to_string( (saved) ? saved_timeStep : timeStep, 5 ), 10 ); + if( this->timer != nullptr ) + { + print_item( " ELA:" ); + print_item( real_to_string( getElapsedTime(), 5 ), 8 ); + } + if( this->timeStep > 0 ) + { + print_item( " T:" ); + print_item( real_to_string( (saved) ? saved_time : time, 5 ), 8 ); + if( (saved) ? saved_timeStep : timeStep > 0 ) { + print_item( " TAU:" ); + print_item( real_to_string( (saved) ? saved_timeStep : timeStep, 5 ), 10 ); + } } const std::string displayed_stage = (saved) ? saved_stage : stage; diff --git a/src/TNL/Solvers/Linear/BICGStab.h b/src/TNL/Solvers/Linear/BICGStab.h index 233590ed7aed38dd46927de100e9fb0c11a9a283..c9d973ad1284bd7d2c9f28cf5b3d609f944714c1 100644 --- a/src/TNL/Solvers/Linear/BICGStab.h +++ b/src/TNL/Solvers/Linear/BICGStab.h @@ -10,46 +10,109 @@ #pragma once -#include "LinearSolver.h" +#include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { +/** + * \brief Iterative solver of linear systems based on the biconjugate gradient stabilized (BICGStab) method. + * + * This solver can be used for nonsymmetric linear systems. + * + * See [Wikipedia](https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method) for more details. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class BICGStab : public LinearSolver< Matrix > { using Base = LinearSolver< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; + using VectorType = typename Traits< Matrix >::VectorType; - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ); + public: - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) override; + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; - bool solve( ConstVectorViewType b, VectorViewType x ) override; + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename Base::DeviceType; -protected: - void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; - void preconditioned_matvec( ConstVectorViewType src, VectorViewType dst ); + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; - void setSize( const VectorViewType& x ); + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; - bool exact_residue = false; + /** + * \brief This is method defines configuration entries for setup of the linear iterative solver. + * + * In addition to config entries defined by \ref IterativeSolver::configSetup, this method + * defines the following: + * + * \e bicgstab-exact-residue - says whether the BiCGstab should compute the exact residue in + * each step (true) or to use a cheap approximation (false). + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ); - typename Traits< Matrix >::VectorType r, r_ast, p, s, Ap, As, M_tmp; + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) override; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); + + void preconditioned_matvec( ConstVectorViewType src, VectorViewType dst ); + + void setSize( const VectorViewType& x ); + + bool exact_residue = false; + + VectorType r, r_ast, p, s, Ap, As, M_tmp; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "BICGStab.hpp" +#include diff --git a/src/TNL/Solvers/Linear/BICGStab.hpp b/src/TNL/Solvers/Linear/BICGStab.hpp index afbdbc39945aa0a07934aed80bc30d3d76760856..9cf51034f4630699217ed27763fc8cf66b137a63 100644 --- a/src/TNL/Solvers/Linear/BICGStab.hpp +++ b/src/TNL/Solvers/Linear/BICGStab.hpp @@ -119,7 +119,6 @@ solve( ConstVectorViewType b, VectorViewType x ) } } - this->refreshSolverMonitor( true ); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/BICGStabL.h b/src/TNL/Solvers/Linear/BICGStabL.h index f30a236342382406e6f40fa14824a813ab7761e5..44f4f2264fff25518f0f498ea2672234d65acea8 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.h +++ b/src/TNL/Solvers/Linear/BICGStabL.h @@ -9,6 +9,26 @@ /* See Copyright Notice in tnl/Copyright */ /* + * TODO: further variations to explore: + * + * [5] Gerard L. G. Sleijpen and Henk A. van der Vorst, "Reliable updated + * residuals in hybrid Bi-CG methods", Computing 56 (2), 141-163 (1996). + * [6] Gerard L. G. Sleijpen and Henk A. van der Vorst, "Maintaining convergence + * properties of BiCGstab methods in finite precision arithmetic", Numerical + * Algorithms 10, 203-223 (1995). + */ + +#pragma once + +#include + +namespace TNL { + namespace Solvers { + namespace Linear { + +/** + * \brief Iterative solver of linear systems based on the BICGStab(l) method. + * * BICGStabL implements an iterative solver for non-symmetric linear systems, * using the BiCGstab(l) algorithm described in [1] and [2]. It is a * generalization of the stabilized biconjugate-gradient (BiCGstab) algorithm @@ -16,37 +36,25 @@ * BiCGstab(2) is a slightly more efficient version of the BiCGstab2 algorithm * by Gutknecht [4], while BiCGstab(l>2) is a further generalization. * - * This code was implemented by: Jakub Klinkovsky - * * [1] Gerard L. G. Sleijpen and Diederik R. Fokkema, "BiCGstab(l) for linear * equations involving unsymmetric matrices with complex spectrum", * Electronic Trans. on Numerical Analysis 1, 11-32 (1993). + * * [2] Gerard L. G. Sleijpen, Henk A. van der Vorst, and Diederik R. Fokkema, * "BiCGstab(l) and other Hybrid Bi-CG Methods", Numerical Algorithms 7, * 75-109 (1994). + * * [3] Henk A. van der Vorst, "Bi-CGSTAB: A fast and smoothly converging variant * of Bi-CG for the solution of nonsymmetric linear systems, SIAM Journal on * scientific and Statistical Computing 13.2, 631-644 (1992). + * * [4] Martin H. Gutknecht, "Variants of BiCGStab for matrices with complex * spectrum", IPS Research Report No. 91-14 (1991). * - * TODO: further variations to explore: + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. * - * [5] Gerard L. G. Sleijpen and Henk A. van der Vorst, "Reliable updated - * residuals in hybrid Bi-CG methods", Computing 56 (2), 141-163 (1996). - * [6] Gerard L. G. Sleijpen and Henk A. van der Vorst, "Maintaining convergence - * properties of BiCGstab methods in finite precision arithmetic", Numerical - * Algorithms 10, 203-223 (1995). + * \tparam Matrix is type of matrix describing the linear system. */ - -#pragma once - -#include "LinearSolver.h" - -namespace TNL { -namespace Solvers { -namespace Linear { - template< typename Matrix > class BICGStabL : public LinearSolver< Matrix > @@ -56,50 +64,99 @@ class BICGStabL // compatibility shortcut using Traits = Linear::Traits< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - // distributed vectors/views - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; - using VectorType = typename Traits::VectorType; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ); - - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) override; - - bool solve( ConstVectorViewType b, VectorViewType x ) override; - -protected: - using DeviceVector = Containers::Vector< RealType, DeviceType, IndexType >; - using HostVector = Containers::Vector< RealType, Devices::Host, IndexType >; - - void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); - - void preconditioned_matvec( ConstVectorViewType src, VectorViewType dst ); - - void setSize( const VectorViewType& x ); - - int ell = 1; - - bool exact_residue = false; - - // matrices (in column-major format) - DeviceVector R, U; - // single vectors (distributed) - VectorType r_ast, M_tmp, res_tmp; - // host-only storage - HostVector T, sigma, g_0, g_1, g_2; - - IndexType size = 0; - IndexType ldSize = 0; + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + */ + using DeviceType = typename Base::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief This is method defines configuration entries for setup of the linear iterative solver. + * + * In addition to config entries defined by \ref IterativeSolver::configSetup, this method + * defines the following: + * + * \e bicgstab-ell - number of Bi-CG iterations before the MR part starts. + * + * \e bicgstab-exact-residue - says whether the BiCGstab should compute the exact residue in + * each step (true) or to use a cheap approximation (false). + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ); + + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) override; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + using VectorType = typename Traits::VectorType; + using DeviceVector = Containers::Vector< RealType, DeviceType, IndexType >; + using HostVector = Containers::Vector< RealType, Devices::Host, IndexType >; + + void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); + + void preconditioned_matvec( ConstVectorViewType src, VectorViewType dst ); + + void setSize( const VectorViewType& x ); + + int ell = 1; + + bool exact_residue = false; + + // matrices (in column-major format) + DeviceVector R, U; + // single vectors (distributed) + VectorType r_ast, M_tmp, res_tmp; + // host-only storage + HostVector T, sigma, g_0, g_1, g_2; + + IndexType size = 0; + IndexType ldSize = 0; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "BICGStabL.hpp" +#include diff --git a/src/TNL/Solvers/Linear/BICGStabL.hpp b/src/TNL/Solvers/Linear/BICGStabL.hpp index 7ae4999c9fad724a4c9be9f87285992b35f24d0e..b6814159016f5b74e91c7c5f8f91100abdf41e17 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.hpp +++ b/src/TNL/Solvers/Linear/BICGStabL.hpp @@ -238,8 +238,6 @@ solve( ConstVectorViewType b, VectorViewType x ) this->setResidue( sigma[ 0 ] / b_norm ); } } - - this->refreshSolverMonitor( true ); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/CG.h b/src/TNL/Solvers/Linear/CG.h index f2cad3f05de87e23a6144ffbe77b666dc25698a8..605c81fb3045ea70313c76771940b8d0ec43b00c 100644 --- a/src/TNL/Solvers/Linear/CG.h +++ b/src/TNL/Solvers/Linear/CG.h @@ -10,35 +10,78 @@ #pragma once -#include "LinearSolver.h" +#include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { +/** + * \brief Iterative solver of linear systems based on the conjugate gradient method. + * + * This solver can be used only for positive-definite linear systems. + * + * See [Wikipedia](https://en.wikipedia.org/wiki/Conjugate_gradient_method) for more details. + * + * \tparam Matrix is type of matrix describing the linear system. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + */ template< typename Matrix > class CG : public LinearSolver< Matrix > { using Base = LinearSolver< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; using VectorType = typename Traits< Matrix >::VectorType; - bool solve( ConstVectorViewType b, VectorViewType x ) override; + public: + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; -protected: - void setSize( const VectorViewType& x ); + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename Base::DeviceType; - VectorType r, p, Ap, z; + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + void setSize( const VectorViewType& x ); + + VectorType r, p, Ap, z; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "CG.hpp" +#include diff --git a/src/TNL/Solvers/Linear/CG.hpp b/src/TNL/Solvers/Linear/CG.hpp index 67442b4b1702ba408f3f9a18470772c11e2b361c..660e3a68d3639767dde7d1b085e393f2e0d185ed 100644 --- a/src/TNL/Solvers/Linear/CG.hpp +++ b/src/TNL/Solvers/Linear/CG.hpp @@ -107,7 +107,6 @@ solve( ConstVectorViewType b, VectorViewType x ) this->setResidue( std::sqrt(s1) / normb ); } - this->refreshSolverMonitor( true ); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/GMRES.h b/src/TNL/Solvers/Linear/GMRES.h index 7314a3a05b8f5c391311b4ba932c7421d1d7b7ee..67c5edc164ecd63bcc0485946cefd42ca30d7ee1 100644 --- a/src/TNL/Solvers/Linear/GMRES.h +++ b/src/TNL/Solvers/Linear/GMRES.h @@ -12,12 +12,37 @@ #pragma once -#include "LinearSolver.h" +#include namespace TNL { -namespace Solvers { -namespace Linear { - + namespace Solvers { + namespace Linear { + +/** + * \brief Iterative solver of linear systems based on the Generalized minimal residual (GMRES) method. + * + * This method can be used for solving of non-semmetric linear systems. This implementation + * offers various methods for the orthogonalization of vectors generated by the Arnoldi algorithm. + * The user may choose one of the following: + * + * \e CGS - Classical Gramm-Schmidt, + * + * \e CGSR - Classical Gramm-Schmidt with reorthogonalization, + * + * \e MGS - Modified Gramm-Schmidt, + * + * \e MGSR - Modified Gramm-Schmidt with reorthogonalization (the default one), + * + * \e CWY - Compact WY form of the Householder reflections). + * + * See [Wikipedia](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method) for more details. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class GMRES : public LinearSolver< Matrix > @@ -25,139 +50,195 @@ class GMRES using Base = LinearSolver< Matrix >; using Traits = Linear::Traits< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - // distributed vectors/views - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; - using VectorType = typename Traits::VectorType; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ); - - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) override; - - bool solve( ConstVectorViewType b, VectorViewType x ) override; - -protected: - // local vectors/views - using ConstDeviceView = typename Traits::ConstLocalViewType; - using DeviceView = typename Traits::LocalViewType; - using DeviceVector = typename Traits::LocalVectorType; - using HostView = typename DeviceView::template Self< RealType, Devices::Host >; - using HostVector = typename DeviceVector::template Self< RealType, Devices::Host >;; - - enum class Variant { CGS, CGSR, MGS, MGSR, CWY }; - -// nvcc allows __cuda_callable__ lambdas only in public methods -#ifdef __NVCC__ -public: -#endif - int orthogonalize_CGS( const int m, const RealType normb, const RealType beta ); -#ifdef __NVCC__ -protected: -#endif - - int orthogonalize_MGS( const int m, const RealType normb, const RealType beta ); - - int orthogonalize_CWY( const int m, const RealType normb, const RealType beta ); - - void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); - - void preconditioned_matvec( VectorViewType w, ConstVectorViewType v ); - -// nvcc allows __cuda_callable__ lambdas only in public methods -#ifdef __NVCC__ -public: -#endif - void hauseholder_generate( const int i, - VectorViewType y_i, - ConstVectorViewType z ); -#ifdef __NVCC__ -protected: -#endif - - void hauseholder_apply_trunc( HostView out, - const int i, + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + */ + using DeviceType = typename Base::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief This is method defines configuration entries for setup of the linear iterative solver. + * + * In addition to config entries defined by \ref IterativeSolver::configSetup, this method + * defines the following: + * + * \e gmres-variant - Algorithm used for the orthogonalization - CGS (Classical Gramm-Schmidt), + * CGSR(Classical Gramm-Schmidt with reorthogonalization), MGS(Modified Gramm-Schmidt), + * MGSR(Modified Gramm-Schmidt with reorthogonalization), CWY(Compact WY form of the Householder reflections). + * + * \e gmres-restarting-min - minimal number of iterations after which the GMRES restarts. + * + * \e gmres-restarting-max - maximal number of iterations after which the GMRES restarts. + * + * \e gmres-restarting-step-min - minimal adjusting step for the adaptivity of the GMRES restarting parameter. + * + * \e gmres-restarting-step-max - maximal adjusting step for the adaptivity of the GMRES restarting parameter. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ); + + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) override; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + using VectorType = typename Traits::VectorType; + // local vectors/views + using ConstDeviceView = typename Traits::ConstLocalViewType; + using DeviceView = typename Traits::LocalViewType; + using DeviceVector = typename Traits::LocalVectorType; + using HostView = typename DeviceView::template Self< RealType, Devices::Host >; + using HostVector = typename DeviceVector::template Self< RealType, Devices::Host >;; + + enum class Variant { CGS, CGSR, MGS, MGSR, CWY }; + + // nvcc allows __cuda_callable__ lambdas only in public methods + #ifdef __NVCC__ + public: + #endif + int orthogonalize_CGS( const int m, const RealType normb, const RealType beta ); + #ifdef __NVCC__ + protected: + #endif + + int orthogonalize_MGS( const int m, const RealType normb, const RealType beta ); + + int orthogonalize_CWY( const int m, const RealType normb, const RealType beta ); + + void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b ); + + void preconditioned_matvec( VectorViewType w, ConstVectorViewType v ); + + // nvcc allows __cuda_callable__ lambdas only in public methods + #ifdef __NVCC__ + public: + #endif + void hauseholder_generate( const int i, VectorViewType y_i, ConstVectorViewType z ); + #ifdef __NVCC__ + protected: + #endif - void hauseholder_cwy( VectorViewType v, - const int i ); - -// nvcc allows __cuda_callable__ lambdas only in public methods -#ifdef __NVCC__ -public: -#endif - void hauseholder_cwy_transposed( VectorViewType z, + void hauseholder_apply_trunc( HostView out, const int i, - ConstVectorViewType w ); -#ifdef __NVCC__ -protected: -#endif - - template< typename Vector > - void update( const int k, - const int m, - const HostVector& H, - const HostVector& s, - DeviceVector& V, - Vector& x ); - - void generatePlaneRotation( RealType& dx, - RealType& dy, - RealType& cs, - RealType& sn ); - - void applyPlaneRotation( RealType& dx, - RealType& dy, - RealType& cs, - RealType& sn ); - - void apply_givens_rotations( const int i, const int m ); - - void setSize( const VectorViewType& x ); - - // Specialized methods to distinguish between normal and distributed matrices - // in the implementation. - template< typename M > - static IndexType getLocalOffset( const M& m ) - { - return 0; - } - - template< typename M > - static IndexType getLocalOffset( const Matrices::DistributedMatrix< M >& m ) - { - return m.getLocalRowRange().getBegin(); - } - - // selected GMRES variant - Variant variant = Variant::CWY; - - // single vectors (distributed) - VectorType r, w, z, _M_tmp; - // matrices (in column-major format) (local) - DeviceVector V, Y; - // (CWY only) duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access - HostVector YL, T; - // host-only storage for Givens rotations and the least squares problem - HostVector cs, sn, H, s; - - IndexType size = 0; - IndexType ldSize = 0; - IndexType localOffset = 0; - int restarting_min = 10; - int restarting_max = 10; - int restarting_step_min = 3; - int restarting_step_max = 3; + VectorViewType y_i, + ConstVectorViewType z ); + + void hauseholder_cwy( VectorViewType v, + const int i ); + + // nvcc allows __cuda_callable__ lambdas only in public methods + #ifdef __NVCC__ + public: + #endif + void hauseholder_cwy_transposed( VectorViewType z, + const int i, + ConstVectorViewType w ); + #ifdef __NVCC__ + protected: + #endif + + template< typename Vector > + void update( const int k, + const int m, + const HostVector& H, + const HostVector& s, + DeviceVector& V, + Vector& x ); + + void generatePlaneRotation( RealType& dx, + RealType& dy, + RealType& cs, + RealType& sn ); + + void applyPlaneRotation( RealType& dx, + RealType& dy, + RealType& cs, + RealType& sn ); + + void apply_givens_rotations( const int i, const int m ); + + void setSize( const VectorViewType& x ); + + // Specialized methods to distinguish between normal and distributed matrices + // in the implementation. + template< typename M > + static IndexType getLocalOffset( const M& m ) + { + return 0; + } + + template< typename M > + static IndexType getLocalOffset( const Matrices::DistributedMatrix< M >& m ) + { + return m.getLocalRowRange().getBegin(); + } + + // selected GMRES variant + Variant variant = Variant::CWY; + + // single vectors (distributed) + VectorType r, w, z, _M_tmp; + // matrices (in column-major format) (local) + DeviceVector V, Y; + // (CWY only) duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access + HostVector YL, T; + // host-only storage for Givens rotations and the least squares problem + HostVector cs, sn, H, s; + + IndexType size = 0; + IndexType ldSize = 0; + IndexType localOffset = 0; + int restarting_min = 10; + int restarting_max = 10; + int restarting_step_min = 3; + int restarting_step_max = 3; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "GMRES.hpp" +#include diff --git a/src/TNL/Solvers/Linear/GMRES.hpp b/src/TNL/Solvers/Linear/GMRES.hpp index 4e6dcb81158f0a93649a2c311ffa2ab1176468ba..9e52d96eeb5053a706260b2a98257b2a02dde3ae 100644 --- a/src/TNL/Solvers/Linear/GMRES.hpp +++ b/src/TNL/Solvers/Linear/GMRES.hpp @@ -31,7 +31,7 @@ configSetup( Config::ConfigDescription& config, const String& prefix ) { LinearSolver< Matrix >::configSetup( config, prefix ); - config.addEntry< String >( prefix + "gmres-variant", "Minimal number of iterations after which the GMRES restarts.", "CWY" ); + config.addEntry< String >( prefix + "gmres-variant", "Algorithm used for the orthogonalization.", "MGSR" ); config.addEntryEnum( "CGS" ); config.addEntryEnum( "CGSR" ); config.addEntryEnum( "MGS" ); @@ -162,7 +162,6 @@ solve( ConstVectorViewType b, VectorViewType x ) if( o_steps < m ) { // exact solution has been reached early update( o_steps, m, H, s, V, x ); - this->refreshSolverMonitor( true ); return this->checkConvergence(); } @@ -179,8 +178,6 @@ solve( ConstVectorViewType b, VectorViewType x ) ++restart_cycles; beta_ratio = beta / beta_old; } - - this->refreshSolverMonitor( true ); return this->checkConvergence(); } @@ -273,8 +270,6 @@ orthogonalize_CGS( const int m, const RealType normb, const RealType beta ) this->setResidue( std::fabs( s[ i + 1 ] ) / normb ); if( ! this->checkNextIteration() ) return i; - else - this->refreshSolverMonitor(); } return m; @@ -347,8 +342,6 @@ orthogonalize_MGS( const int m, const RealType normb, const RealType beta ) this->setResidue( std::fabs( s[ i + 1 ] ) / normb ); if( ! this->checkNextIteration() ) return i; - else - this->refreshSolverMonitor(); } return m; @@ -434,9 +427,7 @@ orthogonalize_CWY( const int m, const RealType normb, const RealType beta ) this->setResidue( std::fabs( s[ i ] ) / normb ); if( i > 0 && ! this->checkNextIteration() ) return i - 1; - else - this->refreshSolverMonitor(); - } + } return m; } diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index f8261efe9ca4ce488977473a93c4b868c2582156..1cc061da98922dda7e855b2915e3ea448c00a7d1 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -8,85 +8,136 @@ #pragma once -#include "LinearSolver.h" - #include +#include #include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { +/** + * \brief Iterative solver of linear systems based on the Jacobi method. + * + * See [Wikipedia](https://en.wikipedia.org/wiki/Jacobi_method) for more details. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class Jacobi : public LinearSolver< Matrix > { using Base = LinearSolver< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ) - { - LinearSolver< Matrix >::configSetup( config, prefix ); - config.addEntry< double >( prefix + "jacobi-omega", "Relaxation parameter of the weighted/damped Jacobi method.", 1.0 ); - } - - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) override - { - if( parameters.checkParameter( prefix + "jacobi-omega" ) ) - this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) ); - if( this->omega <= 0.0 || this->omega > 2.0 ) - { - std::cerr << "Warning: The Jacobi method parameter omega is out of interval (0,2). The value is " << this->omega << " the method will not converge." << std::endl; - } - return LinearSolver< Matrix >::setup( parameters, prefix ); - } - - void setOmega( RealType omega ) - { - this->omega = omega; - } - - RealType getOmega() const - { - return omega; - } - - bool solve( ConstVectorViewType b, VectorViewType x ) override - { - const IndexType size = this->matrix->getRows(); - - Containers::Vector< RealType, DeviceType, IndexType > aux; - aux.setSize( size ); - - this->resetIterations(); - this->setResidue( this->getConvergenceResidue() + 1.0 ); - - RealType bNorm = b.lpNorm( ( RealType ) 2.0 ); - aux = x; - while( this->nextIteration() ) - { - for( IndexType row = 0; row < size; row ++ ) - this->matrix->performJacobiIteration( b, row, x, aux, this->getOmega() ); - for( IndexType row = 0; row < size; row ++ ) - this->matrix->performJacobiIteration( b, row, aux, x, this->getOmega() ); - this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - } - this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor(); - return this->checkConvergence(); - } - -protected: - RealType omega = 0.0; + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + */ + using DeviceType = typename Base::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief This is method defines configuration entries for setup of the linear iterative solver. + * + * In addition to config entries defined by \ref IterativeSolver::configSetup, this method + * defines the following: + * + * \e jacobi-omega - relaxation parameter of the weighted/damped Jacobi method - 1.0 by default. + * \e residue-period - number of iterations between subsequent recomputations of the residue - 4 by default. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, const String& prefix = "" ); + + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, const String& prefix = "" ) override; + + /** + * \brief Setter of the relaxation parameter. + * + * \param omega the relaxation parameter. It is 1 by default. + */ + void setOmega( RealType omega ); + + /** + * \brief Getter of the relaxation parameter. + * + * \return value of the relaxation parameter. + */ + RealType getOmega() const; + + /** + * \brief Set the period for a recomputation of the residue. + * + * \param period number of iterations between subsequent recomputations of the residue. + */ + void setResiduePeriod( IndexType period ); + + /** + * \brief Get the period for a recomputation of the residue. + * + * \return number of iterations between subsequent recomputations of the residue. + */ + IndexType getResiduePerid() const; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + + using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; + + RealType omega = 1.0; + + IndexType residuePeriod = 4; + + VectorType diagonal; + + public: // because nvcc does not accept lambda functions within private or protected methods + void performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + const ConstVectorViewType& in, + VectorViewType& out ) const; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL + +#include \ No newline at end of file diff --git a/src/TNL/Solvers/Linear/Jacobi.hpp b/src/TNL/Solvers/Linear/Jacobi.hpp new file mode 100644 index 0000000000000000000000000000000000000000..ba7ef968ca8e3f7521ba9b672b8189a0885a4992 --- /dev/null +++ b/src/TNL/Solvers/Linear/Jacobi.hpp @@ -0,0 +1,139 @@ +/*************************************************************************** + Jacobi.hpp - description + ------------------- + begin : Dec 22, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +#pragma once + +#include +#include +#include +#include + +namespace TNL { + namespace Solvers { + namespace Linear { + +template< typename Matrix > +void +Jacobi< Matrix >:: +configSetup( Config::ConfigDescription& config, const String& prefix ) +{ + LinearSolver< Matrix >::configSetup( config, prefix ); + config.addEntry< double >( prefix + "jacobi-omega", "Relaxation parameter of the weighted/damped Jacobi method.", 1.0 ); + config.addEntry< int >( prefix + "residue-period", "Number of iterations between subsequent recomputations of the residue.", 4 ); +} + +template< typename Matrix > +bool +Jacobi< Matrix >:: +setup( const Config::ParameterContainer& parameters, const String& prefix ) +{ + if( parameters.checkParameter( prefix + "jacobi-omega" ) ) + this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) ); + if( this->omega <= 0.0 || this->omega > 2.0 ) + { + std::cerr << "Warning: The Jacobi method parameter omega is out of interval (0,2). The value is " << this->omega << " the method will not converge." << std::endl; + } + if( parameters.checkParameter( prefix + "residue-period" ) ) + this->setResiduePeriod( parameters.getParameter< int >( prefix + "residue-period" ) ); + return LinearSolver< Matrix >::setup( parameters, prefix ); +} + +template< typename Matrix > +void +Jacobi< Matrix >:: +setOmega( RealType omega ) +{ + this->omega = omega; +} + +template< typename Matrix > +auto +Jacobi< Matrix >:: +getOmega() const -> RealType +{ + return omega; +} + +template< typename Matrix > +void +Jacobi< Matrix >:: +setResiduePeriod( IndexType period ) +{ + this->residuePeriod = period; +} + +template< typename Matrix > +auto +Jacobi< Matrix >:: +getResiduePerid() const -> IndexType +{ + return this->residuePeriod; +} + +template< typename Matrix > +bool +Jacobi< Matrix >:: +solve( ConstVectorViewType b, VectorViewType x ) +{ + const IndexType size = this->matrix->getRows(); + Containers::Vector< RealType, DeviceType, IndexType > aux; + aux.setSize( size ); + + ///// + // Fetch diagonal elements + this->diagonal.setSize( size ); + auto diagonalView = this->diagonal.getView(); + auto fetch_diagonal = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, const IndexType& columnIdx, const RealType& value ) mutable { + if( columnIdx == rowIdx ) diagonalView[ rowIdx ] = value; + }; + this->matrix->forAllElements( fetch_diagonal ); + + this->resetIterations(); + this->setResidue( this->getConvergenceResidue() + 1.0 ); + + auto bView = b.getView(); + auto xView = x.getView(); + auto auxView = aux.getView(); + RealType bNorm = lpNorm( b, ( RealType ) 2.0 ); + aux = x; + while( this->nextIteration() ) + { + this->performIteration( bView, diagonalView, xView, auxView ); + if( this->getIterations() % this->residuePeriod == 0 ) + this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); + this->currentIteration++; + this->performIteration( bView, diagonalView, auxView, xView ); + if( ( this->getIterations() ) % this->residuePeriod == 0 ) + this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); + + } + this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); + return this->checkConvergence(); +} + +template< typename Matrix > +void +Jacobi< Matrix >:: +performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + const ConstVectorViewType& in, + VectorViewType& out ) const +{ + const RealType omega_ = this->omega; + auto fetch = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, const RealType& value ) { + return value * in[ columnIdx ]; + }; + auto keep = [=] __cuda_callable__ ( IndexType rowIdx, const RealType& value ) mutable { + out[ rowIdx ] = in[ rowIdx ] + omega_ / diagonalView[ rowIdx ] * ( b[ rowIdx ] - value ); + }; + this->matrix->reduceAllRows( fetch, TNL::Plus{}, keep, 0.0 ); +} + + } // namespace Linear + } // namespace Solvers +} // namespace TNL diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index e80bad4b1a98a6cb86142ed1ae61aa6131aead37..f2de9a8f4a8c34cdc4755f365ed2f84f3e39ace3 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -20,58 +20,159 @@ #include namespace TNL { -namespace Solvers { + namespace Solvers { + namespace Linear { + /** - * \brief Namespace for linear system solvers. + * \brief Base class for iterative solvers of systems of linear equations. + * + * To use the linear solver, one needs to first set the matrix of the linear system + * by means of the method \ref LinearSolver::setMatrix. Afterward, one may call + * the method \ref LinearSolver::solve which accepts the right-hand side vector \e b + * and a vector \e x to which the solution will be stored. One may also use appropriate + * preconditioner to speed-up the convergence - see the method + * \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. */ -namespace Linear { - template< typename Matrix > class LinearSolver : public IterativeSolver< typename Matrix::RealType, typename Matrix::IndexType > { -public: - using RealType = typename Matrix::RealType; - using DeviceType = typename Matrix::DeviceType; - using IndexType = typename Matrix::IndexType; - using VectorViewType = typename Traits< Matrix >::VectorViewType; - using ConstVectorViewType = typename Traits< Matrix >::ConstVectorViewType; - using MatrixType = Matrix; - using MatrixPointer = std::shared_ptr< std::add_const_t< MatrixType > >; - using PreconditionerType = Preconditioners::Preconditioner< MatrixType >; - using PreconditionerPointer = std::shared_ptr< std::add_const_t< PreconditionerType > >; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ) - { - IterativeSolver< RealType, IndexType >::configSetup( config, prefix ); - } - - virtual bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) - { - return IterativeSolver< RealType, IndexType >::setup( parameters, prefix ); - } - - void setMatrix( const MatrixPointer& matrix ) - { - this->matrix = matrix; - } - - void setPreconditioner( const PreconditionerPointer& preconditioner ) - { - this->preconditioner = preconditioner; - } - - virtual bool solve( ConstVectorViewType b, VectorViewType x ) = 0; - - virtual ~LinearSolver() {} - -protected: - MatrixPointer matrix = nullptr; - PreconditionerPointer preconditioner = nullptr; -}; + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Matrix::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename Matrix::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Matrix::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Traits< Matrix >::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Traits< Matrix >::ConstVectorViewType; + + /** + * \brief Type of the matrix representing the linear system. + */ + using MatrixType = Matrix; + + /** + * \brief Type of shared pointer to the matrix. + */ + using MatrixPointer = std::shared_ptr< std::add_const_t< MatrixType > >; -} // namespace Linear -} // namespace Solvers + /** + * \brief Type of preconditioner. + */ + using PreconditionerType = Preconditioners::Preconditioner< MatrixType >; + + /** + * \brief Type of shared pointer to the preconditioner. + */ + using PreconditionerPointer = std::shared_ptr< std::add_const_t< PreconditionerType > >; + + /** + * \brief This method defines configuration entries for setup of the linear iterative solver. + * + * See \ref IterativeSolver::configSetup. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ) + { + IterativeSolver< RealType, IndexType >::configSetup( config, prefix ); + } + + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + virtual bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) + { + return IterativeSolver< RealType, IndexType >::setup( parameters, prefix ); + } + + /** + * \brief Set the matrix of the linear system. + * + * \param matrix is a shared pointer to the matrix of the linear system + */ + void setMatrix( const MatrixPointer& matrix ) + { + this->matrix = matrix; + } + + /** + * \brief Set the preconditioner. + * + * \param preconditioner is a shared pointer to preconditioner. + */ + void setPreconditioner( const PreconditionerPointer& preconditioner ) + { + this->preconditioner = preconditioner; + } + + /** + * \brief Method for solving a linear system. + * + * The linear system is defined by the matrix given by the method \ref LinearSolver::setMatrix and + * by the right-hand side vector represented by the vector \e b. The result is stored in the + * vector \e b. The solver can be accelerated with appropriate preconditioner set by the methods + * \ref LinearSolver::setPreconditioner. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + * + * \par Example + * \include Solvers/Linear/IterativeLinearSolverExample.cpp + * \par Output + * \include IterativeLinearSolverExample.out + */ + virtual bool solve( ConstVectorViewType b, VectorViewType x ) = 0; + + /** + * \brief Default destructor. + */ + virtual ~LinearSolver() {} + + protected: + MatrixPointer matrix = nullptr; + PreconditionerPointer preconditioner = nullptr; +}; + } // namespace Linear + } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h index 22bfceb007fa3107813ca069cf7efb861030afb7..a8775ca5f005c7381ce7c99189f1b7511bffe8ed 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h @@ -14,59 +14,162 @@ #include "Preconditioner.h" + namespace TNL { -namespace Solvers { -namespace Linear { -namespace Preconditioners { + namespace Solvers { + namespace Linear { + namespace Preconditioners { +/** + * \brief Diagonal (Jacobi) preconditioner for iterative solvers of linear systems. + * + * See [detailed description]([Netlib](http://netlib.org/linalg/html_templates/node55.html)). + * + * See \ref TNL::Solvers::Linear::Preconditioners::Preconditioner for example of setup with a linear solver. + * + * \tparam Matrix is type of the matrix describing the linear system. + */ template< typename Matrix > class Diagonal : public Preconditioner< Matrix > { -public: - using RealType = typename Matrix::RealType; - using DeviceType = typename Matrix::DeviceType; - using IndexType = typename Matrix::IndexType; - using typename Preconditioner< Matrix >::VectorViewType; - using typename Preconditioner< Matrix >::ConstVectorViewType; - using typename Preconditioner< Matrix >::MatrixPointer; - using VectorType = Containers::Vector< RealType, DeviceType, IndexType >; + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Matrix::RealType; + + /** + * \brief Device where the preconditioner will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename Matrix::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Matrix::IndexType; + + /** + * \brief Type for vector view. + */ + using typename Preconditioner< Matrix >::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using typename Preconditioner< Matrix >::ConstVectorViewType; - virtual void update( const MatrixPointer& matrixPointer ) override; + /** + * \brief Type of the matrix representing the linear system. + */ + using MatrixType = Matrix; - virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; + /** + * \brief Type of shared pointer to the matrix. + */ + using typename Preconditioner< Matrix >::MatrixPointer; -protected: - VectorType diagonal; + /** + * \brief This method updates the preconditioner with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + virtual void update( const MatrixPointer& matrixPointer ) override; + + /** + * \brief This method applies the preconditioner. + * + * \param b is the input vector the preconditioner is applied on. + * \param x is the result of the preconditioning. + */ + virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; + + protected: + + using VectorType = Containers::Vector< RealType, DeviceType, IndexType >; + + VectorType diagonal; }; +/** + * \brief Specialization of the diagonal preconditioner for distributed matrices. + * + * See \ref TNL::Solvers::Linear::Preconditioners::Diagonal + * + * \tparam Matrix is a type of matrix describing the linear system. + */ template< typename Matrix > class Diagonal< Matrices::DistributedMatrix< Matrix > > : public Preconditioner< Matrices::DistributedMatrix< Matrix > > { -public: - using MatrixType = Matrices::DistributedMatrix< Matrix >; - using RealType = typename MatrixType::RealType; - using DeviceType = typename MatrixType::DeviceType; - using IndexType = typename MatrixType::IndexType; - using typename Preconditioner< MatrixType >::VectorViewType; - using typename Preconditioner< MatrixType >::ConstVectorViewType; - using typename Preconditioner< MatrixType >::MatrixPointer; - using VectorType = Containers::Vector< RealType, DeviceType, IndexType >; - using LocalViewType = Containers::VectorView< RealType, DeviceType, IndexType >; - using ConstLocalViewType = Containers::VectorView< std::add_const_t< RealType >, DeviceType, IndexType >; - - virtual void update( const MatrixPointer& matrixPointer ) override; - - virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; - -protected: - VectorType diagonal; + public: + + /** + * \brief Type of the matrix representing the linear system. + */ + using MatrixType = Matrices::DistributedMatrix< Matrix >; + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename MatrixType::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename MatrixType::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename MatrixType::IndexType; + + /** + * \brief Type for vector view. + */ + using typename Preconditioner< MatrixType >::VectorViewType; + + /** + * \brief Type for vector constant view. + */ + using typename Preconditioner< MatrixType >::ConstVectorViewType; + + /** + * \brief Type of shared pointer to the matrix. + */ + using typename Preconditioner< MatrixType >::MatrixPointer; + + /** + * \brief This method updates the preconditioner with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + virtual void update( const MatrixPointer& matrixPointer ) override; + + /** + * \brief This method applies the preconditioner. + * + * \param b is the input vector the preconditioner is applied on. + * \param x is the result of the preconditioning. + */ + virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; + + protected: + using VectorType = Containers::Vector< RealType, DeviceType, IndexType >; + using LocalViewType = Containers::VectorView< RealType, DeviceType, IndexType >; + using ConstLocalViewType = Containers::VectorView< std::add_const_t< RealType >, DeviceType, IndexType >; + + VectorType diagonal; }; -} // namespace Preconditioners -} // namespace Linear -} // namespace Solvers + } // namespace Preconditioners + } // namespace Linear + } // namespace Solvers } // namespace TNL #include "Diagonal.hpp" diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h index b623e86fbd74faada4afed7c0d09c0c872907610..a6d128c394a2040d1ba9ac664736892d1c01b0c5 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h @@ -24,9 +24,9 @@ #endif namespace TNL { -namespace Solvers { -namespace Linear { -namespace Preconditioners { + namespace Solvers { + namespace Linear { + namespace Preconditioners { // implementation template template< typename Matrix, typename Real, typename Device, typename Index > @@ -52,51 +52,101 @@ public: } }; -// actual template to be used by users +/** + * \brief Implementation of a preconditioner based on Incomplete LU. + * + * See [detailed description](https://en.wikipedia.org/wiki/Incomplete_LU_factorization). + * + * \tparam Matrix is type of the matrix describing the linear system. + */ template< typename Matrix > class ILU0 : public ILU0_impl< Matrix, typename Matrix::RealType, typename Matrix::DeviceType, typename Matrix::IndexType > {}; +/** + * \brief Implementation of a preconditioner based in Incomplete LU - specialization for CPU. + * + * See [detailed description](https://en.wikipedia.org/wiki/Incomplete_LU_factorization). + * + * See \ref TNL::Solvers::Linear::Preconditioners::Preconditioner for example of setup with a linear solver. + * + * \tparam Matrix is type of the matrix describing the linear system. + */ template< typename Matrix, typename Real, typename Index > class ILU0_impl< Matrix, Real, Devices::Host, Index > : public Preconditioner< Matrix > { -public: - using RealType = Real; - using DeviceType = Devices::Host; - using IndexType = Index; - using typename Preconditioner< Matrix >::VectorViewType; - using typename Preconditioner< Matrix >::ConstVectorViewType; - using typename Preconditioner< Matrix >::MatrixPointer; - - virtual void update( const MatrixPointer& matrixPointer ) override; - - virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; - -protected: - // The factors L and U are stored separately and the rows of U are reversed. - using CSR = Matrices::SparseMatrix< RealType, DeviceType, IndexType, Matrices::GeneralMatrix, Algorithms::Segments::CSRScalar >; - CSR L, U; - - // Specialized methods to distinguish between normal and distributed matrices - // in the implementation. - template< typename M > - static IndexType getMinColumn( const M& m ) - { - return 0; - } - - template< typename M > - static IndexType getMinColumn( const Matrices::DistributedMatrix< M >& m ) - { - if( m.getRows() == m.getColumns() ) - // square matrix, assume global column indices - return m.getLocalRowRange().getBegin(); - else - // non-square matrix, assume ghost indexing + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = Real; + + /** + * \brief Device where the preconditioner will run on and auxillary data will alloacted on. + */ + using DeviceType = Devices::Host; + + /** + * \brief Type for indexing. + */ + using IndexType = Index; + + /** + * \brief Type for vector view. + */ + using typename Preconditioner< Matrix >::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using typename Preconditioner< Matrix >::ConstVectorViewType; + + /** + * \brief Type of shared pointer to the matrix. + */ + using typename Preconditioner< Matrix >::MatrixPointer; + + /** + * \brief This method updates the preconditioner with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + virtual void update( const MatrixPointer& matrixPointer ) override; + + /** + * \brief This method applies the preconditioner. + * + * \param b is the input vector the preconditioner is applied on. + * \param x is the result of the preconditioning. + */ + virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; + + protected: + // The factors L and U are stored separately and the rows of U are reversed. + using CSR = Matrices::SparseMatrix< RealType, DeviceType, IndexType, Matrices::GeneralMatrix, Algorithms::Segments::CSRScalar >; + CSR L, U; + + // Specialized methods to distinguish between normal and distributed matrices + // in the implementation. + template< typename M > + static IndexType getMinColumn( const M& m ) + { return 0; - } + } + + template< typename M > + static IndexType getMinColumn( const Matrices::DistributedMatrix< M >& m ) + { + if( m.getRows() == m.getColumns() ) + // square matrix, assume global column indices + return m.getLocalRowRange().getBegin(); + else + // non-square matrix, assume ghost indexing + return 0; + } }; template< typename Matrix > diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h index 88f5441bac5b7289532f275ba3a8939ebd2523de..5e09147f98165a93b49712966279bd71c22fc275 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h @@ -27,18 +27,79 @@ template< typename Matrix, typename Real, typename Device, typename Index > class ILUT_impl {}; -// actual template to be used by users +/** + * \brief Implementation of a preconditioner based on Incomplete LU with thresholding. + * + * See [detailed description](https://www-users.cse.umn.edu/~saad/PDF/umsi-92-38.pdf) + * + * See \ref TNL::Solvers::Linear::Preconditioners::Preconditioner for example of setup with a linear solver. + * + * \tparam Matrix is type of the matrix describing the linear system. + */ template< typename Matrix > class ILUT : public ILUT_impl< Matrix, typename Matrix::RealType, typename Matrix::DeviceType, typename Matrix::IndexType > { -public: - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ) - { - config.addEntry< int >( prefix + "ilut-p", "Number of additional non-zero entries to allocate on each row of the factors L and U.", 0 ); - config.addEntry< double >( prefix + "ilut-threshold", "Threshold for droppping small entries.", 1e-4 ); - } + using Base = ILUT_impl< Matrix, typename Matrix::RealType, typename Matrix::DeviceType, typename Matrix::IndexType >; + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Matrix::RealType; + + /** + * \brief Device where the preconditioner will run on and auxillary data will alloacted on. + */ + using DeviceType = Devices::Host; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Matrix::IndexType; + + /** + * \brief Type for vector view. + */ + using typename Preconditioner< Matrix >::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using typename Preconditioner< Matrix >::ConstVectorViewType; + + /** + * \brief Type of shared pointer to the matrix. + */ + using typename Preconditioner< Matrix >::MatrixPointer; + + /** + * \brief This method defines configuration entries for setup of the preconditioner of linear iterative solver. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ) + { + config.addEntry< int >( prefix + "ilut-p", "Number of additional non-zero entries to allocate on each row of the factors L and U.", 0 ); + config.addEntry< double >( prefix + "ilut-threshold", "Threshold for droppping small entries.", 1e-4 ); + } + + /** + * \brief This method updates the preconditioner with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + using Base::update; + + /** + * \brief This method applies the preconditioner. + * + * \param b is the input vector the preconditioner is applied on. + * \param x is the result of the preconditioning. + */ + using Base::solve; }; template< typename Matrix, typename Real, typename Index > @@ -88,6 +149,12 @@ protected: } }; +template< typename Matrix, typename Real, typename Index > +class ILUT_impl< Matrix, Real, Devices::Sequential, Index > +: public ILUT_impl< Matrix, Real, Devices::Host, Index > +{}; + + template< typename Matrix, typename Real, typename Index > class ILUT_impl< Matrix, Real, Devices::Cuda, Index > : public Preconditioner< Matrix > diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h index 672f5167e0b2a7ff87e1551e837884ff7dc8dbdb..4357dcf63ec0e4a2eff9233420cde8d82ac32497 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h @@ -19,47 +19,113 @@ #include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { + namespace Preconditioners { + /** - * \brief Namespace for preconditioners for linear system solvers. + * \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 */ -namespace Preconditioners { - template< typename Matrix > class Preconditioner { -public: - using RealType = typename Matrix::RealType; - using DeviceType = typename Matrix::DeviceType; - using IndexType = typename Matrix::IndexType; - using VectorViewType = typename Traits< Matrix >::VectorViewType; - using ConstVectorViewType = typename Traits< Matrix >::ConstVectorViewType; - using MatrixType = Matrix; - using MatrixPointer = std::shared_ptr< std::add_const_t< MatrixType > >; - - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ) - {} - - virtual bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) - { - return true; - } - - virtual void update( const MatrixPointer& matrixPointer ) - {} - - virtual void solve( ConstVectorViewType b, VectorViewType x ) const - { - throw std::logic_error("The solve() method of a dummy preconditioner should not be called."); - } - - virtual ~Preconditioner() {} + public: + + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Matrix::RealType; + + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + * + * See \ref Devices::Host or \ref Devices::Cuda. + */ + using DeviceType = typename Matrix::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Matrix::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Traits< Matrix >::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Traits< Matrix >::ConstVectorViewType; + + /** + * \brief Type of the matrix representing the linear system. + */ + using MatrixType = Matrix; + + /** + * \brief Type of shared pointer to the matrix. + */ + using MatrixPointer = std::shared_ptr< std::add_const_t< MatrixType > >; + + /** + * \brief This method defines configuration entries for setup of the preconditioner of linear iterative solver. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ) {} + + /** + * \brief Method for setup of the preconditioner of linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + virtual bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) + { + return true; + } + + /** + * \brief This method updates the preconditioner with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + virtual void update( const MatrixPointer& matrixPointer ) + {} + + /** + * \brief This method applies the preconditioner. + * + * \param b is the input vector the preconditioner is applied on. + * \param x is the result of the preconditioning. + */ + virtual void solve( ConstVectorViewType b, VectorViewType x ) const + { + throw std::logic_error("The solve() method of a dummy preconditioner should not be called."); + } + + /** + * \brief Destructor of the preconditioner. + */ + virtual ~Preconditioner() {} }; -} // namespace Preconditioners -} // namespace Linear -} // namespace Solvers + } // namespace Preconditioners + } // namespace Linear + } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/Linear/Preconditioners/_NamespaceDoxy.h b/src/TNL/Solvers/Linear/Preconditioners/_NamespaceDoxy.h new file mode 100644 index 0000000000000000000000000000000000000000..bd29562ff3850ed50d7ac0eb859219ae26b4db83 --- /dev/null +++ b/src/TNL/Solvers/Linear/Preconditioners/_NamespaceDoxy.h @@ -0,0 +1,29 @@ +/*************************************************************************** + _NamespaceDoxy.h - description + ------------------- + begin : Dec 21, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +namespace TNL { + namespace Solvers { + namespace Linear { + /** + * \brief Namespace for preconditioners of linear system solvers. + * + * This namespace contains the following preconditioners for iterative solvers linear systems. + * + * 1. Diagonal - is diagonal or Jacobi preconditioner - see[Netlib](http://netlib.org/linalg/html_templates/node55.html) + * 2. ILU0 - is Incomplete LU preconditioner with the same sparsity pattern as the original matrix - see [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_LU_factorization) + * 3. ILUT - is Incomplete LU preconiditoner with thresholding - see [paper by Y. Saad](https://www-users.cse.umn.edu/~saad/PDF/umsi-92-38.pdf) + */ + namespace Preconditioners { + } // namespace Preconditioners + } // namespace Linear + } // namespace Solvers +} // namespace TNL diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index 8bc1a171694a41f91b7010e138cc8219b91c5e8f..76ff1db97d596e661969ee8d853da739ce9ac2db 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -10,42 +10,135 @@ #pragma once -#include "LinearSolver.h" +#include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { +/** + * \brief Iterative solver of linear systems based on the Successive-overrelaxation (SOR) or Gauss-Seidel method. + * + * See (Wikipedia)[https://en.wikipedia.org/wiki/Successive_over-relaxation] for more details. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class SOR : public LinearSolver< Matrix > { using Base = LinearSolver< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; - static void configSetup( Config::ConfigDescription& config, - const String& prefix = "" ); + public: - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) override; + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; - void setOmega( const RealType& omega ); + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + */ + using DeviceType = typename Base::DeviceType; - const RealType& getOmega() const; + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; - bool solve( ConstVectorViewType b, VectorViewType x ) override; + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief This is method defines configuration entries for setup of the linear iterative solver. + * + * In addition to config entries defined by \ref IterativeSolver::configSetup, this method + * defines the following: + * + * \e sor-omega - relaxation parameter of the weighted/damped Jacobi method. + * + * \param config contains description of configuration parameters. + * \param prefix is a prefix of particular configuration entries. + */ + static void configSetup( Config::ConfigDescription& config, + const String& prefix = "" ); + + /** + * \brief Method for setup of the linear iterative solver based on configuration parameters. + * + * \param parameters contains values of the define configuration entries. + * \param prefix is a prefix of particular configuration entries. + */ + bool setup( const Config::ParameterContainer& parameters, + const String& prefix = "" ) override; + + /** + * \brief Setter of the relaxation parameter. + * + * \param omega the relaxation parameter. It is 1 by default. + */ + void setOmega( const RealType& omega ); + + /** + * \brief Getter of the relaxation parameter. + * + * \return value of the relaxation parameter. + */ + const RealType& getOmega() const; + + /** + * \brief Set the period for a recomputation of the residue. + * + * \param period number of iterations between subsequent recomputations of the residue. + */ + void setResiduePeriod( IndexType period ); + + /** + * \brief Get the period for a recomputation of the residue. + * + * \return number of iterations between subsequent recomputations of the residue. + */ + IndexType getResiduePerid() const; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; + + RealType omega = 1.0; + + IndexType residuePeriod = 4; + + VectorType diagonal; + + public: // because nvcc does not accept lambda functions within private or protected methods + void performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + VectorViewType& x ) const; -protected: - RealType omega = 1.0; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "SOR.hpp" +#include diff --git a/src/TNL/Solvers/Linear/SOR.hpp b/src/TNL/Solvers/Linear/SOR.hpp index 30101c08b831ed78a994d739642660f77f50da94..c7580830b5717502a7bea3bae5f5215056f1cb86 100644 --- a/src/TNL/Solvers/Linear/SOR.hpp +++ b/src/TNL/Solvers/Linear/SOR.hpp @@ -10,7 +10,8 @@ #pragma once -#include "SOR.h" +#include +#include #include namespace TNL { @@ -25,6 +26,7 @@ configSetup( Config::ConfigDescription& config, { LinearSolver< Matrix >::configSetup( config, prefix ); config.addEntry< double >( prefix + "sor-omega", "Relaxation parameter of the SOR method.", 1.0 ); + config.addEntry< int >( prefix + "residue-period", "Says after how many iterations the reside is recomputed.", 4 ); } template< typename Matrix > @@ -39,6 +41,8 @@ setup( const Config::ParameterContainer& parameters, { std::cerr << "Warning: The SOR method parameter omega is out of interval (0,2). The value is " << this->omega << " the method will not converge." << std::endl; } + if( parameters.checkParameter( prefix + "residue-period" ) ) + this->setResiduePeriod( parameters.getParameter< int >( prefix + "residue-period" ) ); return LinearSolver< Matrix >::setup( parameters, prefix ); } @@ -54,28 +58,68 @@ const typename SOR< Matrix > :: RealType& SOR< Matrix > :: getOmega( ) const return this->omega; } +template< typename Matrix > +void +SOR< Matrix >:: +setResiduePeriod( IndexType period ) +{ + this->residuePeriod = period; +} + +template< typename Matrix > +auto +SOR< Matrix >:: +getResiduePerid() const -> IndexType +{ + return this->residuePeriod; +} + template< typename Matrix > bool SOR< Matrix > :: solve( ConstVectorViewType b, VectorViewType x ) { + ///// + // Fetch diagonal elements const IndexType size = this->matrix->getRows(); + this->diagonal.setSize( size ); + auto diagonalView = this->diagonal.getView(); + auto fetch_diagonal = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, const IndexType& columnIdx, const RealType& value ) mutable { + if( columnIdx == rowIdx ) diagonalView[ rowIdx ] = value; + }; + this->matrix->forAllElements( fetch_diagonal ); this->resetIterations(); this->setResidue( this->getConvergenceResidue() + 1.0 ); - RealType bNorm = lpNorm( b, 2.0 ); - + auto bView = b.getView(); + auto xView = x.getView(); + RealType bNorm = lpNorm( b, ( RealType ) 2.0 ); while( this->nextIteration() ) { - for( IndexType row = 0; row < size; row ++ ) - this->matrix->performSORIteration( b, row, x, this->getOmega() ); - this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor(); + this->performIteration( bView, diagonalView, xView ); + if( this->getIterations() % this->residuePeriod == 0 ) + this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); } this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor( true ); return this->checkConvergence(); }; +template< typename Matrix > +void +SOR< Matrix >:: +performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + VectorViewType& x ) const +{ + const RealType omega_ = this->omega; + auto fetch = [=] __cuda_callable__ ( IndexType rowIdx, IndexType columnIdx, const RealType& value ) { + return value * x[ columnIdx ]; + }; + auto keep = [=] __cuda_callable__ ( IndexType rowIdx, const RealType& value ) mutable { + Algorithms::AtomicOperations< DeviceType >::add( x[ rowIdx ], omega_ / diagonalView[ rowIdx ] * ( b[ rowIdx ] - value ) ); + }; + this->matrix->reduceAllRows( fetch, TNL::Plus{}, keep, 0.0 ); +} + } // namespace Linear } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/Linear/TFQMR.h b/src/TNL/Solvers/Linear/TFQMR.h index d267819da7f80d7be624b21c34e9c464e2bd0919..97a2438829318cad0fac0c91fb4af3424fb35de9 100644 --- a/src/TNL/Solvers/Linear/TFQMR.h +++ b/src/TNL/Solvers/Linear/TFQMR.h @@ -10,34 +10,74 @@ #pragma once -#include "LinearSolver.h" +#include namespace TNL { -namespace Solvers { -namespace Linear { + namespace Solvers { + namespace Linear { +/** + * \brief Iterative solver of linear systems based on the Transpose-free quasi-minimal residual (TFQMR) method. + * + * See (Wikipedia)[https://second.wiki/wiki/algoritmo_tfqmr] for more details. + * + * See \ref TNL::Solvers::Linear::IterativeSolver for example of showing how to use the linear solvers. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class TFQMR : public LinearSolver< Matrix > { using Base = LinearSolver< Matrix >; -public: - using RealType = typename Base::RealType; - using DeviceType = typename Base::DeviceType; - using IndexType = typename Base::IndexType; - using VectorViewType = typename Base::VectorViewType; - using ConstVectorViewType = typename Base::ConstVectorViewType; - bool solve( ConstVectorViewType b, VectorViewType x ) override; + public: -protected: - void setSize( const VectorViewType& x ); + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Base::RealType; - typename Traits< Matrix >::VectorType d, r, w, u, v, r_ast, Au, M_tmp; + /** + * \brief Device where the solver will run on and auxillary data will alloacted on. + */ + using DeviceType = typename Base::DeviceType; + + /** + * \brief Type for indexing. + */ + using IndexType = typename Base::IndexType; + + /** + * \brief Type for vector view. + */ + using VectorViewType = typename Base::VectorViewType; + + /** + * \brief Type for constant vector view. + */ + using ConstVectorViewType = typename Base::ConstVectorViewType; + + /** + * \brief Method for solving of a linear system. + * + * See \ref LinearSolver::solve for more details. + * + * \param b vector with the right-hand side of the linear system. + * \param x vector for the solution of the linear system. + * \return true if the solver converged. + * \return false if the solver did not converge. + */ + bool solve( ConstVectorViewType b, VectorViewType x ) override; + + protected: + void setSize( const VectorViewType& x ); + + typename Traits< Matrix >::VectorType d, r, w, u, v, r_ast, Au, M_tmp; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "TFQMR.hpp" +#include diff --git a/src/TNL/Solvers/Linear/TFQMR.hpp b/src/TNL/Solvers/Linear/TFQMR.hpp index 9c2a1f205e6342f5297b833c033418c90b143127..d1277956577dd4fa9dedbe540f4096e9c47f2f0d 100644 --- a/src/TNL/Solvers/Linear/TFQMR.hpp +++ b/src/TNL/Solvers/Linear/TFQMR.hpp @@ -112,11 +112,7 @@ bool TFQMR< Matrix >::solve( ConstVectorViewType b, VectorViewType x ) else { u -= alpha * v; } - - this->refreshSolverMonitor(); } - - this->refreshSolverMonitor( true ); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/_NamespaceDoxy.h b/src/TNL/Solvers/Linear/_NamespaceDoxy.h new file mode 100644 index 0000000000000000000000000000000000000000..5471fc405a5440285e848424f5ece26666aa1b9d --- /dev/null +++ b/src/TNL/Solvers/Linear/_NamespaceDoxy.h @@ -0,0 +1,45 @@ +/*************************************************************************** + _NamespaceDoxy.h - description + ------------------- + begin : Dec 21, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +namespace TNL { + namespace Solvers { + + /** + * \brief Namespace for linear system solvers. + * + * This namespace contains the following algorithms and methods for solution of linear systems. + * + * # Direct methods + * + * # Iterative methods + * + * ## Stationary methods + * 1. Jacobi method - \ref TNL::Solvers::Linear::Jacobi + * 2. Successive-overrelaxation method, SOR - \ref TNL::Solvers::Linear::SOR + * + * ## Krylov subspace methods + * 1. Conjugate gradient method, CG - \ref TNL::Solvers::Linear::CG + * 2. Biconjugate gradient stabilized method, BICGStab - \ref TNL::Solvers::Linear::BICGStab + * 3. BICGStab(l) method - \ref TNL::Solvers::Linear::BICGStabL + * 4. Transpose-free quasi-minimal residual method, TFQMR - \ref TNL::Solvers::Linear::TFQMR + * 5. Generalized minimal residual method, GMERS - \ref TNL::Solvers::Linear::GMRES with various methods of orthogonalization + * 1. [Classical Gramm-Schmidt, CGS](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process) + * 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 + * + */ + namespace Linear { + } // namespace Linear + } // namespace Solvers +} // namespace TNL diff --git a/src/TNL/Solvers/LinearSolverTypeResolver.h b/src/TNL/Solvers/LinearSolverTypeResolver.h index 5615ea7077bad28229b21b8a45c9eb297c5aa3d6..db6fb5121a1b8a4511ce49c57e6539cd070d40ce 100644 --- a/src/TNL/Solvers/LinearSolverTypeResolver.h +++ b/src/TNL/Solvers/LinearSolverTypeResolver.h @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -28,12 +29,18 @@ #include namespace TNL { -namespace Solvers { + namespace Solvers { +/** + * \brief Function returning available linear solvers. + * + * \return container with names of available linear solvers. + */ inline std::vector getLinearSolverOptions() { return { + "jacobi", "sor", "cg", "bicgstab", @@ -46,6 +53,11 @@ getLinearSolverOptions() }; } +/** + * \brief Function returning available linear preconditioners. + * + * \return container with names of available linear preconditioners. + */ inline std::vector getPreconditionerOptions() { @@ -57,10 +69,34 @@ getPreconditionerOptions() }; } +/** + * \brief Function returning shared pointer with linear solver given by its name in a form of a string. + * + * \tparam MatrixType is a type of matrix defining the system of linear equations. + * \param name of the linear solver. The name can be one of the following: + * 1. `jacobi` - for the Jacobi solver - \ref TNL::Solvers::Linear::Jacobi. + * 2. `sor` - for SOR solver - \ref TNL::Solvers::Linear::SOR. + * 3. `cg` - for CG solver - \ref TNL::Solvers::Linear::CG. + * 4. `bicgstab` - for BICGStab solver - \ref TNL::Solvers::Linear::BICGStab. + * 5. `bicgstabl` - for BICGStabL solver - \ref TNL::Solvers::Linear::BICGStabL. + * 6. `gmres` - for GMRES solver - \ref TNL::Solvers::Linear::GMRES. + * 7. `tfqmr` - for TFQMR solver - \ref TNL::Solvers::Linear::TFQMR. + * \return shared pointer with given linear solver. + * + * The following example shows how to use this function: + * + * \includelineno Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverWithRuntimeTypesExample.out + */ template< typename MatrixType > std::shared_ptr< Linear::LinearSolver< MatrixType > > getLinearSolver( std::string name ) { + if( name == "jacobi" ) + return std::make_shared< Linear::Jacobi< MatrixType > >(); if( name == "sor" ) return std::make_shared< Linear::SOR< MatrixType > >(); if( name == "cg" ) @@ -90,6 +126,23 @@ getLinearSolver( std::string name ) return nullptr; } +/** + * \brief Function returning shared pointer with linear preconditioner given by its name in a form of a string. + * + * \tparam MatrixType is a type of matrix defining the system of linear equations. + * \param name of the linear preconditioner. The name can be one of the following: + * 1. `none` - for none preconditioner + * 2. `diagonal` - for diagonal/Jacobi preconditioner - \ref TNL::Solvers::Linear::Preconditioners::Diagonal + * 3. `ilu0` - for ILU(0) preconditioner - \ref TNL::Solvers::Linear::Preconditioners::ILU0 + * 4. `ilut` - for ILU with thresholding preconditioner - \ref TNL::Solvers::Linear::Preconditioners::ILUT + * \return shared pointer with given linear preconditioner. + * + * \includelineno Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverWithRuntimeTypesExample.out + */ template< typename MatrixType > std::shared_ptr< Linear::Preconditioners::Preconditioner< MatrixType > > getPreconditioner( std::string name ) @@ -116,5 +169,5 @@ getPreconditioner( std::string name ) return nullptr; } -} // namespace Solvers + } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/ODE/Euler.hpp b/src/TNL/Solvers/ODE/Euler.hpp index 4532feb4d588a83c1f701add169c3340380c885f..eff39fda59b97ea39a039cf8dedd3ac0a30f81a8 100644 --- a/src/TNL/Solvers/ODE/Euler.hpp +++ b/src/TNL/Solvers/ODE/Euler.hpp @@ -83,8 +83,6 @@ bool Euler< Problem, SolverMonitor > :: solve( DofVectorPointer& _u ) this->resetIterations(); this->setResidue( this->getConvergenceResidue() + 1.0 ); - this -> refreshSolverMonitor(); - /**** * Start the main loop */ @@ -126,17 +124,12 @@ bool Euler< Problem, SolverMonitor > :: solve( DofVectorPointer& _u ) currentTau = this -> getStopTime() - time; //we don't want to keep such tau else this -> tau = currentTau; - this->refreshSolverMonitor(); - /**** * Check stop conditions. */ if( time >= this->getStopTime() || ( this -> getConvergenceResidue() != 0.0 && this->getResidue() < this -> getConvergenceResidue() ) ) - { - this -> refreshSolverMonitor(); return true; - } if( this -> cflCondition != 0.0 ) { diff --git a/src/TNL/Solvers/ODE/Merson.hpp b/src/TNL/Solvers/ODE/Merson.hpp index c48e0484237b5afd2e198da75236d3ae115cecf5..88edd0aded1ce03d40ad6639b191e9f24cfc77d3 100644 --- a/src/TNL/Solvers/ODE/Merson.hpp +++ b/src/TNL/Solvers/ODE/Merson.hpp @@ -108,8 +108,6 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& _u ) this->resetIterations(); this->setResidue( this->getConvergenceResidue() + 1.0 ); - this->refreshSolverMonitor(); - ///// // Start the main loop while( this->checkNextIteration() ) @@ -175,7 +173,6 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& _u ) if( ! this->nextIteration() ) return false; } - this->refreshSolverMonitor(); ///// // Compute the new time step. @@ -195,14 +192,9 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& _u ) // Check stop conditions. if( time >= this->getStopTime() || ( this->getConvergenceResidue() != 0.0 && this->getResidue() < this->getConvergenceResidue() ) ) - { - this->refreshSolverMonitor( true ); return true; - } } - this->refreshSolverMonitor( true ); return this->checkConvergence(); - }; template< typename Problem, typename SolverMonitor > diff --git a/src/TNL/Solvers/SolverMonitor.h b/src/TNL/Solvers/SolverMonitor.h index ff349224b958fa3b299ea0d473ce0d1c46600c56..fca151cec428113e9859ca404e6d8d68a60656b7 100644 --- a/src/TNL/Solvers/SolverMonitor.h +++ b/src/TNL/Solvers/SolverMonitor.h @@ -16,108 +16,143 @@ #include namespace TNL { -namespace Solvers { - + namespace Solvers { + +/** + * \brief Base class for solver monitors. + * + * The solver monitors serve for monitoring a convergence and status of various solvers. + * The solver monitor uses separate thread for monitoring the solver status in preset time period. + */ class SolverMonitor { -public: - SolverMonitor() - : timeout_milliseconds( 500 ), - started( false ), - stopped( false ), - timer( nullptr ) - {} - - virtual void refresh() = 0; - - void setRefreshRate( const int& refreshRate ) - { - timeout_milliseconds = refreshRate; - } - - void setTimer( Timer& timer ) - { - this->timer = &timer; - } - - void runMainLoop() - { - // We need to use both 'started' and 'stopped' to avoid a deadlock - // when the loop thread runs this method delayed after the - // SolverMonitorThread's destructor has already called stopMainLoop() - // from the main thread. - started = true; - - const int timeout_base = 100; - const std::chrono::milliseconds timeout( timeout_base ); - - while( ! stopped ) { - refresh(); - - // make sure to detect changes to refresh rate - int steps = timeout_milliseconds / timeout_base; - if( steps <= 0 ) - steps = 1; - - int i = 0; - while( ! stopped && i++ < steps ) { - std::this_thread::sleep_for( timeout ); - } - } - - // reset to initial state - started = false; - stopped = false; - } + public: - void stopMainLoop() - { - stopped = true; - } + /** + * \brief Basic construct with no arguments + */ + SolverMonitor() + : timeout_milliseconds( 500 ), + started( false ), + stopped( false ), + timer( nullptr ) + {} + + /** + * \brief This abstract method is responsible for printing or visualizing the status of the solver. + */ + virtual void refresh() = 0; + + /** + * \brief Set the time interval between two consecutive calls of \ref SolverMonitor::refresh. + * + * \param refreshRate refresh rate in miliseconds. + */ + void setRefreshRate( const int& refreshRate ) { timeout_milliseconds = refreshRate; } + + /** + * \brief Set a timer object for the solver monitor. + * + * If a timer is set, the monitor can measure real elapsed time since the start of the solver. + * + * \param timer is an instance of \ref TNL::Timer. + */ + void setTimer( Timer& timer ) { this->timer = &timer; } + + /** + * \brief Starts the main loop from which the method \ref SolverMonitor::refresh is called in given time periods. + */ + void runMainLoop() + { + // We need to use both 'started' and 'stopped' to avoid a deadlock + // when the loop thread runs this method delayed after the + // SolverMonitorThread's destructor has already called stopMainLoop() + // from the main thread. + started = true; + + const int timeout_base = 100; + const std::chrono::milliseconds timeout( timeout_base ); + + while( ! stopped ) { + refresh(); + + // make sure to detect changes to refresh rate + int steps = timeout_milliseconds / timeout_base; + if( steps <= 0 ) + steps = 1; + + int i = 0; + while( ! stopped && i++ < steps ) { + std::this_thread::sleep_for( timeout ); + } + } - bool isStopped() const - { - return stopped; - } + // reset to initial state + started = false; + stopped = false; + } -protected: - double getElapsedTime() - { - if( ! timer ) - return 0.0; - return timer->getRealTime(); - } + /** + * \brief Stops the main loop of the monitor. See \ref TNL::SolverMonitor::runMainLoop. + */ + void stopMainLoop() { stopped = true; } + + /** + * \brief Checks whether the main loop was stopped. + * + * \return true if the main loop was stopped. + * \return false if the main loop was not stopped yet. + */ + bool isStopped() const { return stopped; } + + protected: + double getElapsedTime() + { + if( ! timer ) + return 0.0; + return timer->getRealTime(); + } - std::atomic_int timeout_milliseconds; + std::atomic_int timeout_milliseconds; - std::atomic_bool started, stopped; + std::atomic_bool started, stopped; - Timer* timer; + Timer* timer; }; -// a RAII wrapper for launching the SolverMonitor's main loop in a separate thread +/** + * \brief A RAII wrapper for launching the SolverMonitor's main loop in a separate thread. + */ class SolverMonitorThread { public: - SolverMonitorThread( SolverMonitor& solverMonitor ) - : solverMonitor( solverMonitor ), - t( &SolverMonitor::runMainLoop, &solverMonitor ) - {} - - ~SolverMonitorThread() - { - solverMonitor.stopMainLoop(); - if( t.joinable() ) - t.join(); - } + /** + * \brief Constructor with instance of solver monitor. + * + * \param solverMonitor is a reference to an instance of a solver monitor. + */ + SolverMonitorThread( SolverMonitor& solverMonitor ) + : solverMonitor( solverMonitor ), + t( &SolverMonitor::runMainLoop, &solverMonitor ) + {} + + /** + * \brief Destructor. + */ + ~SolverMonitorThread() + { + solverMonitor.stopMainLoop(); + if( t.joinable() ) + t.join(); + } private: - SolverMonitor& solverMonitor; + SolverMonitor& solverMonitor; - std::thread t; + std::thread t; }; -} // namespace Solvers + } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/_NamespaceDoxy.h b/src/TNL/Solvers/_NamespaceDoxy.h new file mode 100644 index 0000000000000000000000000000000000000000..e54373dd593a87dca8d8733e6813140d645e4f76 --- /dev/null +++ b/src/TNL/Solvers/_NamespaceDoxy.h @@ -0,0 +1,21 @@ +/*************************************************************************** + _NamespaceDoxy.h - description + ------------------- + begin : Dec 21, 2021 + copyright : (C) 2021 by Tomas Oberhuber et al. + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once + +namespace TNL { + + /** + * \brief Namespace for solvers. + */ + namespace Solvers { + + } // namespace Solvers +} // namespace TNL