diff --git a/src/TNL/Solvers/Linear/Jacobi.h b/src/TNL/Solvers/Linear/Jacobi.h
index 5016d2330cdbb4b8b7a48625f7f80f189d41fc54..0cdf57845b111c143c387ded596ba679f7ad0f3a 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 <TNL/Containers/Vector.h>
+#include <TNL/Solvers/Linear/LinearSolver.h>
#include <TNL/Solvers/Linear/Utils/LinearResidueGetter.h>
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 <TNL/Solvers/Linear/Jacobi.hpp>
\ No newline at end of file
diff --git a/src/TNL/Solvers/Linear/Jacobi.hpp b/src/TNL/Solvers/Linear/Jacobi.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..eee368973d11aa0b72987edd6fc25e4c0482f5f9
--- /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 <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