Commit 27635cb9 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Fixed CG and LinearResidueGetter for distributed matrices

parent cdb5807d

src/TNL/Solvers/Linear/CG.h

+3 −2
Original line number Original line Diff line number Diff line
@@ -29,13 +29,14 @@ public: 
   using IndexType = typename Base::IndexType;

   using IndexType = typename Base::IndexType;

   using VectorViewType = typename Base::VectorViewType;

   using VectorViewType = typename Base::VectorViewType;

   using ConstVectorViewType = typename Base::ConstVectorViewType;

   using ConstVectorViewType = typename Base::ConstVectorViewType;

   using VectorType = typename Traits< Matrix >::VectorType;





   bool solve( ConstVectorViewType b, VectorViewType x ) override;

   bool solve( ConstVectorViewType b, VectorViewType x ) override;





protected:

protected:

   void setSize( IndexType size );

   void setSize( const VectorViewType& x );





   Containers::Vector< RealType, DeviceType, IndexType >  r, p, Ap, z;

   VectorType r, p, Ap, z;

};

};





} // namespace Linear

} // namespace Linear

src/TNL/Solvers/Linear/CG_impl.h

+6 −6
Original line number Original line Diff line number Diff line
@@ -21,7 +21,7 @@ bool 
CG< Matrix >::

CG< Matrix >::

solve( ConstVectorViewType b, VectorViewType x )

solve( ConstVectorViewType b, VectorViewType x )

{

{

   this->setSize( this->matrix->getRows() );

   this->setSize( x );

   this->resetIterations();

   this->resetIterations();





   RealType alpha, beta, s1, s2;

   RealType alpha, beta, s1, s2;
@@ -113,12 +113,12 @@ solve( ConstVectorViewType b, VectorViewType x ) 




template< typename Matrix >

template< typename Matrix >

void CG< Matrix >::

void CG< Matrix >::

setSize( IndexType size )

setSize( const VectorViewType& x )

{

{

   r.setSize( size );

   r.setLike( x );

   p.setSize( size );

   p.setLike( x );

   Ap.setSize( size );

   Ap.setLike( x );

   z.setSize( size );

   z.setLike( x );

}

}





} // namespace Linear

} // namespace Linear

src/TNL/Solvers/Linear/LinearResidueGetter_impl.h

+4 −1
Original line number Original line Diff line number Diff line
@@ -11,6 +11,7 @@ 
#pragma once

#pragma once





#include <TNL/Solvers/Linear/LinearResidueGetter.h>

#include <TNL/Solvers/Linear/LinearResidueGetter.h>

#include <TNL/Solvers/Linear/Traits.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Containers/Vector.h>





namespace TNL {

namespace TNL {
@@ -28,10 +29,12 @@ getResidue( const Matrix& matrix, 
   using RealType = typename Matrix::RealType;

   using RealType = typename Matrix::RealType;

   using DeviceType = typename Matrix::DeviceType;

   using DeviceType = typename Matrix::DeviceType;

   using IndexType = typename Matrix::IndexType;

   using IndexType = typename Matrix::IndexType;

   using VectorType = typename Traits< Matrix >::VectorType;





   if( bNorm == 0.0 )

   if( bNorm == 0.0 )

      bNorm = lpNorm( b, 2.0 );

      bNorm = lpNorm( b, 2.0 );

   Containers::Vector< RealType, DeviceType, IndexType > v( b.getSize() );

   VectorType v;

   v.setLike( b );

   matrix.vectorProduct( x, v );

   matrix.vectorProduct( x, v );

   return l2Norm( v - b ) / bNorm;

   return l2Norm( v - b ) / bNorm;

}

}

src/TNL/Solvers/Linear/SOR_impl.h

+0 −1
Original line number Original line Diff line number Diff line
@@ -66,7 +66,6 @@ bool SOR< Matrix > :: solve( ConstVectorViewType b, VectorViewType x ) 
   {

   {

      for( IndexType row = 0; row < size; row ++ )

      for( IndexType row = 0; row < size; row ++ )

         this->matrix->performSORIteration( b, row, x, this->getOmega() );

         this->matrix->performSORIteration( b, row, x, this->getOmega() );

      // FIXME: the LinearResidueGetter works only on the host

      this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) );

      this->setResidue( LinearResidueGetter::getResidue( *this->matrix, x, b, bNorm ) );

      this->refreshSolverMonitor();

      this->refreshSolverMonitor();

   }

   }