Splitting Jacobi.h file and writing documentation for the Jacobi method. (9deda633) · Commits · TNL / tnl-dev

src/TNL/Solvers/Linear/Jacobi.h

+90 −68

Original line number	Diff line number	Diff line
		@@ -8,79 +8,99 @@

		#pragma once

		#include "LinearSolver.h"

		#include <TNL/Containers/Vector.h>
		#include <TNL/Solvers/Linear/LinearSolver.h>
		#include <TNL/Solvers/Linear/Utils/LinearResidueGetter.h>

		namespace TNL {
		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:

		/**
		* \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;
		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 );
		}
		/**
		* \brief Type for indexing.
		*/
		using IndexType = typename Base::IndexType;

		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 );
		}
		/**
		* \brief Type for vector view.
		*/
		using VectorViewType = typename Base::VectorViewType;

		void setOmega( RealType omega )
		{
		this->omega = omega;
		}
		/**
		* \brief Type for constant vector view.
		*/
		using ConstVectorViewType = typename Base::ConstVectorViewType;

		RealType getOmega() const
		{
		return omega;
		}
		/**
		* \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 = "" );

		bool solve( ConstVectorViewType b, VectorViewType x ) override
		{
		const IndexType size = this->matrix->getRows();
		/**
		* \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;

		Containers::Vector< RealType, DeviceType, IndexType > aux;
		aux.setSize( size );
		/**
		* \brief Setter of the relaxation parameter.
		*
		* \param omega the relaxation parameter. It is 1 by default.
		*/
		void setOmega( RealType omega );

		this->resetIterations();
		this->setResidue( this->getConvergenceResidue() + 1.0 );
		/**
		* \brief Getter of the relaxation parameter.
		*
		* \return value of the relaxation parameter.
		*/
		RealType getOmega() const;

		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();
		}
		/**
		* \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;
		@@ -89,3 +109,5 @@ protected:
		} // namespace Linear
		} // namespace Solvers
		} // namespace TNL

		#include <TNL/Solvers/Linear/Jacobi.hpp>
		No newline at end of file

src/TNL/Solvers/Linear/Jacobi.hpp

0 → 100644

+87 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		Jacobi.hpp - description
		-------------------
		begin : Dec 22, 2021
		copyright : (C) 2021 by Tomas Oberhuber
		email : tomas.oberhuber@fjfi.cvut.cz
		***************************************************************************/

		#pragma once

		#include <TNL/Containers/Vector.h>
		#include <TNL/Solvers/Linear/LinearSolver.h>
		#include <TNL/Solvers/Linear/Utils/LinearResidueGetter.h>

		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