Reimplementing the Jacobi solver.

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

src/TNL/Solvers/Linear/Jacobi.hpp

@@ -11,6 +11,7 @@
 #include <TNL/Containers/Vector.h>
 #include <TNL/Solvers/Linear/LinearSolver.h>
 #include <TNL/Solvers/Linear/Utils/LinearResidueGetter.h>
+#include <TNL/Functional.h>
 
 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