From ab741a8bbda6627d20502e04892b5ddead475015 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 20 Dec 2021 19:49:03 +0100 Subject: [PATCH 01/54] Deleting whitespaces in File.h. --- src/TNL/File.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/TNL/File.h b/src/TNL/File.h index cef110e16..9048c022e 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 > > -- GitLab From 31b80016b3479be607b742edd289f879db2cc0d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 20 Dec 2021 19:50:39 +0100 Subject: [PATCH 02/54] Writting documentation for IterativeSolver. --- src/TNL/Solvers/IterativeSolver.h | 284 ++++++++++++++++++++++-------- 1 file changed, 215 insertions(+), 69 deletions(-) diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index 97d3f1c56..0015c0dad 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -19,82 +19,228 @@ namespace TNL { 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 \d 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 alowed 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( 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; }; -} // namespace Solvers + } // namespace Solvers } // namespace TNL #include -- GitLab From 2e8f1ee1d2074ff8292997ac1d770b99108ca0bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 20 Dec 2021 19:51:22 +0100 Subject: [PATCH 03/54] Deleting parameter force from IterativeSolver::refreshSolverMonitor. --- src/Benchmarks/ODESolvers/Merson.hpp | 4 ++-- src/TNL/Solvers/IterativeSolver.h | 2 +- src/TNL/Solvers/IterativeSolver.hpp | 4 ++-- src/TNL/Solvers/Linear/BICGStab.hpp | 2 +- src/TNL/Solvers/Linear/BICGStabL.hpp | 2 +- src/TNL/Solvers/Linear/CG.hpp | 2 +- src/TNL/Solvers/Linear/GMRES.hpp | 4 ++-- src/TNL/Solvers/Linear/SOR.hpp | 2 +- src/TNL/Solvers/Linear/TFQMR.hpp | 2 +- src/TNL/Solvers/ODE/Merson.hpp | 4 ++-- 10 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Benchmarks/ODESolvers/Merson.hpp b/src/Benchmarks/ODESolvers/Merson.hpp index 3bb6f9c96..8100feec7 100644 --- a/src/Benchmarks/ODESolvers/Merson.hpp +++ b/src/Benchmarks/ODESolvers/Merson.hpp @@ -222,11 +222,11 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& u ) if( time >= this->getStopTime() || ( this->getConvergenceResidue() != 0.0 && this->getResidue() < this->getConvergenceResidue() ) ) { - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return true; } } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); }; diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index 0015c0dad..08c8dfc36 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -215,7 +215,7 @@ class IterativeSolver /** * \brief Refreshes the solver monitor. */ - void refreshSolverMonitor( bool force = false ); + void refreshSolverMonitor(); protected: Index maxIterations = 1000000000; diff --git a/src/TNL/Solvers/IterativeSolver.hpp b/src/TNL/Solvers/IterativeSolver.hpp index f894f3711..1c6995043 100644 --- a/src/TNL/Solvers/IterativeSolver.hpp +++ b/src/TNL/Solvers/IterativeSolver.hpp @@ -249,13 +249,13 @@ setSolverMonitor( SolverMonitorType& 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 ); } } diff --git a/src/TNL/Solvers/Linear/BICGStab.hpp b/src/TNL/Solvers/Linear/BICGStab.hpp index afbdbc399..be03fa545 100644 --- a/src/TNL/Solvers/Linear/BICGStab.hpp +++ b/src/TNL/Solvers/Linear/BICGStab.hpp @@ -119,7 +119,7 @@ solve( ConstVectorViewType b, VectorViewType x ) } } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/BICGStabL.hpp b/src/TNL/Solvers/Linear/BICGStabL.hpp index 7ae4999c9..5f965d682 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.hpp +++ b/src/TNL/Solvers/Linear/BICGStabL.hpp @@ -239,7 +239,7 @@ solve( ConstVectorViewType b, VectorViewType x ) } } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/CG.hpp b/src/TNL/Solvers/Linear/CG.hpp index 67442b4b1..618e05106 100644 --- a/src/TNL/Solvers/Linear/CG.hpp +++ b/src/TNL/Solvers/Linear/CG.hpp @@ -107,7 +107,7 @@ solve( ConstVectorViewType b, VectorViewType x ) this->setResidue( std::sqrt(s1) / normb ); } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/GMRES.hpp b/src/TNL/Solvers/Linear/GMRES.hpp index 4e6dcb811..3ec79c2af 100644 --- a/src/TNL/Solvers/Linear/GMRES.hpp +++ b/src/TNL/Solvers/Linear/GMRES.hpp @@ -162,7 +162,7 @@ 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 ); + this->refreshSolverMonitor(); return this->checkConvergence(); } @@ -180,7 +180,7 @@ solve( ConstVectorViewType b, VectorViewType x ) beta_ratio = beta / beta_old; } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/SOR.hpp b/src/TNL/Solvers/Linear/SOR.hpp index 30101c08b..959c88516 100644 --- a/src/TNL/Solvers/Linear/SOR.hpp +++ b/src/TNL/Solvers/Linear/SOR.hpp @@ -72,7 +72,7 @@ bool SOR< Matrix > :: solve( ConstVectorViewType b, VectorViewType x ) this->refreshSolverMonitor(); } this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); }; diff --git a/src/TNL/Solvers/Linear/TFQMR.hpp b/src/TNL/Solvers/Linear/TFQMR.hpp index 9c2a1f205..c1644afe6 100644 --- a/src/TNL/Solvers/Linear/TFQMR.hpp +++ b/src/TNL/Solvers/Linear/TFQMR.hpp @@ -116,7 +116,7 @@ bool TFQMR< Matrix >::solve( ConstVectorViewType b, VectorViewType x ) this->refreshSolverMonitor(); } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/ODE/Merson.hpp b/src/TNL/Solvers/ODE/Merson.hpp index c48e04842..4a576c57b 100644 --- a/src/TNL/Solvers/ODE/Merson.hpp +++ b/src/TNL/Solvers/ODE/Merson.hpp @@ -196,11 +196,11 @@ bool Merson< Problem, SolverMonitor >::solve( DofVectorPointer& _u ) if( time >= this->getStopTime() || ( this->getConvergenceResidue() != 0.0 && this->getResidue() < this->getConvergenceResidue() ) ) { - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return true; } } - this->refreshSolverMonitor( true ); + this->refreshSolverMonitor(); return this->checkConvergence(); }; -- GitLab From aeff4024cdadac4882117797dde46827a08b7e6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 20 Dec 2021 20:10:02 +0100 Subject: [PATCH 04/54] Removing useless IterativeSolver::refreshSolverMonitor. --- src/Benchmarks/ODESolvers/Euler.hpp | 7 ------- src/Benchmarks/ODESolvers/Merson.hpp | 7 ------- src/TNL/Solvers/IterativeSolver.h | 2 +- src/TNL/Solvers/IterativeSolver.hpp | 7 +++---- src/TNL/Solvers/Linear/BICGStab.hpp | 1 - src/TNL/Solvers/Linear/BICGStabL.hpp | 2 -- src/TNL/Solvers/Linear/CG.hpp | 1 - src/TNL/Solvers/Linear/GMRES.hpp | 11 +---------- src/TNL/Solvers/Linear/Jacobi.h | 1 - src/TNL/Solvers/Linear/SOR.hpp | 2 -- src/TNL/Solvers/Linear/TFQMR.hpp | 4 ---- src/TNL/Solvers/ODE/Euler.hpp | 7 ------- src/TNL/Solvers/ODE/Merson.hpp | 8 -------- 13 files changed, 5 insertions(+), 55 deletions(-) diff --git a/src/Benchmarks/ODESolvers/Euler.hpp b/src/Benchmarks/ODESolvers/Euler.hpp index 840c0a45a..72fb4867a 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 8100feec7..2fa53ed92 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(); return true; - } } - this->refreshSolverMonitor(); return this->checkConvergence(); }; diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index 08c8dfc36..af3d7f571 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -215,7 +215,7 @@ class IterativeSolver /** * \brief Refreshes the solver monitor. */ - void refreshSolverMonitor(); + //void refreshSolverMonitor(); protected: Index maxIterations = 1000000000; diff --git a/src/TNL/Solvers/IterativeSolver.hpp b/src/TNL/Solvers/IterativeSolver.hpp index 1c6995043..0a78dff89 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,7 +245,7 @@ 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() @@ -257,7 +256,7 @@ refreshSolverMonitor() this->solverMonitor->setResidue( this->getResidue() ); this->solverMonitor->setRefreshRate( this->refreshRate ); } -} +}*/ } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/Linear/BICGStab.hpp b/src/TNL/Solvers/Linear/BICGStab.hpp index be03fa545..9cf51034f 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(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/BICGStabL.hpp b/src/TNL/Solvers/Linear/BICGStabL.hpp index 5f965d682..b68141590 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(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/CG.hpp b/src/TNL/Solvers/Linear/CG.hpp index 618e05106..660e3a68d 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(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/GMRES.hpp b/src/TNL/Solvers/Linear/GMRES.hpp index 3ec79c2af..cf8e8bafa 100644 --- a/src/TNL/Solvers/Linear/GMRES.hpp +++ b/src/TNL/Solvers/Linear/GMRES.hpp @@ -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(); return this->checkConvergence(); } @@ -179,8 +178,6 @@ solve( ConstVectorViewType b, VectorViewType x ) ++restart_cycles; beta_ratio = beta / beta_old; } - - this->refreshSolverMonitor(); 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 f8261efe9..5016d2330 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -79,7 +79,6 @@ public: this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); } this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/Linear/SOR.hpp b/src/TNL/Solvers/Linear/SOR.hpp index 959c88516..c77af2a33 100644 --- a/src/TNL/Solvers/Linear/SOR.hpp +++ b/src/TNL/Solvers/Linear/SOR.hpp @@ -69,10 +69,8 @@ bool SOR< Matrix > :: solve( ConstVectorViewType b, VectorViewType x ) 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->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) ); - this->refreshSolverMonitor(); return this->checkConvergence(); }; diff --git a/src/TNL/Solvers/Linear/TFQMR.hpp b/src/TNL/Solvers/Linear/TFQMR.hpp index c1644afe6..d12779565 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(); return this->checkConvergence(); } diff --git a/src/TNL/Solvers/ODE/Euler.hpp b/src/TNL/Solvers/ODE/Euler.hpp index 4532feb4d..eff39fda5 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 4a576c57b..88edd0ade 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(); return true; - } } - this->refreshSolverMonitor(); return this->checkConvergence(); - }; template< typename Problem, typename SolverMonitor > -- GitLab From 2e4a52c3b4ad48bb66f2fca7b4c4a89eb060153d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 20 Dec 2021 22:26:36 +0100 Subject: [PATCH 05/54] Fixing typos. --- src/TNL/Solvers/IterativeSolver.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index af3d7f571..e943b3fc6 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -17,7 +17,11 @@ #include namespace TNL { -namespace Solvers { + + /** + * \brief Namespace for solvers. + */ + namespace Solvers { /** * \brief Base class for iterative solvers. @@ -135,7 +139,7 @@ class IterativeSolver /** * \brief Sets the residue limit for the divergence criterion. * - * The divergence occurs when the residue \d exceeds the limit. + * The divergence occurs when the residue \b exceeds the limit. * * \param divergenceResidue the residue limit of the divergence. */ @@ -196,7 +200,7 @@ class IterativeSolver bool nextIteration(); /** - * \brief Checks if the solver is alowed to the next iteration. + * \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 -- GitLab From d3b452143d10c56ddb8208a6b4f9a4b0f471790d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:22:22 +0100 Subject: [PATCH 06/54] Writtin documentation for iterative solver. --- src/TNL/Solvers/IterativeSolver.h | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index e943b3fc6..a0e1ccd0d 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -17,10 +17,6 @@ #include namespace TNL { - - /** - * \brief Namespace for solvers. - */ namespace Solvers { /** @@ -231,7 +227,7 @@ class IterativeSolver Real convergenceResidue = 1e-6; // If the current residue is greater than divergenceResidue, the solver is stopped. - Real divergenceResidue = std::numeric_limits< float >::max(); + Real divergenceResidue = std::numeric_limits< Real >::max(); Real currentResidue = 0; -- GitLab From 9288c2bffb8d0be1bf7cb7f023d4cd194972e406 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:22:42 +0100 Subject: [PATCH 07/54] Writtin documentation for iterative solver monitor. --- src/TNL/Solvers/IterativeSolverMonitor.h | 142 ++++++++++++++++------- 1 file changed, 102 insertions(+), 40 deletions(-) diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h index 07dee8de9..02a224703 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.h +++ b/src/TNL/Solvers/IterativeSolverMonitor.h @@ -13,51 +13,113 @@ #include namespace TNL { -namespace Solvers { - + 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. + */ template< typename Real, typename Index> 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; + + IndexType verbose; + + IndexType nodesPerIteration; }; -} // namespace Solvers + } // namespace Solvers } // namespace TNL #include -- GitLab From faa4166ee1faac909ee4aedeffd11d22435537cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:23:06 +0100 Subject: [PATCH 08/54] Editing documentation for linear iterative solver. --- src/TNL/Solvers/Linear/LinearSolver.h | 179 +++++++++++++++++++------- 1 file changed, 133 insertions(+), 46 deletions(-) diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index e80bad4b1..c684d3abd 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -20,58 +20,145 @@ #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. */ -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 representign 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 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 is 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. + */ + 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 Linear + } // namespace Solvers } // namespace TNL -- GitLab From 310e04258c78c3f6a6f191fb863b1ec11e52d078 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:23:33 +0100 Subject: [PATCH 09/54] Writtin documentation for CG method. --- src/TNL/Solvers/Linear/CG.h | 73 +++++++++++++++++++++++++++++-------- 1 file changed, 57 insertions(+), 16 deletions(-) diff --git a/src/TNL/Solvers/Linear/CG.h b/src/TNL/Solvers/Linear/CG.h index f2cad3f05..ecabc0296 100644 --- a/src/TNL/Solvers/Linear/CG.h +++ b/src/TNL/Solvers/Linear/CG.h @@ -10,35 +10,76 @@ #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. + */ 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 -- GitLab From 48503fb846ba2cf875323afb0890551d55471b3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:23:50 +0100 Subject: [PATCH 10/54] Writtin documentation for BICGStab method. --- src/TNL/Solvers/Linear/BICGStab.h | 107 +++++++++++++++++++++++------- 1 file changed, 84 insertions(+), 23 deletions(-) diff --git a/src/TNL/Solvers/Linear/BICGStab.h b/src/TNL/Solvers/Linear/BICGStab.h index 233590ed7..2d1e7d3a1 100644 --- a/src/TNL/Solvers/Linear/BICGStab.h +++ b/src/TNL/Solvers/Linear/BICGStab.h @@ -10,46 +10,107 @@ #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. + * + * \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 -- GitLab From ddc7649817ff3c528ece619ebd89954dbdb28983 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:24:01 +0100 Subject: [PATCH 11/54] Writtin documentation for BICGStab(l) method. --- src/TNL/Solvers/Linear/BICGStabL.h | 168 +++++++++++++++++++++-------- 1 file changed, 122 insertions(+), 46 deletions(-) diff --git a/src/TNL/Solvers/Linear/BICGStabL.h b/src/TNL/Solvers/Linear/BICGStabL.h index f30a23634..682d14370 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.h +++ b/src/TNL/Solvers/Linear/BICGStabL.h @@ -41,12 +41,39 @@ #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 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 + * proposed by van der Vorst [3]. BiCGstab(1) is equivalent to BiCGstab, and + * BiCGstab(2) is a slightly more efficient version of the BiCGstab2 algorithm + * by Gutknecht [4], while BiCGstab(l>2) is a further generalization. + * + * [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). + + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class BICGStabL : public LinearSolver< Matrix > @@ -56,50 +83,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 -- GitLab From 5b58b810038c358ea6c3245e4313b7ff357e41e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:24:51 +0100 Subject: [PATCH 12/54] Fixing documentation for GMRES variant config entry. --- src/TNL/Solvers/Linear/GMRES.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TNL/Solvers/Linear/GMRES.hpp b/src/TNL/Solvers/Linear/GMRES.hpp index cf8e8bafa..5fe58fc02 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 to be used for reorthogonalization.", "CWY" ); config.addEntryEnum( "CGS" ); config.addEntryEnum( "CGSR" ); config.addEntryEnum( "MGS" ); -- GitLab From 1ef9febfd42d7a9d67eacaea72a6f5e332c92a3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:25:06 +0100 Subject: [PATCH 13/54] Writtin documentation for GMRES method. --- src/TNL/Solvers/Linear/GMRES.h | 341 ++++++++++++++++++++------------- 1 file changed, 209 insertions(+), 132 deletions(-) diff --git a/src/TNL/Solvers/Linear/GMRES.h b/src/TNL/Solvers/Linear/GMRES.h index 7314a3a05..7ff51cbb9 100644 --- a/src/TNL/Solvers/Linear/GMRES.h +++ b/src/TNL/Solvers/Linear/GMRES.h @@ -12,12 +12,33 @@ #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. + * + * \tparam Matrix is type of matrix describing the linear system. + */ template< typename Matrix > class GMRES : public LinearSolver< Matrix > @@ -25,139 +46,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 -- GitLab From 43ff25126865843a637cfd03d5b86f5059a3202d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:25:20 +0100 Subject: [PATCH 14/54] Writtin documentation for TFQMR method. --- src/TNL/Solvers/Linear/TFQMR.h | 70 ++++++++++++++++++++++++++-------- 1 file changed, 54 insertions(+), 16 deletions(-) diff --git a/src/TNL/Solvers/Linear/TFQMR.h b/src/TNL/Solvers/Linear/TFQMR.h index d267819da..e9ec68de5 100644 --- a/src/TNL/Solvers/Linear/TFQMR.h +++ b/src/TNL/Solvers/Linear/TFQMR.h @@ -10,34 +10,72 @@ #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. + * + * \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 -- GitLab From 9deda63310793b921f9357d817b52f4fd1da52d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:26:11 +0100 Subject: [PATCH 15/54] Splitting Jacobi.h file and writing documentation for the Jacobi method. --- src/TNL/Solvers/Linear/Jacobi.h | 158 +++++++++++++++++------------- src/TNL/Solvers/Linear/Jacobi.hpp | 87 ++++++++++++++++ 2 files changed, 177 insertions(+), 68 deletions(-) create mode 100644 src/TNL/Solvers/Linear/Jacobi.hpp diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index 5016d2330..0cdf57845 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -8,84 +8,106 @@ #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. + * + * \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 ) ); - 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. + * + * \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 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: + RealType omega = 0.0; }; -} // 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 000000000..eee368973 --- /dev/null +++ b/src/TNL/Solvers/Linear/Jacobi.hpp @@ -0,0 +1,87 @@ +/*************************************************************************** + Jacobi.hpp - description + ------------------- + begin : Dec 22, 2021 + copyright : (C) 2021 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +#pragma once + +#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 ); +} + +template< typename Matrix > +bool +Jacobi< Matrix >:: +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 ); +} + +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 > +bool +Jacobi< Matrix >:: +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 ) ); + return this->checkConvergence(); +} + + } // namespace Linear + } // namespace Solvers +} // namespace TNL -- GitLab From 8ef460b6134e708aa75c6be637d6efbcec54cd04 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:26:28 +0100 Subject: [PATCH 16/54] Writtin documentation for SOR method. --- src/TNL/Solvers/Linear/SOR.h | 107 ++++++++++++++++++++++++++++------- 1 file changed, 86 insertions(+), 21 deletions(-) diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index 8bc1a1716..f11f5b788 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -10,42 +10,107 @@ #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) method. + * + * See (Wikipedia)[https://en.wikipedia.org/wiki/Successive_over-relaxation] for more details. + * + * \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; -protected: - RealType omega = 1.0; + /** + * \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 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: + RealType omega = 1.0; }; -} // namespace Linear -} // namespace Solvers + } // namespace Linear + } // namespace Solvers } // namespace TNL -#include "SOR.hpp" +#include -- GitLab From 566e0bba0d903ee8793f5d54f8dfef8d04a80343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:27:04 +0100 Subject: [PATCH 17/54] Writting documentation for TNL::Solvers namespace. --- src/TNL/Solvers/_NamespaceDoxy.h | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 src/TNL/Solvers/_NamespaceDoxy.h diff --git a/src/TNL/Solvers/_NamespaceDoxy.h b/src/TNL/Solvers/_NamespaceDoxy.h new file mode 100644 index 000000000..e54373dd5 --- /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 -- GitLab From f32bd54e7ec15c203a4a9067202f508b518c7490 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 22 Dec 2021 16:28:11 +0100 Subject: [PATCH 18/54] Writting documentation for TNL::Solvers::Linear namespace. --- src/TNL/Solvers/Linear/_NamespaceDoxy.h | 38 +++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 src/TNL/Solvers/Linear/_NamespaceDoxy.h diff --git a/src/TNL/Solvers/Linear/_NamespaceDoxy.h b/src/TNL/Solvers/Linear/_NamespaceDoxy.h new file mode 100644 index 000000000..ef4147b72 --- /dev/null +++ b/src/TNL/Solvers/Linear/_NamespaceDoxy.h @@ -0,0 +1,38 @@ +/*************************************************************************** + _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 (SOR) method - \ref TNL::Solvers::Linear::SOR + * + * ## Krylov subspace methods + * 1. Conjugate gradient (CG) method - \ref TNL::Solvers::Linear::CG + * 2. Biconjugate gradient stabilised (BICGStab) method - \ref TNL::Solvers::Linear::BICGStab + * 3. Generalized minimal residual (GMRES) method - \ref TNL::Solvers::Linear::GMRES + * 4. Transpose-free quasi-minimal residual (TFQMR) method - \ref TNL::Solvers::Linear::TFQMR + */ + namespace Linear { + } // namespace Linear + } // namespace Solvers +} // namespace TNL -- GitLab From 9d61ceacda3746471ef94e523b73f02705d08190 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 23 Dec 2021 12:55:24 +0100 Subject: [PATCH 19/54] Writting tutorial on iterative linear solvers. --- Documentation/Examples/CMakeLists.txt | 1 + Documentation/Examples/Solvers/CMakeLists.txt | 1 + .../Examples/Solvers/Linear/CMakeLists.txt | 24 ++++++ .../Linear/IterativeLinearSolverExample.cpp | 82 +++++++++++++++++++ .../Linear/IterativeLinearSolverExample.cu | 1 + Documentation/Tutorials/CMakeLists.txt | 1 + .../Tutorials/Solvers/CMakeLists.txt | 1 + .../Tutorials/Solvers/Linear/CMakeLists.txt | 33 ++++++++ .../Solvers/Linear/tutorial_Linear_solvers.md | 38 +++++++++ Documentation/Tutorials/index.md | 5 +- src/TNL/Solvers/Linear/GMRES.hpp | 2 +- src/TNL/Solvers/Linear/LinearSolver.h | 5 ++ src/TNL/Solvers/Linear/SOR.h | 2 +- src/TNL/Solvers/Linear/_NamespaceDoxy.h | 17 ++-- 14 files changed, 204 insertions(+), 9 deletions(-) create mode 100644 Documentation/Examples/Solvers/CMakeLists.txt create mode 100644 Documentation/Examples/Solvers/Linear/CMakeLists.txt create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cu create mode 100644 Documentation/Tutorials/Solvers/CMakeLists.txt create mode 100644 Documentation/Tutorials/Solvers/Linear/CMakeLists.txt create mode 100644 Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 7aa736429..b27168ca0 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 000000000..79b9712aa --- /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 000000000..575542019 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -0,0 +1,24 @@ +set( COMMON_EXAMPLES + IterativeLinearSolverExample +) + +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 000000000..186d67cce --- /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 ); + std::cout << "Vector b = " << b << std::endl; + + /*** + * Solver the linear system + */ + using LinearSolver = TNL::Solvers::Linear::GMRES< MatrixType >; + LinearSolver solver; + solver.setMatrix( matrix_ptr ); + x = 0.0; + 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::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/IterativeLinearSolverExample.cu b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cu new file mode 120000 index 000000000..17884b901 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverExample.cpp \ No newline at end of file diff --git a/Documentation/Tutorials/CMakeLists.txt b/Documentation/Tutorials/CMakeLists.txt index 05ed1f33c..e30aab1e2 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 000000000..79b9712aa --- /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 000000000..71facac2f --- /dev/null +++ b/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt @@ -0,0 +1,33 @@ +IF( BUILD_CUDA ) + CUDA_ADD_EXECUTABLE( ArrayAllocation ArrayAllocation.cu ) + ADD_CUSTOM_COMMAND( COMMAND ArrayAllocation > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayAllocation.out OUTPUT ArrayAllocation.out ) + CUDA_ADD_EXECUTABLE( ArrayIO ArrayIO.cu ) + ADD_CUSTOM_COMMAND( COMMAND ArrayIO > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayIO.out OUTPUT ArrayIO.out ) + CUDA_ADD_EXECUTABLE( ArrayView-1 ArrayView-1.cu ) + ADD_CUSTOM_COMMAND( COMMAND ArrayView-1 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayView-1.out OUTPUT ArrayView-1.out ) + CUDA_ADD_EXECUTABLE( ArrayView-2 ArrayView-2.cu ) + ADD_CUSTOM_COMMAND( COMMAND ArrayView-2 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayView-2.out OUTPUT ArrayView-2.out ) + CUDA_ADD_EXECUTABLE( ArrayViewForElements ArrayViewForElements.cu ) + ADD_CUSTOM_COMMAND( COMMAND ArrayViewForElements > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayViewForElements.out OUTPUT ArrayViewForElements.out ) + CUDA_ADD_EXECUTABLE( contains contains.cu ) + ADD_CUSTOM_COMMAND( COMMAND contains > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/contains.out OUTPUT contains.out ) + CUDA_ADD_EXECUTABLE( ElementsAccessing-1 ElementsAccessing-1.cu ) + ADD_CUSTOM_COMMAND( COMMAND ElementsAccessing-1 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ElementsAccessing-1.out OUTPUT ElementsAccessing-1.out ) + CUDA_ADD_EXECUTABLE( ElementsAccessing-2 ElementsAccessing-2.cu ) + ADD_CUSTOM_COMMAND( COMMAND ElementsAccessing-2 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ElementsAccessing-2.out OUTPUT ElementsAccessing-2.out ) + ADD_EXECUTABLE( StaticArrayExample StaticArrayExample.cpp ) + ADD_CUSTOM_COMMAND( COMMAND StaticArrayExample > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticArrayExample.out OUTPUT StaticArrayExample.out ) +ENDIF() + +IF( BUILD_CUDA ) +ADD_CUSTOM_TARGET( TutorialsArrays-cuda ALL DEPENDS + ArrayAllocation.out + ArrayIO.out + ArrayView-1.out + ArrayView-2.out + contains.out + ElementsAccessing-1.out + ElementsAccessing-2.out + ArrayViewForElements.out + StaticArrayExample.out ) +ENDIF() 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 000000000..2df2f7319 --- /dev/null +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -0,0 +1,38 @@ +\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) - CPU only currently +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 Krylov subspace methods can be combined with the following precoditioners: + +1. Jacobi +2. ILU - CPU only currently + +# Iterative solvers of linear systems + +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 + +The result looks as follows: + +\include IterativeLinearSolverExample.out \ No newline at end of file diff --git a/Documentation/Tutorials/index.md b/Documentation/Tutorials/index.md index 031de3fae..900641f1b 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/TNL/Solvers/Linear/GMRES.hpp b/src/TNL/Solvers/Linear/GMRES.hpp index 5fe58fc02..9e52d96ee 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", "Algorithm to be used for reorthogonalization.", "CWY" ); + config.addEntry< String >( prefix + "gmres-variant", "Algorithm used for the orthogonalization.", "MGSR" ); config.addEntryEnum( "CGS" ); config.addEntryEnum( "CGSR" ); config.addEntryEnum( "MGS" ); diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index c684d3abd..4f6dc22b9 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -146,6 +146,11 @@ class LinearSolver * \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; diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index f11f5b788..c5fc4c8aa 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -17,7 +17,7 @@ namespace TNL { namespace Linear { /** - * \brief Iterative solver of linear systems based on the Successive-overrelaxation (SOR) method. + * \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. * diff --git a/src/TNL/Solvers/Linear/_NamespaceDoxy.h b/src/TNL/Solvers/Linear/_NamespaceDoxy.h index ef4147b72..5471fc405 100644 --- a/src/TNL/Solvers/Linear/_NamespaceDoxy.h +++ b/src/TNL/Solvers/Linear/_NamespaceDoxy.h @@ -24,13 +24,20 @@ namespace TNL { * * ## Stationary methods * 1. Jacobi method - \ref TNL::Solvers::Linear::Jacobi - * 2. Successive-overrelaxation (SOR) method - \ref TNL::Solvers::Linear::SOR + * 2. Successive-overrelaxation method, SOR - \ref TNL::Solvers::Linear::SOR * * ## Krylov subspace methods - * 1. Conjugate gradient (CG) method - \ref TNL::Solvers::Linear::CG - * 2. Biconjugate gradient stabilised (BICGStab) method - \ref TNL::Solvers::Linear::BICGStab - * 3. Generalized minimal residual (GMRES) method - \ref TNL::Solvers::Linear::GMRES - * 4. Transpose-free quasi-minimal residual (TFQMR) method - \ref TNL::Solvers::Linear::TFQMR + * 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 -- GitLab From c5b307bb78450d7ecad8e0a0b9a2741acdc6f932 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Sun, 26 Dec 2021 19:47:14 +0100 Subject: [PATCH 20/54] Reimplementing the Jacobi solver. --- src/TNL/Solvers/Linear/Jacobi.h | 17 ++++++++-- src/TNL/Solvers/Linear/Jacobi.hpp | 52 +++++++++++++++++++++++++------ 2 files changed, 58 insertions(+), 11 deletions(-) diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index 0cdf57845..cebc84472 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -61,7 +61,8 @@ class Jacobi * 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. + * \e jacobi-omega - relaxation parameter of the weighted/damped Jacobi method - 1.0 by default. + * \e residue-period - says after how many iterations the reside is recomputed - 4 by default. * * \param config contains description of configuration parameters. * \param prefix is a prefix of particular configuration entries. @@ -103,7 +104,19 @@ class Jacobi bool solve( ConstVectorViewType b, VectorViewType x ) override; protected: - RealType omega = 0.0; + + using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; + + void performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + const ConstVectorViewType& in, + VectorViewType& out ) const; + + RealType omega = 1.0; + + IndexType residuePeriod = 4; + + VectorType diagonal; }; } // namespace Linear diff --git a/src/TNL/Solvers/Linear/Jacobi.hpp b/src/TNL/Solvers/Linear/Jacobi.hpp index eee368973..9a900fb6f 100644 --- a/src/TNL/Solvers/Linear/Jacobi.hpp +++ b/src/TNL/Solvers/Linear/Jacobi.hpp @@ -11,6 +11,7 @@ #include #include #include +#include namespace TNL { namespace Solvers { @@ -23,12 +24,13 @@ 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", "Says after how many iterations the reside is recomputed.", 4 ); } template< typename Matrix > bool Jacobi< Matrix >:: -setup( const Config::ParameterContainer& parameters, const String& prefix ) override +setup( const Config::ParameterContainer& parameters, const String& prefix ) { if( parameters.checkParameter( prefix + "jacobi-omega" ) ) this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) ); @@ -58,30 +60,62 @@ getOmega() const -> RealType template< typename Matrix > bool Jacobi< Matrix >:: -solve( ConstVectorViewType b, VectorViewType x ) override +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 ); - RealType bNorm = b.lpNorm( ( RealType ) 2.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() ) { - 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->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 -- GitLab From 45d5810f4aabc1b865583f3fd4196efcc94e2c6a Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Sun, 26 Dec 2021 22:01:57 +0100 Subject: [PATCH 21/54] Added explicit conversion to RealType in SparseMatrixView::forElements. --- src/TNL/Matrices/SparseMatrixView.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TNL/Matrices/SparseMatrixView.hpp b/src/TNL/Matrices/SparseMatrixView.hpp index 2af26efd5..d7ed73815 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 ] ); } -- GitLab From 0038a2ce2a3887e3e720fc0ddd27fb1b0cdb4533 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Sun, 26 Dec 2021 22:02:50 +0100 Subject: [PATCH 22/54] Fixing type in iterative linear solver example. --- .../Examples/Solvers/Linear/IterativeLinearSolverExample.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp index 186d67cce..229aab369 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp @@ -57,15 +57,15 @@ void iterativeLinearSolverExample() Vector x( size, 1.0 ); Vector b( size ); matrix_ptr->vectorProduct( x, b ); + x = 0.0; std::cout << "Vector b = " << b << std::endl; /*** - * Solver the linear system + * Solve the linear system */ using LinearSolver = TNL::Solvers::Linear::GMRES< MatrixType >; LinearSolver solver; solver.setMatrix( matrix_ptr ); - x = 0.0; solver.solve( b, x ); std::cout << "Vector x = " << x << std::endl; } -- GitLab From e29a0c9a43f5f1caba5fe653de06f30522e53684 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Sun, 26 Dec 2021 22:03:43 +0100 Subject: [PATCH 23/54] Added default values of template parameters of IterativeSolverMonitor. --- src/TNL/Solvers/IterativeSolverMonitor.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h index 02a224703..401e34c23 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.h +++ b/src/TNL/Solvers/IterativeSolverMonitor.h @@ -21,7 +21,8 @@ namespace TNL { * \tparam Real is a type of the floating-point arithmetics. * \tparam Index is an indexing type. */ -template< typename Real, typename Index> +template< typename Real = double, + typename Index = int > class IterativeSolverMonitor : public SolverMonitor { public: -- GitLab From b274befe35ab57c2dce1d6d6dccff5927e781ef7 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Sun, 26 Dec 2021 22:04:34 +0100 Subject: [PATCH 24/54] Added example of iterative solver of linear systems with solver monitor. --- .../Examples/Solvers/Linear/CMakeLists.txt | 1 + ...terativeLinearSolverWithMonitorExample.cpp | 102 ++++++++++++++++++ ...IterativeLinearSolverWithMonitorExample.cu | 1 + .../Solvers/Linear/tutorial_Linear_solvers.md | 10 +- 4 files changed, 113 insertions(+), 1 deletion(-) create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cu diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt index 575542019..6cc75f08e 100644 --- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -1,5 +1,6 @@ set( COMMON_EXAMPLES IterativeLinearSolverExample + IterativeLinearSolverWithMonitorExample ) if( BUILD_CUDA ) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp new file mode 100644 index 000000000..ac3a8f83e --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp @@ -0,0 +1,102 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template< typename Device > +void iterativeLinearSolverExample() +{ + /*** + * Set the following matrix (dots represent zero matrix elements): + * + * / 2 -1 . . . \ + * | -1 2 -1 . . | + * | . -1 2 -1. . | + * | . . -1 2 -1 | + * \ . . . -1 2 / + */ + 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 ); + + + /*** + * Setup monitor of the iterative solver + */ + using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >; + IterativeSolverMonitorType monitor; + TNL::Solvers::SolverMonitorThread t(monitor); + monitor.setRefreshRate(100); // refresh rate = timeout in milliseconds + monitor.setVerbose(1); + monitor.setStage( "Jacobi stage:" ); + TNL::Timer timer; + monitor.setTimer( timer ); + solver.setSolverMonitor(monitor); + solver.setOmega( 0.001 ); + //std::this_thread::sleep_for(std::chrono::milliseconds(3000)); // wait for 3 seconds to let the monitor doing something + timer.start(); + 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 000000000..c0474305b --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithMonitorExample.cpp \ No newline at end of file diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 2df2f7319..24516565c 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -35,4 +35,12 @@ All iterative solvers for linear systems can be found in the namespace \ref TNL: The result looks as follows: -\include IterativeLinearSolverExample.out \ No newline at end of file +\include IterativeLinearSolverExample.out + + + +\includelineno Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp + +The result looks as follows: + +\include IterativeLinearSolverWithMonitorExample.out \ No newline at end of file -- GitLab From 7c9c0e27c9ea37746e7abbb5044dae2923d671c7 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Mon, 27 Dec 2021 11:28:24 +0100 Subject: [PATCH 25/54] Small fixes in the Jacobi linear solver. --- src/TNL/Solvers/Linear/Jacobi.h | 16 +++++++++++++++- src/TNL/Solvers/Linear/Jacobi.hpp | 20 +++++++++++++++++++- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index cebc84472..1b940355f 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -62,7 +62,7 @@ class Jacobi * defines the following: * * \e jacobi-omega - relaxation parameter of the weighted/damped Jacobi method - 1.0 by default. - * \e residue-period - says after how many iterations the reside is recomputed - 4 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. @@ -91,6 +91,20 @@ class Jacobi */ 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. * diff --git a/src/TNL/Solvers/Linear/Jacobi.hpp b/src/TNL/Solvers/Linear/Jacobi.hpp index 9a900fb6f..ba7ef968c 100644 --- a/src/TNL/Solvers/Linear/Jacobi.hpp +++ b/src/TNL/Solvers/Linear/Jacobi.hpp @@ -24,7 +24,7 @@ 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", "Says after how many iterations the reside is recomputed.", 4 ); + config.addEntry< int >( prefix + "residue-period", "Number of iterations between subsequent recomputations of the residue.", 4 ); } template< typename Matrix > @@ -38,6 +38,8 @@ setup( const Config::ParameterContainer& parameters, const String& prefix ) { 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 ); } @@ -57,6 +59,22 @@ 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 >:: -- GitLab From 6bcf11d2b60246429ffd3d6a4cf1bb428f940245 Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Mon, 27 Dec 2021 11:29:01 +0100 Subject: [PATCH 26/54] Reimplementing SOR solver - should be working even on GPU now. --- src/TNL/Solvers/Linear/SOR.h | 24 ++++++++++++++ src/TNL/Solvers/Linear/SOR.hpp | 58 ++++++++++++++++++++++++++++++---- 2 files changed, 76 insertions(+), 6 deletions(-) diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index c5fc4c8aa..1517821bc 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -93,6 +93,20 @@ class SOR */ 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. * @@ -106,7 +120,17 @@ class SOR bool solve( ConstVectorViewType b, VectorViewType x ) override; protected: + using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; + + void performIteration( const ConstVectorViewType& b, + const ConstVectorViewType& diagonalView, + VectorViewType& x ) const; + RealType omega = 1.0; + + IndexType residuePeriod = 4; + + VectorType diagonal; }; } // namespace Linear diff --git a/src/TNL/Solvers/Linear/SOR.hpp b/src/TNL/Solvers/Linear/SOR.hpp index c77af2a33..c7580830b 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,26 +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->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 ) ); 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 -- GitLab From fa7bd33f3e718ba3fb8d6540066f23b805554d8f Mon Sep 17 00:00:00 2001 From: Tomas Oberhuber Date: Mon, 27 Dec 2021 11:30:21 +0100 Subject: [PATCH 27/54] Deleting deprecated method performSORIteration from matrices. --- src/TNL/Matrices/DenseMatrix.h | 6 ----- src/TNL/Matrices/DenseMatrix.hpp | 22 ------------------- src/TNL/Matrices/DenseMatrixView.h | 6 ----- src/TNL/Matrices/DenseMatrixView.hpp | 22 ------------------- src/TNL/Matrices/LambdaMatrix.h | 7 ------ src/TNL/Matrices/LambdaMatrix.hpp | 16 -------------- src/TNL/Matrices/MultidiagonalMatrix.h | 7 ------ src/TNL/Matrices/MultidiagonalMatrix.hpp | 23 -------------------- src/TNL/Matrices/MultidiagonalMatrixView.h | 7 ------ src/TNL/Matrices/MultidiagonalMatrixView.hpp | 22 ------------------- src/TNL/Matrices/SparseMatrix.h | 6 ----- src/TNL/Matrices/SparseMatrix.hpp | 19 ---------------- src/TNL/Matrices/SparseMatrixView.h | 6 ----- src/TNL/Matrices/SparseMatrixView.hpp | 17 --------------- src/TNL/Matrices/TridiagonalMatrix.h | 7 ------ src/TNL/Matrices/TridiagonalMatrix.hpp | 21 ------------------ src/TNL/Matrices/TridiagonalMatrixView.h | 7 ------ src/TNL/Matrices/TridiagonalMatrixView.hpp | 22 ------------------- 18 files changed, 243 deletions(-) diff --git a/src/TNL/Matrices/DenseMatrix.h b/src/TNL/Matrices/DenseMatrix.h index 3d0e2d62d..d46b844e6 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 a42421aa7..ae3f7ee6b 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 89a2219b3..712a565ae 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 3a44269d1..336fcbb9a 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/LambdaMatrix.h b/src/TNL/Matrices/LambdaMatrix.h index 2b788ed52..0278a00ef 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 730c8b5b1..e348a6dba 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 3a74b0b44..d56f1fd06 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 244188831..6ece15b4c 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 869ddb973..a0e13874d 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 550158e4d..65293b809 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 fa3441b01..76d47b50d 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 93cd9b173..b7df383b0 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 eaf692ac6..704d8c7ed 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 d7ed73815..3d0ce7342 100644 --- a/src/TNL/Matrices/SparseMatrixView.hpp +++ b/src/TNL/Matrices/SparseMatrixView.hpp @@ -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 45e6d132f..bb47f6812 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 93ebcd8b3..80236b7de 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 c8e0ecdca..3420f6815 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 c67073381..64da40624 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, -- GitLab From b9fb9066d5162f5646b1877477e5a71815de3b9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 17:55:03 +0100 Subject: [PATCH 28/54] Fixing CMakeList in linear solvers tutorial. --- .../Tutorials/Solvers/Linear/CMakeLists.txt | 33 ------------------- 1 file changed, 33 deletions(-) diff --git a/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt b/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt index 71facac2f..e69de29bb 100644 --- a/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt +++ b/Documentation/Tutorials/Solvers/Linear/CMakeLists.txt @@ -1,33 +0,0 @@ -IF( BUILD_CUDA ) - CUDA_ADD_EXECUTABLE( ArrayAllocation ArrayAllocation.cu ) - ADD_CUSTOM_COMMAND( COMMAND ArrayAllocation > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayAllocation.out OUTPUT ArrayAllocation.out ) - CUDA_ADD_EXECUTABLE( ArrayIO ArrayIO.cu ) - ADD_CUSTOM_COMMAND( COMMAND ArrayIO > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayIO.out OUTPUT ArrayIO.out ) - CUDA_ADD_EXECUTABLE( ArrayView-1 ArrayView-1.cu ) - ADD_CUSTOM_COMMAND( COMMAND ArrayView-1 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayView-1.out OUTPUT ArrayView-1.out ) - CUDA_ADD_EXECUTABLE( ArrayView-2 ArrayView-2.cu ) - ADD_CUSTOM_COMMAND( COMMAND ArrayView-2 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayView-2.out OUTPUT ArrayView-2.out ) - CUDA_ADD_EXECUTABLE( ArrayViewForElements ArrayViewForElements.cu ) - ADD_CUSTOM_COMMAND( COMMAND ArrayViewForElements > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ArrayViewForElements.out OUTPUT ArrayViewForElements.out ) - CUDA_ADD_EXECUTABLE( contains contains.cu ) - ADD_CUSTOM_COMMAND( COMMAND contains > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/contains.out OUTPUT contains.out ) - CUDA_ADD_EXECUTABLE( ElementsAccessing-1 ElementsAccessing-1.cu ) - ADD_CUSTOM_COMMAND( COMMAND ElementsAccessing-1 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ElementsAccessing-1.out OUTPUT ElementsAccessing-1.out ) - CUDA_ADD_EXECUTABLE( ElementsAccessing-2 ElementsAccessing-2.cu ) - ADD_CUSTOM_COMMAND( COMMAND ElementsAccessing-2 > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/ElementsAccessing-2.out OUTPUT ElementsAccessing-2.out ) - ADD_EXECUTABLE( StaticArrayExample StaticArrayExample.cpp ) - ADD_CUSTOM_COMMAND( COMMAND StaticArrayExample > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticArrayExample.out OUTPUT StaticArrayExample.out ) -ENDIF() - -IF( BUILD_CUDA ) -ADD_CUSTOM_TARGET( TutorialsArrays-cuda ALL DEPENDS - ArrayAllocation.out - ArrayIO.out - ArrayView-1.out - ArrayView-2.out - contains.out - ElementsAccessing-1.out - ElementsAccessing-2.out - ArrayViewForElements.out - StaticArrayExample.out ) -ENDIF() -- GitLab From a78c169b0f9f831b4c5f63f66e03e83eea679b9e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 17:55:49 +0100 Subject: [PATCH 29/54] Removing performSORIteration from PyTNL matrix interaface. --- src/Python/pytnl/tnl/SparseMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Python/pytnl/tnl/SparseMatrix.h b/src/Python/pytnl/tnl/SparseMatrix.h index aa0ea3394..a8bfc3e9d 100644 --- a/src/Python/pytnl/tnl/SparseMatrix.h +++ b/src/Python/pytnl/tnl/SparseMatrix.h @@ -149,7 +149,7 @@ 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 >) + //.def("performSORIteration", &Matrix::template performSORIteration< VectorType, VectorType >) // TODO: export for more types .def("assign", []( Matrix& matrix, const Matrix& other ) -> Matrix& { return matrix = other; -- GitLab From 1e0dc11bbfecd8e2a8c2517db88ff4e85518a8f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 17:56:13 +0100 Subject: [PATCH 30/54] Editing documentation of GMRES. --- src/TNL/Solvers/Linear/GMRES.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TNL/Solvers/Linear/GMRES.h b/src/TNL/Solvers/Linear/GMRES.h index 7ff51cbb9..c2f0c32ef 100644 --- a/src/TNL/Solvers/Linear/GMRES.h +++ b/src/TNL/Solvers/Linear/GMRES.h @@ -36,7 +36,7 @@ namespace TNL { * \e CWY - Compact WY form of the Householder reflections). * * See [Wikipedia](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method) for more details. - * + * * \tparam Matrix is type of matrix describing the linear system. */ template< typename Matrix > -- GitLab From 86ae79d4aff285e877b80e894f3c8a39377eab51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 20:04:50 +0100 Subject: [PATCH 31/54] Fixes of Jacobi and SOR solver to be acepted by nvcc. --- src/TNL/Solvers/Linear/Jacobi.h | 11 ++++++----- src/TNL/Solvers/Linear/SOR.h | 10 ++++++---- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index 1b940355f..18afea204 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -121,16 +121,17 @@ class Jacobi using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; - void performIteration( const ConstVectorViewType& b, - const ConstVectorViewType& diagonalView, - const ConstVectorViewType& in, - VectorViewType& out ) const; - 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 diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index 1517821bc..a64a4a509 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -122,15 +122,17 @@ class SOR protected: using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >; - void performIteration( const ConstVectorViewType& b, - const ConstVectorViewType& diagonalView, - VectorViewType& x ) const; - 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; + }; } // namespace Linear -- GitLab From bb3324240a0e50d8e8cecfa34d7aa4419954af0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 20:36:51 +0100 Subject: [PATCH 32/54] Writting documentation on solver monitor. --- src/TNL/Solvers/SolverMonitor.h | 203 +++++++++++++++++++------------- 1 file changed, 119 insertions(+), 84 deletions(-) diff --git a/src/TNL/Solvers/SolverMonitor.h b/src/TNL/Solvers/SolverMonitor.h index ff349224b..4f25e862e 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 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 mian 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 -- GitLab From 36f368c73898d93db44a9f63774ff8310dfa2fcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Mon, 27 Dec 2021 20:48:57 +0100 Subject: [PATCH 33/54] Iterative solver monitor does not print elpased time, time and tau if one does not set timer and time step respectively. --- src/TNL/Solvers/IterativeSolverMonitor.hpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/TNL/Solvers/IterativeSolverMonitor.hpp b/src/TNL/Solvers/IterativeSolverMonitor.hpp index a1fa8c454..996f47f95 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; -- GitLab From d82fd0bdfb15aa1aac0dd7c1a06cd3a1666a55d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Tue, 28 Dec 2021 12:59:00 +0100 Subject: [PATCH 34/54] Deleting performSORIteration method from the distributed matrix. --- src/TNL/Matrices/DistributedMatrix.h | 7 ------- src/TNL/Matrices/DistributedMatrix_impl.h | 12 ------------ 2 files changed, 19 deletions(-) diff --git a/src/TNL/Matrices/DistributedMatrix.h b/src/TNL/Matrices/DistributedMatrix.h index 250521f8e..cdc6387cc 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 498868509..e66e86ad0 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 -- GitLab From 59a4d824d4029e4878998b4bdd2850b7ff31491a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 13:54:18 +0100 Subject: [PATCH 35/54] Added sequential device to atomic operations. --- src/TNL/Algorithms/AtomicOperations.h | 15 +++++++++++++++ src/TNL/Atomic.h | 12 ++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/TNL/Algorithms/AtomicOperations.h b/src/TNL/Algorithms/AtomicOperations.h index dbc1d18a7..4833e33ce 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 0d78a660f..0d0318778 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,17 @@ public: } }; +template< typename T > +class Atomic< T, Devices::Sequential > : public Atomic< T, Devices::Host > +{ + public: + + using Atomic; + using operator=; + using fetch_max; + using fetch_min; +}; + template< typename T > class Atomic< T, Devices::Cuda > { -- GitLab From e97cba150f066ebaf7087714b074553bd590f864 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 13:55:13 +0100 Subject: [PATCH 36/54] Small fixes and typo fix in documentation of solver monitor. --- src/TNL/Solvers/Linear/LinearSolver.h | 1 - src/TNL/Solvers/SolverMonitor.h | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index 4f6dc22b9..8408be911 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -163,7 +163,6 @@ class LinearSolver MatrixPointer matrix = nullptr; PreconditionerPointer preconditioner = nullptr; }; - } // namespace Linear } // namespace Solvers } // namespace TNL diff --git a/src/TNL/Solvers/SolverMonitor.h b/src/TNL/Solvers/SolverMonitor.h index 4f25e862e..fca151cec 100644 --- a/src/TNL/Solvers/SolverMonitor.h +++ b/src/TNL/Solvers/SolverMonitor.h @@ -93,7 +93,7 @@ class SolverMonitor } /** - * \brief Stops the main loop of the monitor. See \ref SolverMonitor::runMainLoop. + * \brief Stops the main loop of the monitor. See \ref TNL::SolverMonitor::runMainLoop. */ void stopMainLoop() { stopped = true; } @@ -101,7 +101,7 @@ class SolverMonitor * \brief Checks whether the main loop was stopped. * * \return true if the main loop was stopped. - * \return false if the mian loop was not stopped yet. + * \return false if the main loop was not stopped yet. */ bool isStopped() const { return stopped; } -- GitLab From 8a3cf685de6eeec2023216dc5400d93c1d0bd991 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 13:59:12 +0100 Subject: [PATCH 37/54] Added TODO to iterative solver monitor. --- src/TNL/Solvers/IterativeSolverMonitor.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h index 401e34c23..651f38cde 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.h +++ b/src/TNL/Solvers/IterativeSolverMonitor.h @@ -115,6 +115,7 @@ class IterativeSolverMonitor : public SolverMonitor IndexType iterations, saved_iterations, iterations_before_refresh; + // TODO: move verbose to SolverMonitor IndexType verbose; IndexType nodesPerIteration; -- GitLab From 50b88716b7dcd58ab3409a3fa20a7c10e1cf50bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 14:00:06 +0100 Subject: [PATCH 38/54] Small changes in example on linear iterative solvers with monitor. --- .../Linear/IterativeLinearSolverExample.cpp | 6 ++--- ...terativeLinearSolverWithMonitorExample.cpp | 23 ++++++++----------- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp index 229aab369..fe4d42b14 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp @@ -4,7 +4,7 @@ #include #include #include -#include +#include template< typename Device > void iterativeLinearSolverExample() @@ -63,7 +63,7 @@ void iterativeLinearSolverExample() /*** * Solve the linear system */ - using LinearSolver = TNL::Solvers::Linear::GMRES< MatrixType >; + using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >; LinearSolver solver; solver.setMatrix( matrix_ptr ); solver.solve( b, x ); @@ -73,7 +73,7 @@ void iterativeLinearSolverExample() int main( int argc, char* argv[] ) { std::cout << "Solving linear system on host: " << std::endl; - iterativeLinearSolverExample< TNL::Devices::Host >(); + iterativeLinearSolverExample< TNL::Devices::Sequential >(); #ifdef HAVE_CUDA std::cout << "Solving linear system on CUDA device: " << std::endl; diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp index ac3a8f83e..a37da76bb 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp @@ -4,7 +4,7 @@ #include #include #include -#include +#include #include #include @@ -14,11 +14,11 @@ void iterativeLinearSolverExample() /*** * Set the following matrix (dots represent zero matrix elements): * - * / 2 -1 . . . \ - * | -1 2 -1 . . | - * | . -1 2 -1. . | - * | . . -1 2 -1 | - * \ . . . -1 2 / + * / 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 >; @@ -68,23 +68,18 @@ void iterativeLinearSolverExample() 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 t(monitor); - monitor.setRefreshRate(100); // refresh rate = timeout in milliseconds + TNL::Solvers::SolverMonitorThread monitorThread(monitor); + monitor.setRefreshRate(10); // refresh rate in milliseconds monitor.setVerbose(1); monitor.setStage( "Jacobi stage:" ); - TNL::Timer timer; - monitor.setTimer( timer ); solver.setSolverMonitor(monitor); - solver.setOmega( 0.001 ); - //std::this_thread::sleep_for(std::chrono::milliseconds(3000)); // wait for 3 seconds to let the monitor doing something - timer.start(); solver.solve( b, x ); monitor.stopMainLoop(); std::cout << "Vector x = " << x << std::endl; -- GitLab From 29bbaef605830d75bdc0b872cc059de3cc681ac3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 14:00:35 +0100 Subject: [PATCH 39/54] Added example on linear iterative solvers with monitor and timer. --- .../Examples/Solvers/Linear/CMakeLists.txt | 1 + .../IterativeLinearSolverWithTimerExample.cpp | 101 ++++++++++++++++++ .../IterativeLinearSolverWithTimerExample.cu | 1 + 3 files changed, 103 insertions(+) create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cu diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt index 6cc75f08e..a07ba232a 100644 --- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -1,6 +1,7 @@ set( COMMON_EXAMPLES IterativeLinearSolverExample IterativeLinearSolverWithMonitorExample + IterativeLinearSolverWithTimerExample ) if( BUILD_CUDA ) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp new file mode 100644 index 000000000..278ccbad0 --- /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 000000000..4c497154d --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithTimerExample.cpp \ No newline at end of file -- GitLab From 0400f1583ce05d68ba967b7f8cd2cea5a27ff680 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 14:01:08 +0100 Subject: [PATCH 40/54] Writting tutorial on iterative linear solvers with monitor and timer. --- .../Solvers/Linear/tutorial_Linear_solvers.md | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 24516565c..96f689e99 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -22,7 +22,7 @@ Solvers of linear systems are one of the most important algorithms in scientific 4. Modified Gramm-Schmidt with reorthogonalization, MGSR 5. Compact WY form of the Householder reflections, CWY -The Krylov subspace methods can be combined with the following precoditioners: +The Krylov subspace methods can be combined with the following preconditioners: 1. Jacobi 2. ILU - CPU only currently @@ -33,14 +33,44 @@ All iterative solvers for linear systems can be found in the namespace \ref TNL: \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 - +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 IterativeLinearSolverWithMonitorExample.out \ No newline at end of file +\include IterativeLinearSolverWithTimerExample.out -- GitLab From 099ae2fba3d2e3512043add7b41ca74541ad1317 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= Date: Tue, 28 Dec 2021 09:18:32 +0100 Subject: [PATCH 41/54] Removed performSORIteration from DistributedMatrix --- src/Python/pytnl/tnl/SparseMatrix.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Python/pytnl/tnl/SparseMatrix.h b/src/Python/pytnl/tnl/SparseMatrix.h index a8bfc3e9d..2c66c6a98 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; -- GitLab From c6746c8ccbc2a35d37188d325eb121facdd31684 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 18:28:36 +0100 Subject: [PATCH 42/54] Fix of using of using in Atomic.h --- src/TNL/Atomic.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/TNL/Atomic.h b/src/TNL/Atomic.h index 0d0318778..ff8d419a1 100644 --- a/src/TNL/Atomic.h +++ b/src/TNL/Atomic.h @@ -100,12 +100,13 @@ public: template< typename T > class Atomic< T, Devices::Sequential > : public Atomic< T, Devices::Host > { + using Base = Atomic< T, Devices::Host >; public: - using Atomic; - using operator=; - using fetch_max; - using fetch_min; + using Base::Atomic; + using Base::operator=; + using Base::fetch_max; + using Base::fetch_min; }; template< typename T > -- GitLab From f02c6d5ba3ea5a5936947a83faf1960cf5f0ebdf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 18:29:12 +0100 Subject: [PATCH 43/54] Fix of typo in LinearSolver. --- src/TNL/Solvers/Linear/LinearSolver.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index 8408be911..858f08189 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -69,7 +69,7 @@ class LinearSolver using ConstVectorViewType = typename Traits< Matrix >::ConstVectorViewType; /** - * \brief Type of the matrix representign the linear system. + * \brief Type of the matrix representing the linear system. */ using MatrixType = Matrix; @@ -89,7 +89,7 @@ class LinearSolver using PreconditionerPointer = std::shared_ptr< std::add_const_t< PreconditionerType > >; /** - * \brief This is method defines configuration entries for setup of the linear iterative solver. + * \brief This method defines configuration entries for setup of the linear iterative solver. * * See \ref IterativeSolver::configSetup. * -- GitLab From 097584e45c6865bbbac4a505602b6ecdc81c5a9e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 18:29:50 +0100 Subject: [PATCH 44/54] Writting documentation on preconditioners of linear solvers. --- .../Solvers/Linear/tutorial_Linear_solvers.md | 6 + .../Solvers/Linear/Preconditioners/Diagonal.h | 174 ++++++++++++++---- src/TNL/Solvers/Linear/Preconditioners/ILU0.h | 123 +++++++++---- src/TNL/Solvers/Linear/Preconditioners/ILUT.h | 75 +++++++- .../Linear/Preconditioners/Preconditioner.h | 129 +++++++++---- .../Linear/Preconditioners/_NamespaceDoxy.h | 29 +++ 6 files changed, 418 insertions(+), 118 deletions(-) create mode 100644 src/TNL/Solvers/Linear/Preconditioners/_NamespaceDoxy.h diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 96f689e99..c58b28273 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -29,6 +29,8 @@ The Krylov subspace methods can be combined with the following preconditioners: # 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 @@ -55,6 +57,8 @@ 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 @@ -74,3 +78,5 @@ The only changes happen on lines 83-85 where we create an instance of TNL timer The result looks as follows: \include IterativeLinearSolverWithTimerExample.out + +## Setup with preconditioner \ No newline at end of file diff --git a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h index 22bfceb00..a25bd9040 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h @@ -15,58 +15,160 @@ #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)) + * + * \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: - virtual void update( const MatrixPointer& matrixPointer ) override; + /** + * \brief Floating point type used for computations. + */ + using RealType = typename Matrix::RealType; - virtual void solve( ConstVectorViewType b, VectorViewType x ) const override; + /** + * \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; -protected: - VectorType diagonal; + /** + * \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 the matrix representing the linear system. + */ + using MatrixType = Matrix; + + /** + * \brief Type of shared pointer to the matrix. + */ + using typename Preconditioner< Matrix >::MatrixPointer; + + /** + * \brief This methods update the preconditione 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 methods 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 This methods update the preconditione with respect to given matrix. + * + * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + */ + using typename Preconditioner< MatrixType >::MatrixPointer; + + /** + * \brief This methods update the preconditione 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 methods 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 b623e86fb..4b3586448 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,98 @@ 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). + * + * \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 methods update the preconditione 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 methods 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 88f5441ba..843adaec5 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h @@ -27,18 +27,77 @@ 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) + * + * \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 methods update the preconditione 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 methods 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 > diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h index 672f5167e..43b429da7 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h @@ -19,47 +19,104 @@ #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. */ -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 methods update the preconditione 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 methods 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 000000000..bd29562ff --- /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 -- GitLab From 57b6e2208bfdaddbc2b9070538a76efa7036340b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 18:40:21 +0100 Subject: [PATCH 45/54] Small fixes in comments in examples on iterative solvers for linear systems. --- .../Solvers/Linear/IterativeLinearSolverExample.cpp | 4 ++-- .../Linear/IterativeLinearSolverWithMonitorExample.cpp | 6 +++--- .../Linear/IterativeLinearSolverWithTimerExample.cpp | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp index fe4d42b14..222a8d866 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverExample.cpp @@ -52,7 +52,7 @@ void iterativeLinearSolverExample() std::cout << *matrix_ptr << std::endl; /*** - * Set the right-hand side vector + * Set the right-hand side vector. */ Vector x( size, 1.0 ); Vector b( size ); @@ -61,7 +61,7 @@ void iterativeLinearSolverExample() std::cout << "Vector b = " << b << std::endl; /*** - * Solve the linear system + * Solve the linear system. */ using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >; LinearSolver solver; diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp index a37da76bb..92cf9c840 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp @@ -54,7 +54,7 @@ void iterativeLinearSolverExample() std::cout << *matrix_ptr << std::endl; /*** - * Set the right-hand side vector + * Set the right-hand side vector. */ Vector x( size, 1.0 ); Vector b( size ); @@ -63,7 +63,7 @@ void iterativeLinearSolverExample() std::cout << "Vector b = " << b << std::endl; /*** - * Setup solver of the linear system + * Setup solver of the linear system. */ using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >; LinearSolver solver; @@ -71,7 +71,7 @@ void iterativeLinearSolverExample() solver.setOmega( 0.0005 ); /*** - * Setup monitor of the iterative solver + * Setup monitor of the iterative solver. */ using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >; IterativeSolverMonitorType monitor; diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp index 278ccbad0..a9ae53c1d 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp @@ -55,7 +55,7 @@ void iterativeLinearSolverExample() std::cout << *matrix_ptr << std::endl; /*** - * Set the right-hand side vector + * Set the right-hand side vector. */ Vector x( size, 1.0 ); Vector b( size ); @@ -64,7 +64,7 @@ void iterativeLinearSolverExample() std::cout << "Vector b = " << b << std::endl; /*** - * Setup solver of the linear system + * Setup solver of the linear system. */ using LinearSolver = TNL::Solvers::Linear::Jacobi< MatrixType >; LinearSolver solver; @@ -72,7 +72,7 @@ void iterativeLinearSolverExample() solver.setOmega( 0.0005 ); /*** - * Setup monitor of the iterative solver + * Setup monitor of the iterative solver. */ using IterativeSolverMonitorType = TNL::Solvers::IterativeSolverMonitor< double, int >; IterativeSolverMonitorType monitor; -- GitLab From 8b1748bdff1b6317a76ab4b71acac862e0229b6c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 21:08:38 +0100 Subject: [PATCH 46/54] Writting documentation on linear preconditioners. --- .../Examples/Solvers/Linear/CMakeLists.txt | 1 + ...eLinearSolverWithPreconditionerExample.cpp | 87 +++++++++++++++++++ ...veLinearSolverWithPreconditionerExample.cu | 1 + .../Solvers/Linear/tutorial_Linear_solvers.md | 12 ++- src/TNL/Solvers/IterativeSolverMonitor.h | 18 ++++ src/TNL/Solvers/Linear/LinearSolver.h | 10 +++ src/TNL/Solvers/Linear/Preconditioners/ILU0.h | 1 + .../Linear/Preconditioners/Preconditioner.h | 9 ++ 8 files changed, 138 insertions(+), 1 deletion(-) create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt index a07ba232a..0d2354617 100644 --- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -2,6 +2,7 @@ set( COMMON_EXAMPLES IterativeLinearSolverExample IterativeLinearSolverWithMonitorExample IterativeLinearSolverWithTimerExample + IterativeLinearSolverWithPreconditionerExample ) if( BUILD_CUDA ) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp new file mode 100644 index 000000000..3a2141074 --- /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 000000000..dc6a12a58 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithPreconditionerExample.cpp \ No newline at end of file diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index c58b28273..25e709ae6 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -79,4 +79,14 @@ The result looks as follows: \include IterativeLinearSolverWithTimerExample.out -## Setup with preconditioner \ No newline at end of file +## Setup with preconditioner + +Preconditioners of iterative solvers can significantly improve the performance of the solver. In the case of the linear systems, they are used mainly with the Krylov subspace methods. Preconditioners cannot be used with the starionary methods (\ref TNL::Solvers::Linear::Jacobi and \ref TNL::Solvers::Linear::SOR). The following example shows how to setup an iterative solver of linear systems with preconditioning. + +\includelineno Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp + +In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen on the lines (68-72). Similar to the matrix of the linear system, the preconditioner is passed to the solver by the means of smart shared pointer (\ref std::shared_ptr). The instance is created on the lines 68 and 69. Next we just need to connect the solver with the preconditioner (line 72, \ref TNL::Solvers::Linear::LinearSolver). + +The result looks as follows: + +\include IterativeLinearSolverWithPreconditionerExample.out diff --git a/src/TNL/Solvers/IterativeSolverMonitor.h b/src/TNL/Solvers/IterativeSolverMonitor.h index 651f38cde..8f42a6b28 100644 --- a/src/TNL/Solvers/IterativeSolverMonitor.h +++ b/src/TNL/Solvers/IterativeSolverMonitor.h @@ -20,6 +20,24 @@ namespace TNL { * * \tparam Real is a type of the floating-point arithmetics. * \tparam Index is an indexing type. + * + * The following example shows how to use the iterative solver monitor for monitoring + * convergence of linear iterative solver: + * + * \includelineno Solvers/Linear/IterativeLinearSolverWithMonitorExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverWithMonitorExample.out + * + * The following example shows how to employ timer (\ref TNL::Timer) to the monitor + * of iterative solvers: + * + * \includelineno Solvers/Linear/IterativeLinearSolverWithTimerExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverWithTimerExample.out */ template< typename Real = double, typename Index = int > diff --git a/src/TNL/Solvers/Linear/LinearSolver.h b/src/TNL/Solvers/Linear/LinearSolver.h index 858f08189..f2de9a8f4 100644 --- a/src/TNL/Solvers/Linear/LinearSolver.h +++ b/src/TNL/Solvers/Linear/LinearSolver.h @@ -34,6 +34,16 @@ namespace TNL { * \ref LinearSolver::setPreconditioner. * * \tparam Matrix is type of matrix representing the linear system. + * + * The following example demonstrates the use iterative linear solvers: + * + * \includelineno Solvers/Linear/IterativeLinearSolverExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverExample.out + * + * See also \ref TNL::Solvers::IterativeSolverMonitor for monitoring of iterative solvers. */ template< typename Matrix > class LinearSolver diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h index 4b3586448..6ce3fef51 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h @@ -76,6 +76,7 @@ class ILU0_impl< Matrix, Real, Devices::Host, Index > : public Preconditioner< Matrix > { public: + /** * \brief Floating point type used for computations. */ diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h index 43b429da7..e0dc721e3 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h @@ -27,6 +27,15 @@ namespace TNL { * \brief Base class for preconditioners of of iterative solvers of linear systems. * * \tparam Matrix is type of matrix describing the linear system. + * + * The following example shows how to setup an iterative solver of linear systems with + * preconditioning: + * + * \includelineno Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp + * + * The result looks as follows: + * + * \include IterativeLinearSolverWithPreconditionerExample.out */ template< typename Matrix > class Preconditioner -- GitLab From fd1386cbf5d7b26998faf6fa8d3bccec5936b1e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Wed, 29 Dec 2021 21:24:54 +0100 Subject: [PATCH 47/54] Fixes in tutorial and example on linear preconditioners. --- .../Linear/IterativeLinearSolverWithPreconditionerExample.cpp | 2 +- .../Tutorials/Solvers/Linear/tutorial_Linear_solvers.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp index 3a2141074..035dbb6ce 100644 --- a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp @@ -66,7 +66,7 @@ void iterativeLinearSolverExample() */ using LinearSolver = TNL::Solvers::Linear::TFQMR< MatrixType >; using Preconditioner = TNL::Solvers::Linear::Preconditioners::Diagonal< MatrixType >; - auto preconditioner_ptr = std::make_shared< Preconditioner >; + auto preconditioner_ptr = std::make_shared< Preconditioner >(); preconditioner_ptr->update( matrix_ptr ); LinearSolver solver; solver.setMatrix( matrix_ptr ); diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 25e709ae6..4832f4f12 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -85,7 +85,7 @@ Preconditioners of iterative solvers can significantly improve the performance o \includelineno Solvers/Linear/IterativeLinearSolverWithPreconditionerExample.cpp -In this example, we solve the same problem as in all other examples in this section. The only differences concerning the preconditioner happen on the lines (68-72). Similar to the matrix of the linear system, the preconditioner is passed to the solver by the means of smart shared pointer (\ref std::shared_ptr). The instance is created on the lines 68 and 69. Next we just need to connect the solver with the preconditioner (line 72, \ref TNL::Solvers::Linear::LinearSolver). +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: -- GitLab From f87849867eccfe51ddd2f9d5dbbcaf3498a36915 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:07:50 +0100 Subject: [PATCH 48/54] Added example to documentation of iterative linear solvers. --- src/TNL/Solvers/Linear/BICGStab.h | 2 ++ src/TNL/Solvers/Linear/BICGStabL.h | 4 +++- src/TNL/Solvers/Linear/CG.h | 2 ++ src/TNL/Solvers/Linear/GMRES.h | 4 ++++ src/TNL/Solvers/Linear/Jacobi.h | 2 ++ 5 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/TNL/Solvers/Linear/BICGStab.h b/src/TNL/Solvers/Linear/BICGStab.h index 2d1e7d3a1..c9d973ad1 100644 --- a/src/TNL/Solvers/Linear/BICGStab.h +++ b/src/TNL/Solvers/Linear/BICGStab.h @@ -23,6 +23,8 @@ namespace TNL { * * 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 > diff --git a/src/TNL/Solvers/Linear/BICGStabL.h b/src/TNL/Solvers/Linear/BICGStabL.h index 682d14370..1814ac314 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.h +++ b/src/TNL/Solvers/Linear/BICGStabL.h @@ -71,7 +71,9 @@ namespace TNL { * * [4] Martin H. Gutknecht, "Variants of BiCGStab for matrices with complex * spectrum", IPS Research Report No. 91-14 (1991). - + * + * 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 > diff --git a/src/TNL/Solvers/Linear/CG.h b/src/TNL/Solvers/Linear/CG.h index ecabc0296..605c81fb3 100644 --- a/src/TNL/Solvers/Linear/CG.h +++ b/src/TNL/Solvers/Linear/CG.h @@ -24,6 +24,8 @@ namespace TNL { * 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 diff --git a/src/TNL/Solvers/Linear/GMRES.h b/src/TNL/Solvers/Linear/GMRES.h index c2f0c32ef..67c5edc16 100644 --- a/src/TNL/Solvers/Linear/GMRES.h +++ b/src/TNL/Solvers/Linear/GMRES.h @@ -37,6 +37,10 @@ namespace TNL { * * 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 > diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h index 18afea204..1cc061da9 100644 --- a/src/TNL/Solvers/Linear/Jacobi.h +++ b/src/TNL/Solvers/Linear/Jacobi.h @@ -21,6 +21,8 @@ namespace TNL { * * 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 > -- GitLab From 860952ee829f07b8080cc0ec9084c9b7caf0fbeb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:08:00 +0100 Subject: [PATCH 49/54] Added example to documentation of iterative linear solvers. --- src/TNL/Solvers/Linear/SOR.h | 2 ++ src/TNL/Solvers/Linear/TFQMR.h | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/TNL/Solvers/Linear/SOR.h b/src/TNL/Solvers/Linear/SOR.h index a64a4a509..76ff1db97 100644 --- a/src/TNL/Solvers/Linear/SOR.h +++ b/src/TNL/Solvers/Linear/SOR.h @@ -21,6 +21,8 @@ namespace TNL { * * 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 > diff --git a/src/TNL/Solvers/Linear/TFQMR.h b/src/TNL/Solvers/Linear/TFQMR.h index e9ec68de5..97a243882 100644 --- a/src/TNL/Solvers/Linear/TFQMR.h +++ b/src/TNL/Solvers/Linear/TFQMR.h @@ -21,6 +21,8 @@ namespace TNL { * * 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 > -- GitLab From fdcb91b7069633d6963e3d114cf3a981da95ee01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:08:30 +0100 Subject: [PATCH 50/54] Added example to documentation of linear preconditioners. --- src/TNL/Solvers/Linear/Preconditioners/Diagonal.h | 4 +++- src/TNL/Solvers/Linear/Preconditioners/ILU0.h | 2 ++ src/TNL/Solvers/Linear/Preconditioners/ILUT.h | 8 ++++++++ 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h index a25bd9040..2e2762406 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h @@ -22,7 +22,9 @@ namespace TNL { /** * \brief Diagonal (Jacobi) preconditioner for iterative solvers of linear systems. * - * See [detailed description]([Netlib](http://netlib.org/linalg/html_templates/node55.html)) + * 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. */ diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h index 6ce3fef51..d38357c0e 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h @@ -69,6 +69,8 @@ class ILU0 * * 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 > diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h index 843adaec5..f20600389 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h @@ -32,6 +32,8 @@ class ILUT_impl * * 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 > @@ -147,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 > -- GitLab From 00103113c71bc68bc0e53c014f95dfb0c02d4300 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:09:22 +0100 Subject: [PATCH 51/54] Added example on getLinearSolver and getPreconditioner. --- .../Examples/Solvers/Linear/CMakeLists.txt | 1 + ...iveLinearSolverWithRuntimeTypesExample.cpp | 84 +++++++++++++++++++ ...tiveLinearSolverWithRuntimeTypesExample.cu | 1 + 3 files changed, 86 insertions(+) create mode 100644 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp create mode 120000 Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu diff --git a/Documentation/Examples/Solvers/Linear/CMakeLists.txt b/Documentation/Examples/Solvers/Linear/CMakeLists.txt index 0d2354617..97c77e731 100644 --- a/Documentation/Examples/Solvers/Linear/CMakeLists.txt +++ b/Documentation/Examples/Solvers/Linear/CMakeLists.txt @@ -3,6 +3,7 @@ set( COMMON_EXAMPLES IterativeLinearSolverWithMonitorExample IterativeLinearSolverWithTimerExample IterativeLinearSolverWithPreconditionerExample + IterativeLinearSolverWithRuntimeTypesExample ) if( BUILD_CUDA ) diff --git a/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cpp new file mode 100644 index 000000000..cbe7707cf --- /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 000000000..ab40ffd37 --- /dev/null +++ b/Documentation/Examples/Solvers/Linear/IterativeLinearSolverWithRuntimeTypesExample.cu @@ -0,0 +1 @@ +IterativeLinearSolverWithRuntimeTypesExample.cpp \ No newline at end of file -- GitLab From 558a6a99d11f5015f91f784739b180348c4467af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:09:56 +0100 Subject: [PATCH 52/54] Writting documentation on getLinearSolver and getPreconditioner. --- src/TNL/Solvers/LinearSolverTypeResolver.h | 57 +++++++++++++++++++++- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/src/TNL/Solvers/LinearSolverTypeResolver.h b/src/TNL/Solvers/LinearSolverTypeResolver.h index 5615ea707..db6fb5121 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 -- GitLab From 2be4c35aad069ab73b843ff415f5c777ecf93349 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 13:10:16 +0100 Subject: [PATCH 53/54] Writting tutorial on getLinearSolver and getPreconditioner. --- .../Solvers/Linear/tutorial_Linear_solvers.md | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 4832f4f12..0c2723aee 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -2,7 +2,6 @@ [TOC] - # Introduction Solvers of linear systems are one of the most important algorithms in scientific computations. TNL offers the followiing iterative methods: @@ -22,10 +21,12 @@ Solvers of linear systems are one of the most important algorithms in scientific 4. Modified Gramm-Schmidt with reorthogonalization, MGSR 5. Compact WY form of the Householder reflections, CWY -The Krylov subspace methods can be combined with the following preconditioners: +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. Jacobi -2. ILU - CPU only currently +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 @@ -90,3 +91,15 @@ In this example, we solve the same problem as in all other examples in this sect 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 -- GitLab From 449d8097cf969fc8558fd3a209c5c0f06eb7ed4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= Date: Thu, 30 Dec 2021 16:16:24 +0100 Subject: [PATCH 54/54] Small fixes in documentation of lineat solvers. --- .../Solvers/Linear/tutorial_Linear_solvers.md | 2 +- src/TNL/Solvers/Linear/BICGStabL.h | 23 +------------------ .../Solvers/Linear/Preconditioners/Diagonal.h | 13 +++++------ src/TNL/Solvers/Linear/Preconditioners/ILU0.h | 4 ++-- src/TNL/Solvers/Linear/Preconditioners/ILUT.h | 4 ++-- .../Linear/Preconditioners/Preconditioner.h | 4 ++-- 6 files changed, 14 insertions(+), 36 deletions(-) diff --git a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md index 0c2723aee..b79dbbe1b 100644 --- a/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md +++ b/Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md @@ -8,7 +8,7 @@ Solvers of linear systems are one of the most important algorithms in scientific 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) - CPU only currently + 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) diff --git a/src/TNL/Solvers/Linear/BICGStabL.h b/src/TNL/Solvers/Linear/BICGStabL.h index 1814ac314..44f4f2264 100644 --- a/src/TNL/Solvers/Linear/BICGStabL.h +++ b/src/TNL/Solvers/Linear/BICGStabL.h @@ -9,27 +9,6 @@ /* See Copyright Notice in tnl/Copyright */ /* - * 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 - * proposed by van der Vorst [3]. BiCGstab(1) is equivalent to BiCGstab, and - * 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: * * [5] Gerard L. G. Sleijpen and Henk A. van der Vorst, "Reliable updated @@ -119,7 +98,7 @@ class BICGStabL * 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). * diff --git a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h index 2e2762406..a8775ca5f 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Diagonal.h @@ -14,6 +14,7 @@ #include "Preconditioner.h" + namespace TNL { namespace Solvers { namespace Linear { @@ -72,14 +73,14 @@ class Diagonal using typename Preconditioner< Matrix >::MatrixPointer; /** - * \brief This methods update the preconditione with respect to given matrix. + * \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 methods applies the preconditioner. + * \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. @@ -139,21 +140,19 @@ class Diagonal< Matrices::DistributedMatrix< Matrix > > using typename Preconditioner< MatrixType >::ConstVectorViewType; /** - * \brief This methods update the preconditione with respect to given matrix. - * - * \param matrixPointer smart pointer (\ref std::shared_ptr) to matrix the preconditioner is related to. + * \brief Type of shared pointer to the matrix. */ using typename Preconditioner< MatrixType >::MatrixPointer; /** - * \brief This methods update the preconditione with respect to given matrix. + * \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 methods applies the preconditioner. + * \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. diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h index d38357c0e..a6d128c39 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILU0.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILU0.h @@ -110,14 +110,14 @@ class ILU0_impl< Matrix, Real, Devices::Host, Index > using typename Preconditioner< Matrix >::MatrixPointer; /** - * \brief This methods update the preconditione with respect to given matrix. + * \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 methods applies the preconditioner. + * \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. diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h index f20600389..5e09147f9 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/ILUT.h +++ b/src/TNL/Solvers/Linear/Preconditioners/ILUT.h @@ -87,14 +87,14 @@ class ILUT } /** - * \brief This methods update the preconditione with respect to given matrix. + * \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 methods applies the preconditioner. + * \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. diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h index e0dc721e3..4357dcf63 100644 --- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h +++ b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h @@ -101,7 +101,7 @@ class Preconditioner } /** - * \brief This methods update the preconditione with respect to given matrix. + * \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. */ @@ -109,7 +109,7 @@ class Preconditioner {} /** - * \brief This methods applies the preconditioner. + * \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. -- GitLab