Skip to content

Commit

Permalink
Converts PCPG solver to use dense matrix traits
Browse files Browse the repository at this point in the history
  • Loading branch information
@hkthorn
hkthorn committed Oct 17, 2023
1 parent b1b7e86 commit 801ff16
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 84 deletions.
143 changes: 65 additions & 78 deletions packages/belos/src/BelosPCPGIter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,9 @@
#include "BelosStatusTest.hpp"
#include "BelosOperatorTraits.hpp"
#include "BelosMultiVecTraits.hpp"
#include "BelosDenseMatTraits.hpp"
#include "BelosCGIteration.hpp"

#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"
Expand All @@ -66,7 +66,7 @@
\class Belos::PCPGIter
\brief This class implements the PCPG iteration, where a
single-std::vector Krylov subspace is constructed. The documentation
single-vector Krylov subspace is constructed. The documentation
refers to blocks, but note that at this point, all blocks have unit
dimension.
Expand All @@ -82,7 +82,7 @@ namespace Belos {
*
* The structure is utilized by initialize() and getState().
*/
template <class ScalarType, class MV>
template <class ScalarType, class MV, class DM>
struct PCPGIterState {
/*! \brief The current dimension of the reduction.
*
Expand Down Expand Up @@ -115,7 +115,7 @@ namespace Belos {
*
* The \c curDim by \c curDim D = diag(P'*AP) = U' * C
*/
Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
Teuchos::RCP<const DM> D;

PCPGIterState() : curDim(0),
prevUdim(0),
Expand All @@ -128,16 +128,17 @@ namespace Belos {

//@}

template<class ScalarType, class MV, class OP>
class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
template<class ScalarType, class MV, class OP, class DM>
class PCPGIter : virtual public Iteration<ScalarType,MV,OP,DM> {

public:

//
// Convenience typedefs
//
typedef MultiVecTraits<ScalarType,MV> MVT;
typedef MultiVecTraits<ScalarType,MV,DM> MVT;
typedef OperatorTraits<ScalarType,MV,OP> OPT;
typedef DenseMatTraits<ScalarType,DM> DMT;
typedef Teuchos::ScalarTraits<ScalarType> SCT;
typedef typename SCT::magnitudeType MagnitudeType;

Expand All @@ -153,8 +154,8 @@ namespace Belos {
*/
PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
const Teuchos::RCP<OutputManager<ScalarType> > &printer,
const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
Teuchos::ParameterList &params );

//! Destructor.
Expand Down Expand Up @@ -203,14 +204,14 @@ namespace Belos {
* \note For any pointer in \c newstate which directly points to the multivectors in
* the solver, the data is not (supposed to be) copied.
*/
void initialize(PCPGIterState<ScalarType,MV>& newstate);
void initialize(PCPGIterState<ScalarType,MV,DM>& newstate);

/*! \brief Initialize the solver with the initial vectors from the linear problem.
* An exception is thrown if initialzed is called and newstate.R is null.
*/
void initialize()
{
PCPGIterState<ScalarType,MV> empty;
PCPGIterState<ScalarType,MV,DM> empty;
initialize(empty);
}

Expand All @@ -221,8 +222,8 @@ namespace Belos {
* \returns A PCPGIterState object containing const pointers to the current
* solver state.
*/
PCPGIterState<ScalarType,MV> getState() const {
PCPGIterState<ScalarType,MV> state;
PCPGIterState<ScalarType,MV,DM> getState() const {
PCPGIterState<ScalarType,MV,DM> state;
state.Z = Z_; // CG state
state.P = P_;
state.AP = AP_;
Expand Down Expand Up @@ -316,8 +317,8 @@ namespace Belos {
//
const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
const Teuchos::RCP<OutputManager<ScalarType> > om_;
const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > stest_;
const Teuchos::RCP<OrthoManager<ScalarType,MV,DM> > ortho_;

//
// Algorithmic parameters
Expand Down Expand Up @@ -375,16 +376,16 @@ namespace Belos {
//
// Projected matrices
// D_ : Diagonal matrix of pivots D = P'AP
Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
Teuchos::RCP<DM> D_;
};

//////////////////////////////////////////////////////////////////////////////////////////////////
// Constructor.
template<class ScalarType, class MV, class OP>
PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
template<class ScalarType, class MV, class OP, class DM>
PCPGIter<ScalarType,MV,OP,DM>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
const Teuchos::RCP<OutputManager<ScalarType> > &printer,
const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
Teuchos::ParameterList &params ):
lp_(problem),
om_(printer),
Expand Down Expand Up @@ -417,8 +418,8 @@ namespace Belos {

//////////////////////////////////////////////////////////////////////////////////////////////////
// Set the block size and adjust as necessary
template <class ScalarType, class MV, class OP>
void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
template<class ScalarType, class MV, class OP, class DM>
void PCPGIter<ScalarType,MV,OP,DM>::setSize( int savedBlocks )
{
// allocate space only; perform no computation
// Any change in size invalidates the state of the solver as implemented here.
Expand All @@ -437,8 +438,8 @@ namespace Belos {

//////////////////////////////////////////////////////////////////////////////////////////////////
// Enable the reuse of a single solver object for completely different linear systems
template <class ScalarType, class MV, class OP>
void PCPGIter<ScalarType,MV,OP>::resetState()
template<class ScalarType, class MV, class OP, class DM>
void PCPGIter<ScalarType,MV,OP,DM>::resetState()
{
stateStorageInitialized_ = false;
initialized_ = false;
Expand All @@ -449,8 +450,8 @@ namespace Belos {

//////////////////////////////////////////////////////////////////////////////////////////////////
// Setup the state storage. Called by either initialize or, if savedBlocks_ changes, setSize.
template <class ScalarType, class MV, class OP>
void PCPGIter<ScalarType,MV,OP>::setStateSize ()
template<class ScalarType, class MV, class OP, class DM>
void PCPGIter<ScalarType,MV,OP,DM>::setStateSize ()
{
if (!stateStorageInitialized_) {

Expand Down Expand Up @@ -514,14 +515,14 @@ namespace Belos {
}
if (keepDiagonal_) {
if (D_ == Teuchos::null) {
D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
D_ = DMT::Create();
}
if (initDiagonal_) {
D_->shape( newsd, newsd );
DMT::Reshape( *D_, newsd, newsd, true );
}
else {
if (D_->numRows() < newsd || D_->numCols() < newsd) {
D_->shapeUninitialized( newsd, newsd );
DMT::Reshape( *D_, newsd, newsd, false );
}
}
}
Expand All @@ -533,8 +534,8 @@ namespace Belos {

//////////////////////////////////////////////////////////////////////////////////////////////////
// Initialize the iteration object
template <class ScalarType, class MV, class OP>
void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV>& newstate)
template<class ScalarType, class MV, class OP, class DM>
void PCPGIter<ScalarType,MV,OP,DM>::initialize(PCPGIterState<ScalarType,MV,DM>& newstate)
{

TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
Expand Down Expand Up @@ -621,8 +622,8 @@ namespace Belos {

//////////////////////////////////////////////////////////////////////////////////////////////////
// Iterate until the status test informs us we should stop.
template <class ScalarType, class MV, class OP>
void PCPGIter<ScalarType,MV,OP>::iterate()
template<class ScalarType, class MV, class OP, class DM>
void PCPGIter<ScalarType,MV,OP,DM>::iterate()
{
//
// Allocate/initialize data structures
Expand All @@ -633,10 +634,9 @@ namespace Belos {
const bool debug = false;

// Allocate memory for scalars.
Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
ScalarType alpha, beta, rHz_old;
Teuchos::RCP<DM> pAp = DMT::Create(1,1);
Teuchos::RCP<DM> rHz = DMT::Create(1,1);

if( iter_ != 0 )
std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl; //DMD
Expand All @@ -645,12 +645,12 @@ namespace Belos {
std::vector<int> prevInd;
Teuchos::RCP<const MV> Uprev;
Teuchos::RCP<const MV> Cprev;
Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
Teuchos::RCP<DM> CZ;

if( prevUdim_ ){
prevInd.resize( prevUdim_ );
for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
CZ.reshape( prevUdim_ , 1 );
CZ = DMT::Create( prevUdim_ , 1 );
Uprev = MVT::CloneView(*U_, prevInd);
Cprev = MVT::CloneView(*C_, prevInd);
}
Expand All @@ -677,14 +677,9 @@ namespace Belos {
Teuchos::RCP<MV> P;
curind[0] = curDim_ - 1; // column = dimension - 1
P = MVT::CloneViewNonConst(*U_,curind);
MVT::MvTransMv( one, *Cprev, *P, CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P ); // P -= U*(C'Z)
MVT::MvTransMv( one, *Cprev, *P, *CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *P ); // P -= U*(C'Z)

if( debug ){
MVT::MvTransMv( one, *Cprev, *P, CZ );
std::cout << " Input CZ post ortho " << std::endl;
CZ.print( std::cout );
}
if( curDim_ == savedBlocks_ ){
std::vector<int> zero_index(1);
zero_index[0] = 0;
Expand All @@ -694,7 +689,8 @@ namespace Belos {
}

// Compute first <r,z> a.k.a. rHz
MVT::MvTransMv( one, *R_, *Z_, rHz );
MVT::MvTransMv( one, *R_, *Z_, *rHz );
DMT::SyncDeviceToHost( *rHz );

////////////////////////////////////////////////////////////////
// iterate until the status test is satisfied
Expand All @@ -713,51 +709,52 @@ namespace Belos {
P = MVT::CloneView(*U_,curind);
AP = MVT::CloneViewNonConst(*C_,curind);
lp_->applyOp( *P, *AP );
MVT::MvTransMv( one, *P, *AP, pAp );
MVT::MvTransMv( one, *P, *AP, *pAp );
}else{
if( prevUdim_ + iter_ == savedBlocks_ ){
AP = MVT::CloneViewNonConst(*C_,curind);
lp_->applyOp( *P_, *AP );
MVT::MvTransMv( one, *P_, *AP, pAp );
MVT::MvTransMv( one, *P_, *AP, *pAp );
}else{
lp_->applyOp( *P_, *AP_ );
MVT::MvTransMv( one, *P_, *AP_, pAp );
MVT::MvTransMv( one, *P_, *AP_, *pAp );
}
}
DMT::SyncDeviceToHost( *pAp );

if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
(*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
DMT::Value( *D_, iter_-1 ,iter_-1 ) = DMT::ValueConst(*pAp,0,0);

// positive pAp required
TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, CGPositiveDefiniteFailure,
TEUCHOS_TEST_FOR_EXCEPTION( DMT::ValueConst(*pAp,0,0) <= zero, CGPositiveDefiniteFailure,
"Belos::PCPGIter::iterate(): non-positive value for p^H*A*p encountered!" );

// alpha := <R_,Z_> / <P,AP>
alpha(0,0) = rHz(0,0) / pAp(0,0);
alpha = DMT::ValueConst(*rHz,0,0) / DMT::ValueConst(*pAp,0,0);

// positive alpha required
TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, CGPositiveDefiniteFailure,
TEUCHOS_TEST_FOR_EXCEPTION( alpha <= zero, CGPositiveDefiniteFailure,
"Belos::PCPGIter::iterate(): non-positive value for alpha encountered!" );

// solution update x += alpha * P
if( curDim_ < savedBlocks_ ){
MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
MVT::MvAddMv( one, *cur_soln_vec, alpha, *P, *cur_soln_vec );
}else{
MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
}
//lp_->updateSolution(); ... does nothing.
//
// The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
//
rHz_old(0,0) = rHz(0,0);
rHz_old = DMT::ValueConst(*rHz,0,0);
//
// residual update R_ := R_ - alpha * AP
//
if( prevUdim_ + iter_ <= savedBlocks_ ){
MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
MVT::MvAddMv( one, *R_, -alpha, *AP, *R_ );
AP = Teuchos::null;
}else{
MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
}
//
// update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
Expand All @@ -768,23 +765,19 @@ namespace Belos {
Z_ = R_;
}
//
MVT::MvTransMv( one, *R_, *Z_, rHz );
MVT::MvTransMv( one, *R_, *Z_, *rHz );
DMT::SyncDeviceToHost( *rHz );
//
beta(0,0) = rHz(0,0) / rHz_old(0,0);
beta = DMT::ValueConst(*rHz,0,0) / rHz_old;
//
if( curDim_ < savedBlocks_ ){
curDim_++; // update basis dim
curind[0] = curDim_ - 1;
Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
MVT::MvAddMv( one, *Z_, beta, *P, *Pnext );
if( prevUdim_ ){ // Deflate seed space
MVT::MvTransMv( one, *Cprev, *Z_, CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
if( debug ){
std::cout << " Check CZ " << std::endl;
MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
CZ.print( std::cout );
}
MVT::MvTransMv( one, *Cprev, *Z_, *CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *Pnext ); // Pnext -= U*(C'Z)
}
P = Teuchos::null;
if( curDim_ == savedBlocks_ ){
Expand All @@ -794,16 +787,10 @@ namespace Belos {
}
Pnext = Teuchos::null;
}else{
MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
if( prevUdim_ ){ // Deflate seed space
MVT::MvTransMv( one, *Cprev, *Z_, CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ ); // P_ -= U*(C'Z)

if( debug ){
std::cout << " Check CZ " << std::endl;
MVT::MvTransMv( one, *Cprev, *P_, CZ );
CZ.print( std::cout );
}
MVT::MvTransMv( one, *Cprev, *Z_, *CZ );
MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *P_ ); // P_ -= U*(C'Z)
}
}
// CGB: 5/26/2010
Expand Down
Loading

0 comments on commit 801ff16

Please sign in to comment.