Commit 0a561076 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Refactoring GMRES variants

parent e7c80efe

src/TNL/Solvers/Linear/CMakeLists.txt

+0 −2
Original line number Diff line number Diff line
@@ -6,8 +6,6 @@ SET( headers BICGStab.h 
             BICGStabL_impl.h

             CG.h

             CG_impl.h

             CWYGMRES.h

             CWYGMRES_impl.h

             GMRES.h

             GMRES_impl.h

             Jacobi.h

src/TNL/Solvers/Linear/CWYGMRES.h

deleted100644 → 0
+0 −115
Original line number Diff line number Diff line 
/***************************************************************************

                          CWYGMRES.h  -  description

                             -------------------

    begin                : May 13, 2016

    copyright            : (C) 2016 by Tomas Oberhuber et al.

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



// Implemented by: Jakub Klinkovsky



#pragma once



#include "LinearSolver.h"



#include <TNL/Containers/Vector.h>



namespace TNL {

namespace Solvers {

namespace Linear {



template< typename Matrix >

class CWYGMRES

: public LinearSolver< Matrix >

{

   using Base = LinearSolver< Matrix >;

public:

   using RealType = typename Base::RealType;

   using DeviceType = typename Base::DeviceType;

   using IndexType = typename Base::IndexType;

   using VectorViewType = typename Base::VectorViewType;

   using ConstVectorViewType = typename Base::ConstVectorViewType;



   String getType() const;



   static void configSetup( Config::ConfigDescription& config,

                            const String& prefix = "" );



   bool setup( const Config::ParameterContainer& parameters,

               const String& prefix = "" ) override;



   void setRestarting( IndexType rest );



   bool solve( ConstVectorViewType b, VectorViewType x ) override;



protected:

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

   using HostVector = Containers::Vector< RealType, Devices::Host, IndexType >;



   void hauseholder_generate( DeviceVector& Y,

                              HostVector& T,

                              const int& i,

                              DeviceVector& w );



   void hauseholder_apply_trunc( HostVector& out,

                                 DeviceVector& Y,

                                 HostVector& T,

                                 const int& i,

                                 DeviceVector& w );



   void hauseholder_cwy( DeviceVector& w,

                         DeviceVector& Y,

                         HostVector& T,

                         const int& i );



   void hauseholder_cwy_transposed( DeviceVector& w,

                                    DeviceVector& Y,

                                    HostVector& T,

                                    const int& i,

                                    DeviceVector& z );



   template< typename Vector >

   void update( IndexType k,

                IndexType m,

                const HostVector& H,

                const HostVector& s,

                DeviceVector& v,

                Vector& x );



   void generatePlaneRotation( RealType& dx,

                               RealType& dy,

                               RealType& cs,

                               RealType& sn );



   void applyPlaneRotation( RealType& dx,

                            RealType& dy,

                            RealType& cs,

                            RealType& sn );





   void setSize( IndexType _size, IndexType m );



   // single vectors

   DeviceVector r, z, w, _M_tmp;

   // matrices (in column-major format)

   DeviceVector V, Y;

   // duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access

   HostVector YL, T;

   // host-only storage for Givens rotations and the least squares problem

   HostVector cs, sn, H, s;



   IndexType size = 0;

   IndexType ldSize = 0;

   IndexType restarting_min = 10;

   IndexType restarting_max = 10;

   IndexType restarting_step_min = 3;

   IndexType restarting_step_max = 3;

};



} // namespace Linear

} // namespace Solvers

} // namespace TNL



#include "CWYGMRES_impl.h"

src/TNL/Solvers/Linear/CWYGMRES_impl.h

deleted100644 → 0
+0 −633

File deleted.

Preview size limit exceeded, changes collapsed.

src/TNL/Solvers/Linear/GMRES.h

+68 −18
Original line number Diff line number Diff line 
/***************************************************************************

                          GMRES.h  -  description

                             -------------------

    begin                : 2007/07/31

    copyright            : (C) 2007 by Tomas Oberhuber

    begin                : May 13, 2016

    copyright            : (C) 2016 by Tomas Oberhuber et al.

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



// Implemented by: Jakub Klinkovsky



#pragma once



#include "LinearSolver.h"
@@ -38,18 +40,54 @@ public: 
   bool setup( const Config::ParameterContainer& parameters,

               const String& prefix = "" ) override;



   void setRestarting( IndexType rest );



   bool solve( ConstVectorViewType b, VectorViewType x ) override;



protected:

   template< typename VectorT >

   void update( IndexType k,

                IndexType m,

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

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

                Containers::Vector< RealType, DeviceType, IndexType >& v,

                VectorT& x );

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

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

   using HostView = Containers::VectorView< RealType, Devices::Host, IndexType >;

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

   using HostVector = Containers::Vector< RealType, Devices::Host, IndexType >;



   enum class Variant { MGS, MGSR, CWY };



   int orthogonalize_MGS( const int m, const RealType normb, const RealType beta );



   int orthogonalize_CWY( const int m, const RealType normb, const RealType beta );



   void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b );



   void preconditioned_matvec( DeviceVector& w, ConstDeviceView v );



   void hauseholder_generate( DeviceVector& Y,

                              HostVector& T,

                              const int i,

                              DeviceVector& w );



   void hauseholder_apply_trunc( HostView out,

                                 DeviceVector& Y,

                                 HostVector& T,

                                 const int i,

                                 DeviceVector& w );



   void hauseholder_cwy( DeviceVector& w,

                         DeviceVector& Y,

                         HostVector& T,

                         const int i );



   void hauseholder_cwy_transposed( DeviceVector& w,

                                    DeviceVector& Y,

                                    HostVector& T,

                                    const int i,

                                    DeviceVector& z );



   template< typename Vector >

   void update( const int k,

                const int m,

                const HostVector& H,

                const HostVector& s,

                DeviceVector& v,

                Vector& x );



   void generatePlaneRotation( RealType& dx,

                               RealType& dy,
@@ -61,17 +99,29 @@ protected: 
                            RealType& cs,

                            RealType& sn );



   void apply_givens_rotations( const int i, const int m );





   void setSize( const IndexType size );



   void setSize( IndexType _size, IndexType m );

   // selected GMRES variant

   Variant variant = Variant::CWY;



   Containers::Vector< RealType, DeviceType, IndexType > _r, w, _v, _M_tmp;

   Containers::Vector< RealType, Devices::Host, IndexType > _s, _cs, _sn, _H;

   // single vectors

   DeviceVector r, w, z, _M_tmp;

   // matrices (in column-major format)

   DeviceVector V, Y;

   // (CWY only) duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access

   HostVector YL, T;

   // host-only storage for Givens rotations and the least squares problem

   HostVector cs, sn, H, s;



   IndexType size = 0;

   IndexType restarting_min = 10;

   IndexType restarting_max = 10;

   IndexType restarting_step_min = 3;

   IndexType restarting_step_max = 3;

   IndexType ldSize = 0;

   int restarting_min = 10;

   int restarting_max = 10;

   int restarting_step_min = 3;

   int restarting_step_max = 3;

};



} // namespace Linear

src/TNL/Solvers/Linear/GMRES_impl.h

+525 −217

File changed.

Preview size limit exceeded, changes collapsed.