Commit 2c6da83e authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Optimized smart pointers in linear solvers

parent 3c352e0e

src/TNL/Solvers/Linear/BICGStab.h

+9 −17
Original line number Diff line number Diff line
@@ -12,6 +12,7 @@ 


#include <math.h>

#include <TNL/Object.h>

#include <TNL/SharedPointer.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Containers/SharedVector.h>

#include <TNL/Solvers/Linear/Preconditioners/Dummy.h>
@@ -38,12 +39,9 @@ class BICGStab : public Object, 
   typedef typename Matrix::DeviceType DeviceType;

   typedef Matrix MatrixType;

   typedef Preconditioner PreconditionerType;

   typedef SharedPointer< MatrixType, DeviceType > MatrixPointer;

   // TODO: make this 'typedef SharedPointer< const MatrixType, DeviceType > ConstMatrixPointer;'

   typedef SharedPointer< const MatrixType, DeviceType, true > MatrixPointer;





   public:



   BICGStab();



   String getType() const;
@@ -54,19 +52,13 @@ class BICGStab : public Object, 
   bool setup( const Config::ParameterContainer& parameters,

               const String& prefix = "" );



   void setMatrix( MatrixPointer& matrix );

   void setMatrix( MatrixPointer matrix );



   void setPreconditioner( const PreconditionerType& preconditioner );



#ifdef HAVE_NOT_CXX11

   template< typename VectorPointer,

             typename ResidueGetter  >

   bool solve( const VectorPointer& b, VectorPointer& x );

#else

   template< typename VectorPointer,

             typename ResidueGetter = LinearResidueGetter< MatrixPointer, VectorPointer >  >

   bool solve( const VectorPointer& b, VectorPointer& x );

#endif

   template< typename Vector,

             typename ResidueGetter = LinearResidueGetter< Matrix, Vector >  >

   bool solve( const Vector& b, Vector& x );



   ~BICGStab();

src/TNL/Solvers/Linear/BICGStab_impl.h

+11 −11
Original line number Diff line number Diff line
@@ -60,7 +60,7 @@ setup( const Config::ParameterContainer& parameters, 


template< typename Matrix,

          typename Preconditioner >

void BICGStab< Matrix, Preconditioner >::setMatrix( MatrixPointer& matrix )

void BICGStab< Matrix, Preconditioner >::setMatrix( MatrixPointer matrix )

{

   this->matrix = matrix;

}
@@ -74,8 +74,8 @@ void BICGStab< Matrix, Preconditioner > :: setPreconditioner( const Precondition 


template< typename Matrix,

          typename Preconditioner >

   template< typename VectorPointer, typename ResidueGetter >

bool BICGStab< Matrix, Preconditioner >::solve( const VectorPointer& b, VectorPointer& x )

   template< typename Vector, typename ResidueGetter >

bool BICGStab< Matrix, Preconditioner >::solve( const Vector& b, Vector& x )

{

   if( ! this->setSize( matrix->getRows() ) ) return false;
@@ -86,7 +86,7 @@ bool BICGStab< Matrix, Preconditioner >::solve( const VectorPointer& b, VectorPo 
   // r_0 = b - A x_0, p_0 = r_0

   // r^ast_0 = r_0



   this->matrix -> vectorProduct( *x, r );

   this->matrix -> vectorProduct( x, r );



   //if( bNorm == 0.0 ) bNorm = 1.0;
@@ -107,10 +107,10 @@ bool BICGStab< Matrix, Preconditioner >::solve( const VectorPointer& b, VectorPo 
   }

   else*/

   {

      r.addVector( *b, 1.0, -1.0 );

      r.addVector( b, 1.0, -1.0 );

      p = r_ast = r;

      rho = r.scalarProduct( r_ast );

      bNorm = b->lpNorm( 2.0 );

      bNorm = b.lpNorm( 2.0 );

   }



   while( this->nextIteration() )
@@ -156,7 +156,7 @@ bool BICGStab< Matrix, Preconditioner >::solve( const VectorPointer& b, VectorPo 
      /****

       * x_{j+1} = x_j + alpha_j * p_j + omega_j * s_j

       */

      x->addVectors( p, alpha, s, omega );

      x.addVectors( p, alpha, s, omega );

      

      /****

       * r_{j+1} = s_j - omega_j * A * s_j
@@ -182,18 +182,18 @@ bool BICGStab< Matrix, Preconditioner >::solve( const VectorPointer& b, VectorPo 
      residue /= bNorm;

      this->setResidue( residue );

      if( this->getIterations() % 10 == 0 )

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

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

   }

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

   this->refreshSolverMonitor( true );

   return this->checkConvergence();

};

}



template< typename Matrix,

          typename Preconditioner >

BICGStab< Matrix, Preconditioner > :: ~BICGStab()

{

};

}



template< typename Matrix,

          typename Preconditioner >
@@ -212,7 +212,7 @@ bool BICGStab< Matrix, Preconditioner > :: setSize( IndexType size ) 
   }

   return true;



};

}



} // namespace Linear

} // namespace Solvers

src/TNL/Solvers/Linear/CG.h

+8 −14
Original line number Diff line number Diff line
@@ -12,6 +12,7 @@ 


#include <math.h>

#include <TNL/Object.h>

#include <TNL/SharedPointer.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Containers/SharedVector.h>

#include <TNL/Solvers/Linear/Preconditioners/Dummy.h>
@@ -37,9 +38,8 @@ class CG : public Object, 
   typedef typename Matrix::DeviceType DeviceType;

   typedef Matrix MatrixType;

   typedef Preconditioner PreconditionerType;

   typedef SharedPointer< MatrixType, DeviceType > MatrixPointer;

   typedef SharedPointer< PreconditionerType, DeviceType > PreconditionerPointer;

   // TODO: make this const

   typedef SharedPointer< const MatrixType, DeviceType, true > MatrixPointer;

   typedef SharedPointer< const PreconditionerType, DeviceType, true > PreconditionerPointer;





   CG();
@@ -52,19 +52,13 @@ class CG : public Object, 
   bool setup( const Config::ParameterContainer& parameters,

               const String& prefix = "" );



   void setMatrix( MatrixPointer& matrix );

   void setMatrix( MatrixPointer matrix );



   void setPreconditioner( const PreconditionerType& preconditioner );



#ifdef HAVE_NOT_CXX11

   template< typename Vector,

             typename ResidueGetter >

             typename ResidueGetter = LinearResidueGetter< Matrix, Vector >  >

   bool solve( const Vector& b, Vector& x );

#else

   template< typename VectorPointer,

             typename ResidueGetter = LinearResidueGetter< Matrix, VectorPointer >  >

   bool solve( const VectorPointer& b, VectorPointer& x );

#endif



   ~CG();

src/TNL/Solvers/Linear/CG_impl.h

+9 −9
Original line number Diff line number Diff line
@@ -52,7 +52,7 @@ setup( const Config::ParameterContainer& parameters, 


template< typename Matrix,

          typename Preconditioner >

void CG< Matrix, Preconditioner >::setMatrix( MatrixPointer& matrix )

void CG< Matrix, Preconditioner >::setMatrix( MatrixPointer matrix )

{

   this->matrix = matrix;

}
@@ -66,10 +66,10 @@ void CG< Matrix, Preconditioner > :: setPreconditioner( const Preconditioner& pr 


template< typename Matrix,

          typename Preconditioner >

   template< typename VectorPointer, typename ResidueGetter >

   template< typename Vector, typename ResidueGetter >

bool

CG< Matrix, Preconditioner >::

solve( const VectorPointer& b, VectorPointer& x )

solve( const Vector& b, Vector& x )

{

   if( ! this->setSize( matrix->getRows() ) ) return false;
@@ -77,13 +77,13 @@ solve( const VectorPointer& b, VectorPointer& x ) 
   this->setResidue( this->getConvergenceResidue() + 1.0 );



   RealType alpha, beta, s1, s2;

   RealType bNorm = b->lpNorm( ( RealType ) 2.0 );

   RealType bNorm = b.lpNorm( ( RealType ) 2.0 );



   /****

    * r_0 = b - A x_0, p_0 = r_0

    */

   this->matrix->vectorProduct( *x, r );

   r. addVector( *b, 1.0, -1.0 );

   this->matrix->vectorProduct( x, r );

   r. addVector( b, 1.0, -1.0 );

   p = r;



   while( this->nextIteration() )
@@ -105,7 +105,7 @@ solve( const VectorPointer& b, VectorPointer& x ) 
      /****

       * 2. x_{j+1} = x_j + \alpha_j p_j

       */

      x->addVector( p, alpha );

      x.addVector( p, alpha );

      

      /****

       * 3. r_{j+1} = r_j - \alpha_j A * p_j
@@ -136,9 +136,9 @@ solve( const VectorPointer& b, VectorPointer& x ) 
      new_r.swap( r );

 

      if( this->getIterations() % 10 == 0 )

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

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

   }

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

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

   this->refreshSolverMonitor( true );

   return this->checkConvergence();

};

src/TNL/Solvers/Linear/GMRES.h

+6 −13
Original line number Diff line number Diff line
@@ -12,12 +12,12 @@ 


#include <math.h>

#include <TNL/Object.h>

#include <TNL/SharedPointer.h>

#include <TNL/Containers/Vector.h>

#include <TNL/Containers/SharedVector.h>

#include <TNL/Solvers/Linear/Preconditioners/Dummy.h>

#include <TNL/Solvers/IterativeSolver.h>

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

#include <TNL/SharedPointer.h>



namespace TNL {

namespace Solvers {
@@ -38,8 +38,7 @@ class GMRES : public Object, 
   typedef typename Matrix :: DeviceType DeviceType;

   typedef Matrix MatrixType;

   typedef Preconditioner PreconditionerType;

   typedef SharedPointer< MatrixType, DeviceType > MatrixPointer;

   // TODO: make this: typedef SharedPointer< const MatrixType, DeviceType > ConstMatrixPointer;

   typedef SharedPointer< const MatrixType, DeviceType, true > MatrixPointer;



   GMRES();
@@ -53,19 +52,13 @@ class GMRES : public Object, 


   void setRestarting( IndexType rest );



   void setMatrix( MatrixPointer& matrix );

   void setMatrix( MatrixPointer matrix );



   void setPreconditioner( const PreconditionerType& preconditioner );



#ifdef HAVE_NOT_CXX11

   template< typename VectorPointer,

             typename ResidueGetter >

   template< typename Vector,

             typename ResidueGetter = LinearResidueGetter< Matrix, Vector >  >

   bool solve( const Vector& b, Vector& x );

#else

   template< typename VectorPointer,

             typename ResidueGetter = LinearResidueGetter< Matrix, VectorPointer >  >

   bool solve( const VectorPointer& b, VectorPointer& x );

#endif



   ~GMRES();