Writting documentation on preconditioners of linear solvers. (097584e4) · Commits · TNL / tnl-dev

Documentation/Tutorials/Solvers/Linear/tutorial_Linear_solvers.md

+6 −0

Original line number	Diff line number	Diff line
		@@ -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

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

+138 −36

Original line number	Diff line number	Diff line
		@@ -19,48 +19,150 @@ 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:

		/**
		* \brief Floating point type used for computations.
		*/
		using RealType = typename Matrix::RealType;

		/**
		* \brief Device where the preconditioner will run on and auxillary data will alloacted on.
		*
		* See \ref Devices::Host or \ref Devices::Cuda.
		*/
		using DeviceType = typename Matrix::DeviceType;

		/**
		* \brief Type for indexing.
		*/
		using IndexType = typename Matrix::IndexType;

		/**
		* \brief Type for vector view.
		*/
		using typename Preconditioner< Matrix >::VectorViewType;

		/**
		* \brief Type for constant vector view.
		*/
		using typename Preconditioner< Matrix >::ConstVectorViewType;

		/**
		* \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;
		using VectorType = Containers::Vector< RealType, DeviceType, IndexType >;

		/**
		* \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:

		/**
		* \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;
		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 >;

		/**
		* \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;
		};

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

+85 −38

Original line number	Diff line number	Diff line
		@@ -52,26 +52,73 @@ 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:
		/**
		* \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:

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

+67 −8

Original line number	Diff line number	Diff line
		@@ -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 >
		{
		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 >

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

+93 −36

Original line number	Diff line number	Diff line
		@@ -21,41 +21,98 @@
		namespace TNL {
		namespace Solvers {
		namespace Linear {
		/**
		* \brief Namespace for preconditioners for linear system solvers.
		*/
		namespace Preconditioners {

		/**
		* \brief Base class for preconditioners of of iterative solvers of linear systems.
		*
		* \tparam Matrix is type of matrix describing the linear system.
		*/
		template< typename Matrix >
		class 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 = "" )
		{}
		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() {}
		};