Commit ed489170 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Style changes in GMRES - no modifications of the algorithms

parent 4c2bccf2

src/TNL/Containers/Algorithms/cuda-reduction_impl.h

+0 −3
Original line number Diff line number Diff line
@@ -39,8 +39,6 @@ namespace Algorithms { 
 */

const int minGPUReductionDataSize = 256;//65536; //16384;//1024;//256;



//static CudaReductionBuffer cudaReductionBuffer( 8 * minGPUReductionDataSize );



#ifdef HAVE_CUDA



template< typename Operation, int blockSize >
@@ -72,7 +70,6 @@ typename Operation::IndexType reduceOnCudaDevice( Operation& operation, 
   // create reference to the reduction buffer singleton and set default size

   CudaReductionBuffer & cudaReductionBuffer = CudaReductionBuffer::getInstance( 8 * minGPUReductionDataSize );

 

   //CudaReductionBuffer cudaReductionBuffer( 8 * minGPUReductionDataSize );

   if( ! cudaReductionBuffer.setSize( gridSize.x * sizeof( ResultType ) ) )

      return false;

   output = cudaReductionBuffer.template getData< ResultType >();

src/TNL/Solvers/Linear/GMRES.h

+20 −22
Original line number Diff line number Diff line
@@ -14,7 +14,6 @@ 
#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>
@@ -27,12 +26,12 @@ template< typename Matrix, 
          typename Preconditioner = Preconditioners::Dummy< typename Matrix :: RealType,

                                                            typename Matrix :: DeviceType,

                                                            typename Matrix :: IndexType> >

class GMRES : public Object,

class GMRES

: public Object,

  public IterativeSolver< typename Matrix :: RealType,

                          typename Matrix :: IndexType >

{

public:



   typedef typename Matrix::RealType RealType;

   typedef typename Matrix::IndexType IndexType;

   typedef typename Matrix::DeviceType DeviceType;
@@ -64,7 +63,6 @@ class GMRES : public Object, 
   ~GMRES();



protected:



   template< typename VectorT >

   void update( IndexType k,

                IndexType m,
@@ -109,10 +107,10 @@ namespace TNL { 
namespace Solvers {

namespace Linear {

   

extern template class GMRES< Matrices::CSR< float,  Devices::Host, int > >;

/*extern template class GMRES< Matrices::CSR< float,  Devices::Host, int > >;

extern template class GMRES< Matrices::CSR< double, Devices::Host, int > >;

extern template class GMRES< Matrices::CSR< float,  Devices::Host, long int > >;

extern template class GMRES< Matrices::CSR< double, Devices::Host, long int > >;

extern template class GMRES< Matrices::CSR< double, Devices::Host, long int > >;*/



/*extern template class GMRES< Ellpack< float,  Devices::Host, int > >;

extern template class GMRES< Ellpack< double, Devices::Host, int > >;

src/TNL/Solvers/Linear/GMRES_impl.h

+92 −78
Original line number Diff line number Diff line
@@ -16,11 +16,19 @@ namespace Linear { 


template< typename Matrix,

           typename Preconditioner >

GMRES< Matrix, Preconditioner > :: GMRES()

GMRES< Matrix, Preconditioner >::

GMRES()

: size( 0 ),

  restarting( 10 )

{

};

}



template< typename Matrix,

          typename Preconditioner >

GMRES< Matrix, Preconditioner >::

~GMRES()

{

}



template< typename Matrix,

          typename Preconditioner >
@@ -58,23 +66,29 @@ setup( const Config::ParameterContainer& parameters, 


template< typename Matrix,

          typename Preconditioner >

void GMRES< Matrix, Preconditioner > :: setRestarting( IndexType rest )

void

GMRES< Matrix, Preconditioner >::

setRestarting( IndexType rest )

{

   if( size != 0 )

      setSize( size, rest );

   restarting = rest;

};

}



template< typename Matrix,

          typename Preconditioner >

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

void

GMRES< Matrix, Preconditioner >::

setMatrix( const MatrixPointer& matrix )

{

   this->matrix = matrix;

}



template< typename Matrix,

          typename Preconditioner >

void GMRES< Matrix, Preconditioner > :: setPreconditioner( const PreconditionerPointer& preconditioner )

void

GMRES< Matrix, Preconditioner >::

setPreconditioner( const PreconditionerPointer& preconditioner )

{

   this->preconditioner = preconditioner;

}
@@ -82,8 +96,11 @@ void GMRES< Matrix, Preconditioner > :: setPreconditioner( const PreconditionerP 
template< typename Matrix,

          typename Preconditioner >

   template< typename Vector, typename ResidueGetter >

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

bool

GMRES< Matrix, Preconditioner >::

solve( const Vector& b, Vector& x )

{

   Assert( matrix, std::cerr << "No matrix was set in GMRES. Call setMatrix() before solve()." << std::endl );

   if( restarting <= 0 )

   {

      std::cerr << "I have wrong value for the restarting of the GMRES solver. It is set to " << restarting
@@ -107,7 +124,6 @@ bool GMRES< Matrix, Preconditioner >::solve( const Vector& b, Vector& x ) 
   RealType *M_tmp = _M_tmp.getData();



   RealType normb( 0.0 ), beta( 0.0 );

   //T normb( 0.0 ), beta( 0.0 ); does not work with openmp yet

   /****

    * 1. Solve r from M r = b - A x_0

    */
@@ -118,20 +134,16 @@ bool GMRES< Matrix, Preconditioner >::solve( const Vector& b, Vector& x ) 


      matrix -> vectorProduct( x, _M_tmp );

      _M_tmp.addVector( b, ( RealType ) 1.0, -1.0 );

      /*for( i = 0; i < size; i ++ )

         M_tmp[ i ] = b[ i ] - M_tmp[ i ];*/



      this->preconditioner->solve( _M_tmp, _r );

      beta = _r. lpNorm( ( RealType ) 2.0 );

   }

   else

   {

      matrix -> vectorProduct( x, _r );

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

      _r.addVector( b, ( RealType ) 1.0, -1.0 );

      beta = _r. lpNorm( ( RealType ) 2.0 );

      //cout << "x = " << x << std::endl;

   }

   beta = _r.lpNorm( ( RealType ) 2.0 );

 

   //cout << "norm b = " << normb << std::endl;

   //cout << " beta = " << beta << std::endl;
@@ -142,7 +154,7 @@ bool GMRES< Matrix, Preconditioner >::solve( const Vector& b, Vector& x ) 
   this->resetIterations();

   this->setResidue( beta / normb );



   Containers::SharedVector< RealType, DeviceType, IndexType > vi, vk;

   Containers::Vector< RealType, DeviceType, IndexType > vi, vk;

   while( this->checkNextIteration() )

   {

      const IndexType m = restarting;
@@ -283,18 +295,14 @@ bool GMRES< Matrix, Preconditioner >::solve( const Vector& b, Vector& x ) 
   }

   this->refreshSolverMonitor( true );

   return this->checkConvergence();

};



template< typename Matrix,

          typename Preconditioner >

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

{

};

}



template< typename Matrix,

          typename Preconditioner >

   template< typename VectorT >

void GMRES< Matrix, Preconditioner > :: update( IndexType k,

void

GMRES< Matrix, Preconditioner >::

update( IndexType k,

        IndexType m,

        const Containers::Vector< RealType, Devices::Host, IndexType >& H,

        const Containers::Vector< RealType, Devices::Host, IndexType >& s,
@@ -317,17 +325,19 @@ void GMRES< Matrix, Preconditioner > :: update( IndexType k, 
         y[ j ] -= H[ j + i * ( m + 1 ) ] * y[ i ];

   }



   Containers::SharedVector< RealType, DeviceType, IndexType > vi;

   Containers::Vector< RealType, DeviceType, IndexType > vi;

   for( i = 0; i <= k; i++)

   {

      vi.bind( &( v.getData()[ i * this->size ] ), x.getSize() );

      x.addVector( vi, y[ i ] );

   }

};

}



template< typename Matrix,

          typename Preconditioner >

void GMRES< Matrix, Preconditioner > :: generatePlaneRotation( RealType &dx,

void

GMRES< Matrix, Preconditioner >::

generatePlaneRotation( RealType& dx,

                       RealType& dy,

                       RealType& cs,

                       RealType& sn )
@@ -350,11 +360,13 @@ void GMRES< Matrix, Preconditioner > :: generatePlaneRotation( RealType &dx, 
         cs = 1.0 / std::sqrt( 1.0 + temp * temp );

         sn = temp * cs;

      }

};

}



template< typename Matrix,

          typename Preconditioner >

void GMRES< Matrix, Preconditioner > :: applyPlaneRotation( RealType &dx,

void

GMRES< Matrix, Preconditioner >::

applyPlaneRotation( RealType& dx,

                    RealType& dy,

                    RealType& cs,

                    RealType& sn )
@@ -362,11 +374,13 @@ void GMRES< Matrix, Preconditioner > :: applyPlaneRotation( RealType &dx, 
   RealType temp  =  cs * dx + sn * dy;

   dy =  cs * dy - sn * dx;

   dx = temp;

};

}



template< typename Matrix,

          typename Preconditioner >

bool GMRES< Matrix, Preconditioner > :: setSize( IndexType _size, IndexType m )

bool

GMRES< Matrix, Preconditioner >::

setSize( IndexType _size, IndexType m )

{

   if( size == _size && restarting == m ) return true;

   size = _size;
@@ -384,7 +398,7 @@ bool GMRES< Matrix, Preconditioner > :: setSize( IndexType _size, IndexType m ) 
      return false;

   }

   return true;

};

}



} // namespace Linear

} // namespace Solvers