Skip to content

Commit

Permalink
Finished converting GMRES-variants to dense matrix traits
Browse files Browse the repository at this point in the history 
The pseudo-block and polynomial-preconditioned GMRES variants
are integrated with dense matrix traits.  All GMRES variants have
been cleaned up to remove default templates where appropriate.
  • Loading branch information
@hkthorn
hkthorn committed Feb 28, 2024
1 parent 0e89333 commit bc36ff2
Show file tree
Hide file tree
Showing 12 changed files with 788 additions and 162 deletions.
3 changes: 1 addition & 2 deletions packages/belos/src/BelosBlockFGmresIter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@
#include "BelosOperatorTraits.hpp"
#include "BelosMultiVecTraits.hpp"
#include "BelosDenseMatTraits.hpp"
#include "BelosTeuchosDenseAdapter.hpp"

#include "Teuchos_BLAS.hpp"
#include "Teuchos_ScalarTraits.hpp"
Expand All @@ -81,7 +80,7 @@

namespace Belos {

template<class ScalarType, class MV, class OP, class DM = Teuchos::SerialDenseMatrix<int, ScalarType>>
template<class ScalarType, class MV, class OP, class DM>
class BlockFGmresIter : virtual public GmresIteration<ScalarType,MV,OP,DM> {

public:
Expand Down
9 changes: 4 additions & 5 deletions packages/belos/src/BelosBlockGCRODRSolMgr.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,6 @@
#include "BelosGmresIteration.hpp"
#include "BelosBlockGCRODRIter.hpp"
#include "BelosBlockGmresIter.hpp"
#include "BelosBlockFGmresIter.hpp"
#include "BelosStatusTestMaxIters.hpp"
#include "BelosStatusTestGenResNorm.hpp"
#include "BelosStatusTestCombo.hpp"
Expand Down Expand Up @@ -296,7 +295,7 @@ class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP, DM> {
// "AugKryl" indicates it is specialized for building a recycle space from the augmented Krylov subspace

// Functions which control the building of a recycle space
void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter);
void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP,DM> > block_gmres_iter);
void buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > gcrodr_iter);

// Recycling with Harmonic Ritz Vectors
Expand Down Expand Up @@ -1205,7 +1204,7 @@ class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP, DM> {
}

template<class ScalarType, class MV, class OP>
void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP,DM> > block_gmres_iter){

ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
Expand Down Expand Up @@ -1972,8 +1971,8 @@ ReturnType BlockGCRODRSolMgr<ScalarType,MV,OP>::solve() {
primeList.set("Num Blocks",numBlocks_-1);
}
//Create Block GMRES iteration object to perform one cycle of GMRES
Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP,DM> > block_gmres_iter;
block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP,DM>(problem_,printer_,outputTest_,ortho_,primeList) );

// MLP: ADD LOGIC TO DEAL WITH USER ASKING TO GENERATE A LARGER SPACE THAN dim AS IN HEIDI'S BlockGmresSolMgr CODE (DIDN'T WE ALREADY DO THIS SOMEWHERE?)
block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
Expand Down
3 changes: 1 addition & 2 deletions packages/belos/src/BelosBlockGmresIter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@
#include "BelosOperatorTraits.hpp"
#include "BelosMultiVecTraits.hpp"
#include "BelosDenseMatTraits.hpp"
#include "BelosTeuchosDenseAdapter.hpp"

#include "Teuchos_BLAS.hpp"
#include "Teuchos_LAPACK.hpp"
Expand All @@ -82,7 +81,7 @@

namespace Belos {

template<class ScalarType, class MV, class OP, class DM = Teuchos::SerialDenseMatrix<int, ScalarType>>
template<class ScalarType, class MV, class OP, class DM>
class BlockGmresIter : virtual public GmresIteration<ScalarType,MV,OP,DM> {

public:
Expand Down
13 changes: 6 additions & 7 deletions packages/belos/src/BelosGmresIteration.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@
#include "BelosTypes.hpp"
#include "BelosIteration.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_SerialDenseVector.hpp"

namespace Belos {

Expand Down Expand Up @@ -103,7 +102,7 @@ namespace Belos {
*
* This struct is utilized by PseudoBlockGmresIter::initialize() and PseudoBlockGmresIter::getState().
*/
template <class ScalarType, class MV>
template <class ScalarType, class MV, class DM = Teuchos::SerialDenseMatrix<int, ScalarType>>
struct PseudoBlockGmresIterState {

typedef Teuchos::ScalarTraits<ScalarType> SCT;
Expand All @@ -121,14 +120,14 @@ namespace Belos {
* The \c curDim by \c curDim leading submatrix of H is the
* projection of problem->getOperator() by the first \c curDim vectors in V.
*/
std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > H;
std::vector<Teuchos::RCP<const DM> > H;
/*! \brief The current upper-triangular matrix from the QR reduction of H. */
std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > R;
std::vector<Teuchos::RCP<const DM> > R;
/*! \brief The current right-hand side of the least squares system RY = Z. */
std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > Z;
std::vector<Teuchos::RCP<const DM> > Z;
/*! \brief The current Given's rotation coefficients. */
std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > sn;
std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > > cs;
std::vector<Teuchos::RCP<const std::vector<ScalarType> > > sn;
std::vector<Teuchos::RCP<const std::vector<MagnitudeType> > > cs;

PseudoBlockGmresIterState() : curDim(0), V(0),
H(0), R(0), Z(0),
Expand Down
31 changes: 17 additions & 14 deletions packages/belos/src/BelosGmresPolyOp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,12 @@ namespace Belos {

// Create a shell class for the MV, inherited off MultiVec<> that will operate with the GmresPolyOp.
template <class ScalarType, class MV>
class GmresPolyMv : public MultiVec< ScalarType >
class GmresPolyMv : public MultiVec< ScalarType, Teuchos::SerialDenseMatrix<int,ScalarType> >
{
public:

typedef Teuchos::SerialDenseMatrix<int,ScalarType> DM;

GmresPolyMv ( const Teuchos::RCP<MV>& mv_in )
: mv_(mv_in)
{}
Expand Down Expand Up @@ -136,7 +138,7 @@ namespace Belos {
int GetNumberVecs () const { return MVT::GetNumberVecs( *mv_ ); }
void MvTimesMatAddMv (const ScalarType alpha,
const MultiVec<ScalarType>& A,
const Teuchos::SerialDenseMatrix<int,ScalarType>& B, const ScalarType beta)
const DM& B, const ScalarType beta)
{
const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
MVT::MvTimesMatAddMv( alpha, *(A_in.getConstMV()), B, beta, *mv_ );
Expand All @@ -149,7 +151,7 @@ namespace Belos {
}
void MvScale ( const ScalarType alpha ) { MVT::MvScale( *mv_, alpha ); }
void MvScale ( const std::vector<ScalarType>& alpha ) { MVT::MvScale( *mv_, alpha ); }
void MvTransMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, Teuchos::SerialDenseMatrix<int,ScalarType>& B) const
void MvTransMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, DM& B) const
{
const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
MVT::MvTransMv( alpha, *(A_in.getConstMV()), *mv_, B );
Expand All @@ -174,7 +176,7 @@ namespace Belos {

private:

typedef MultiVecTraits<ScalarType,MV> MVT;
typedef MultiVecTraits<ScalarType,MV,DM> MVT;

Teuchos::RCP<MV> mv_;

Expand Down Expand Up @@ -293,6 +295,7 @@ namespace Belos {
typedef Teuchos::ScalarTraits<ScalarType> SCT ;
typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
typedef Teuchos::ScalarTraits<MagnitudeType> MCT ;
typedef Teuchos::SerialDenseMatrix<int,ScalarType> DM;

// Default polynomial parameters
static constexpr int maxDegree_default_ = 25;
Expand All @@ -314,7 +317,7 @@ namespace Belos {
Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcpFromRef(std::cout);

// Orthogonalization manager.
Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > ortho_;

// Current polynomial parameters
MagnitudeType polyTol_ = DefaultSolverParameters::polyTol;
Expand Down Expand Up @@ -551,28 +554,28 @@ namespace Belos {
// Create orthogonalization manager if we need to.
if (ortho_.is_null()) {
params_->set("Orthogonalization", orthoType_);
Belos::OrthoManagerFactory<ScalarType, MV, OP> factory;
Belos::OrthoManagerFactory<ScalarType, MV, OP, DM> factory;
Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null

ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
}

// Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP,DM> > maxItrTst =
Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP,DM>( maxDegree_ ) );

// Implicit residual test, using the native residual to determine if convergence was achieved.
Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polyTol_ ) );
Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP,DM> > convTst =
Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP,DM>( polyTol_ ) );
convTst->defineScaleForm( convertStringToScaleType("Norm of RHS"), Belos::TwoNorm );

// Convergence test that stops the iteration when either are satisfied.
Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP,DM> > polyTest =
Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP,DM>( StatusTestCombo<ScalarType,MV,OP,DM>::OR, maxItrTst, convTst ) );

// Create Gmres iteration object to perform one cycle of Gmres.
Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP,DM> > gmres_iter;
gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP,DM>(newProblem,printer_,polyTest,ortho_,polyList) );

// Create the first block in the current Krylov basis (residual).
Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
Expand Down
13 changes: 7 additions & 6 deletions packages/belos/src/BelosGmresPolySolMgr.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,6 +155,7 @@ class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP,DM> {
typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;

typedef Teuchos::ScalarTraits<MagnitudeType> MTS;
typedef Teuchos::SerialDenseMatrix<int,ScalarType> dm_t; // Internally, we use a different SDM
typedef Belos::GmresPolyOp<ScalarType,MV,OP> gmres_poly_t;
typedef Belos::GmresPolyMv<ScalarType,MV> gmres_poly_mv_t;

Expand Down Expand Up @@ -674,8 +675,8 @@ ReturnType GmresPolySolMgr<ScalarType,MV,OP,DM>::solve ()

// Then the polynomial will be used as an operator for an outer solver.
// Use outer solver parameter list passed in a sublist.
Belos::GenericSolverFactory<ScalarType, MultiVec<ScalarType>, Operator<ScalarType> > factory;
RCP<SolverManager<ScalarType, MultiVec<ScalarType>, Operator<ScalarType> > > solver = factory.create( outerSolverType_, outerParams_ );
Belos::GenericSolverFactory<ScalarType, MultiVec<ScalarType,dm_t>, Operator<ScalarType> > factory;
RCP<SolverManager<ScalarType, MultiVec<ScalarType,dm_t>, Operator<ScalarType> > > solver = factory.create( outerSolverType_, outerParams_ );
TEUCHOS_TEST_FOR_EXCEPTION( solver == Teuchos::null, std::invalid_argument,
"Belos::GmresPolySolMgr::solve(): Selected solver is not valid.");

Expand All @@ -684,8 +685,8 @@ ReturnType GmresPolySolMgr<ScalarType,MV,OP,DM>::solve ()
RCP<gmres_poly_mv_t> new_lhs = rcp( new gmres_poly_mv_t( problem_->getLHS() ) );
RCP<gmres_poly_mv_t> new_rhs = rcp( new gmres_poly_mv_t( rcp_const_cast<MV>( problem_->getRHS() ) ) );
RCP<gmres_poly_t> A = rcp( new gmres_poly_t( problem_ ) ); // This just performs problem_->applyOp
RCP<LinearProblem<ScalarType,MultiVec<ScalarType>,Operator<ScalarType> > > newProblem =
rcp( new LinearProblem<ScalarType,MultiVec<ScalarType>,Operator<ScalarType> >( A, new_lhs, new_rhs ) );
RCP<LinearProblem<ScalarType,MultiVec<ScalarType,dm_t>,Operator<ScalarType> > > newProblem =
rcp( new LinearProblem<ScalarType,MultiVec<ScalarType,dm_t>,Operator<ScalarType> >( A, new_lhs, new_rhs ) );
std::string solverLabel = label_ + ": Hybrid Gmres";
newProblem->setLabel(solverLabel);

Expand All @@ -711,8 +712,8 @@ ReturnType GmresPolySolMgr<ScalarType,MV,OP,DM>::solve ()
else if (hasOuterSolver_) {

// There is no polynomial, just create the outer solver with the outerSolverType_ and outerParams_.
Belos::SolverFactory<ScalarType, MV, OP, DM> factory;
RCP<SolverManager<ScalarType, MV, OP, DM> > solver = factory.create( outerSolverType_, outerParams_ );
Belos::SolverFactory<ScalarType, MV, OP, dm_t> factory;
RCP<SolverManager<ScalarType, MV, OP, dm_t> > solver = factory.create( outerSolverType_, outerParams_ );
TEUCHOS_TEST_FOR_EXCEPTION( solver == Teuchos::null, std::invalid_argument,
"Belos::GmresPolySolMgr::solve(): Selected solver is not valid.");

Expand Down
Loading

0 comments on commit bc36ff2

Please sign in to comment.