Commit 8ae818f1 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Updating linear solvers and preconditioners for distributed matrices/vectors

parent 63dda575

src/Benchmarks/LinearSolvers/tnl-benchmark-linear-solvers.h

+46 −46
Original line number Diff line number Diff line
@@ -99,17 +99,17 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 
   benchmarkPreconditionerUpdate< Diagonal >( benchmark, parameters, cudaMatrixPointer );

#endif



   benchmark.setOperation("GMRES (Jacobi)");

   benchmarkSolver< GMRES, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< GMRES, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif



   benchmark.setOperation("CWYGMRES (Jacobi)");

   benchmarkSolver< CWYGMRES, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< CWYGMRES, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//   benchmark.setOperation("GMRES (Jacobi)");

//   benchmarkSolver< GMRES, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< GMRES, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



//   benchmark.setOperation("CWYGMRES (Jacobi)");

//   benchmarkSolver< CWYGMRES, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< CWYGMRES, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



   benchmark.setOperation("TFQMR (Jacobi)");

   benchmarkSolver< TFQMR, Diagonal >( benchmark, parameters, matrixPointer, x0, b );
@@ -125,11 +125,11 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 


   for( int ell = 1; ell <= ell_max; ell++ ) {

      parameters.template setParameter< int >( "bicgstab-ell", ell );

      benchmark.setOperation("BiCGstab(" + String(ell) + ") (Jacobi)");

      benchmarkSolver< BICGStabL, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

      benchmarkSolver< BICGStabL, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//      benchmark.setOperation("BiCGstab(" + String(ell) + ") (Jacobi)");

//      benchmarkSolver< BICGStabL, Diagonal >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//      benchmarkSolver< BICGStabL, Diagonal >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif

   }
@@ -139,17 +139,17 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 
   benchmarkPreconditionerUpdate< ILU0 >( benchmark, parameters, cudaMatrixPointer );

#endif



   benchmark.setOperation("GMRES (ILU0)");

   benchmarkSolver< GMRES, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< GMRES, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//   benchmark.setOperation("GMRES (ILU0)");

//   benchmarkSolver< GMRES, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< GMRES, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



   benchmark.setOperation("CWYGMRES (ILU0)");

   benchmarkSolver< CWYGMRES, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< CWYGMRES, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//   benchmark.setOperation("CWYGMRES (ILU0)");

//   benchmarkSolver< CWYGMRES, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< CWYGMRES, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



   benchmark.setOperation("TFQMR (ILU0)");

   benchmarkSolver< TFQMR, ILU0 >( benchmark, parameters, matrixPointer, x0, b );
@@ -165,11 +165,11 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 


   for( int ell = 1; ell <= ell_max; ell++ ) {

      parameters.template setParameter< int >( "bicgstab-ell", ell );

      benchmark.setOperation("BiCGstab(" + String(ell) + ") (ILU0)");

      benchmarkSolver< BICGStabL, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

      benchmarkSolver< BICGStabL, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//      benchmark.setOperation("BiCGstab(" + String(ell) + ") (ILU0)");

//      benchmarkSolver< BICGStabL, ILU0 >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//      benchmarkSolver< BICGStabL, ILU0 >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif

   }
@@ -179,17 +179,17 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 
   benchmarkPreconditionerUpdate< ILUT >( benchmark, parameters, cudaMatrixPointer );

#endif



   benchmark.setOperation("GMRES (ILUT)");

   benchmarkSolver< GMRES, ILUT >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< GMRES, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//   benchmark.setOperation("GMRES (ILUT)");

//   benchmarkSolver< GMRES, ILUT >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< GMRES, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



   benchmark.setOperation("CWYGMRES (ILUT)");

   benchmarkSolver< CWYGMRES, ILUT >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

   benchmarkSolver< CWYGMRES, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//   benchmark.setOperation("CWYGMRES (ILUT)");

//   benchmarkSolver< CWYGMRES, ILUT >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//   benchmarkSolver< CWYGMRES, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif



   benchmark.setOperation("TFQMR (ILUT)");

   benchmarkSolver< TFQMR, ILUT >( benchmark, parameters, matrixPointer, x0, b );
@@ -205,11 +205,11 @@ benchmarkIterativeSolvers( Benchmark& benchmark, 


   for( int ell = 1; ell <= ell_max; ell++ ) {

      parameters.template setParameter< int >( "bicgstab-ell", ell );

      benchmark.setOperation("BiCGstab(" + String(ell) + ") (ILUT)");

      benchmarkSolver< BICGStabL, ILUT >( benchmark, parameters, matrixPointer, x0, b );

#ifdef HAVE_CUDA

      benchmarkSolver< BICGStabL, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

#endif

//      benchmark.setOperation("BiCGstab(" + String(ell) + ") (ILUT)");

//      benchmarkSolver< BICGStabL, ILUT >( benchmark, parameters, matrixPointer, x0, b );

//#ifdef HAVE_CUDA

//      benchmarkSolver< BICGStabL, ILUT >( benchmark, parameters, cudaMatrixPointer, cuda_x0, cuda_b );

//#endif

   }

}

src/TNL/Solvers/Linear/BICGStab.h

+2 −4
Original line number Diff line number Diff line
@@ -12,8 +12,6 @@ 


#include "LinearSolver.h"



#include <TNL/Containers/Vector.h>



namespace TNL {

namespace Solvers {

namespace Linear {
@@ -41,11 +39,11 @@ public: 
   bool solve( ConstVectorViewType b, VectorViewType x ) override;



protected:

   void setSize( IndexType size );

   void setSize( const VectorViewType& x );



   bool exact_residue = false;



   Containers::Vector< RealType, DeviceType, IndexType > r, r_ast, p, s, Ap, As, M_tmp;

   typename Traits< Matrix >::VectorType r, r_ast, p, s, Ap, As, M_tmp;

};



} // namespace Linear

src/TNL/Solvers/Linear/BICGStab_impl.h

+9 −9
Original line number Diff line number Diff line
@@ -48,7 +48,7 @@ setup( const Config::ParameterContainer& parameters, 
template< typename Matrix >

bool BICGStab< Matrix >::solve( ConstVectorViewType b, VectorViewType x )

{

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

   this->setSize( x );



   RealType alpha, beta, omega, aux, rho, rho_old, b_norm;
@@ -161,15 +161,15 @@ bool BICGStab< Matrix >::solve( ConstVectorViewType b, VectorViewType x ) 
}



template< typename Matrix >

void BICGStab< Matrix > :: setSize( IndexType size )

void BICGStab< Matrix > :: setSize( const VectorViewType& x )

{

   r.setSize( size );

   r_ast.setSize( size );

   p.setSize( size );

   s.setSize( size );

   Ap.setSize( size );

   As.setSize( size );

   M_tmp.setSize( size );

   r.setLike( x );

   r_ast.setLike( x );

   p.setLike( x );

   s.setLike( x );

   Ap.setLike( x );

   As.setLike( x );

   M_tmp.setLike( x );

}



} // namespace Linear

src/TNL/Solvers/Linear/LinearSolver.h

+4 −3
Original line number Diff line number Diff line
@@ -17,9 +17,10 @@ 


#include <TNL/Solvers/IterativeSolver.h>

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

#include <TNL/Containers/VectorView.h>

#include <TNL/Pointers/SharedPointer.h>



#include "Traits.h"



namespace TNL {

namespace Solvers {

namespace Linear {
@@ -32,8 +33,8 @@ public: 
   using RealType = typename Matrix::RealType;

   using DeviceType = typename Matrix::DeviceType;

   using IndexType = typename Matrix::IndexType;

   using VectorViewType = Containers::VectorView< RealType, DeviceType, IndexType >;

   using ConstVectorViewType = Containers::VectorView< typename std::add_const< RealType >::type, DeviceType, IndexType >;

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

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

   using MatrixType = Matrix;

   using MatrixPointer = Pointers::SharedPointer< typename std::add_const< MatrixType >::type >;

   using PreconditionerType = Preconditioners::Preconditioner< MatrixType >;

src/TNL/Solvers/Linear/Preconditioners/Diagonal.h

+29 −0
Original line number Diff line number Diff line
@@ -47,6 +47,35 @@ protected: 
   VectorType diagonal;

};



template< typename Matrix, typename Communicator >

class Diagonal< DistributedContainers::DistributedMatrix< Matrix, Communicator > >

: public Preconditioner< DistributedContainers::DistributedMatrix< Matrix, Communicator > >

{

public:

   using MatrixType = DistributedContainers::DistributedMatrix< Matrix, Communicator >;

   using RealType = typename MatrixType::RealType;

   using DeviceType = typename MatrixType::DeviceType;

   using IndexType = typename MatrixType::IndexType;

   using typename Preconditioner< MatrixType >::VectorViewType;

   using typename Preconditioner< MatrixType >::ConstVectorViewType;

   using typename Preconditioner< MatrixType >::MatrixPointer;

   using VectorType = Containers::Vector< RealType, DeviceType, IndexType >;

   using LocalVectorViewType = Containers::VectorView< RealType, DeviceType, IndexType >;

   using ConstLocalVectorViewType = Containers::VectorView< typename std::add_const< RealType >::type, DeviceType, IndexType >;



   virtual void update( const MatrixPointer& matrixPointer ) override;



   virtual void solve( ConstVectorViewType b, VectorViewType x ) const override;



   String getType() const

   {

      return String( "Diagonal" );

   }



protected:

   VectorType diagonal;

};



} // namespace Preconditioners

} // namespace Linear

} // namespace Solvers