42 #ifndef BELOS_PCPG_ITER_HPP 43 #define BELOS_PCPG_ITER_HPP 59 #include "Teuchos_BLAS.hpp" 60 #include "Teuchos_LAPACK.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_ScalarTraits.hpp" 64 #include "Teuchos_ParameterList.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 87 template <
class ScalarType,
class MV>
120 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
D;
124 R(Teuchos::null), Z(Teuchos::null),
125 P(Teuchos::null), AP(Teuchos::null),
126 U(Teuchos::null), C(Teuchos::null),
184 template<
class ScalarType,
class MV,
class OP>
194 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
211 Teuchos::ParameterList ¶ms );
315 if (!initialized_)
return 0;
321 if (!initialized_)
return 0;
344 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
345 "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
349 void setSize(
int savedBlocks );
370 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
371 const Teuchos::RCP<OutputManager<ScalarType> > om_;
372 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
373 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
391 bool stateStorageInitialized_;
421 Teuchos::RCP<MV> AP_;
431 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
436 template<
class ScalarType,
class MV,
class OP>
441 Teuchos::ParameterList ¶ms ):
448 stateStorageInitialized_(false),
449 keepDiagonal_(false),
450 initDiagonal_(false),
457 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Saved Blocks"), std::invalid_argument,
458 "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
459 int rb = Teuchos::getParameter<int>(params,
"Saved Blocks");
462 keepDiagonal_ = params.get(
"Keep Diagonal",
false);
465 initDiagonal_ = params.get(
"Initialize Diagonal",
false);
473 template <
class ScalarType,
class MV,
class OP>
479 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument,
"Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
481 if ( savedBlocks_ != savedBlocks) {
482 stateStorageInitialized_ =
false;
483 savedBlocks_ = savedBlocks;
484 initialized_ =
false;
493 template <
class ScalarType,
class MV,
class OP>
496 stateStorageInitialized_ =
false;
497 initialized_ =
false;
505 template <
class ScalarType,
class MV,
class OP>
508 if (!stateStorageInitialized_) {
511 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
512 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
513 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
520 int newsd = savedBlocks_ ;
525 if (Z_ == Teuchos::null) {
526 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
529 if (P_ == Teuchos::null) {
530 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
533 if (AP_ == Teuchos::null) {
534 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
538 if (C_ == Teuchos::null) {
541 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
542 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
543 "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
544 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
545 "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
551 Teuchos::RCP<const MV> tmp = C_;
555 if (U_ == Teuchos::null) {
556 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
557 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
558 "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
564 Teuchos::RCP<const MV> tmp = U_;
569 if (D_ == Teuchos::null) {
570 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
573 D_->shape( newsd, newsd );
576 if (D_->numRows() < newsd || D_->numCols() < newsd) {
577 D_->shapeUninitialized( newsd, newsd );
582 stateStorageInitialized_ =
true;
589 template <
class ScalarType,
class MV,
class OP>
593 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
594 "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
598 std::string errstr(
"Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
600 if (newstate.
R != Teuchos::null){
603 if (newstate.
U == Teuchos::null){
609 prevUdim_ = newstate.
curDim;
610 if (newstate.
C == Teuchos::null){
611 std::vector<int> index(prevUdim_);
612 for (
int i=0; i< prevUdim_; ++i)
617 lp_->apply( *Ukeff, *Ckeff );
619 curDim_ = prevUdim_ ;
623 if (!stateStorageInitialized_)
630 newstate.
curDim = curDim_;
634 std::vector<int> zero_index(1);
636 if ( lp_->getLeftPrec() != Teuchos::null ) {
637 lp_->applyLeftPrec( *R_, *Z_ );
644 std::vector<int> nextind(1);
645 nextind[0] = curDim_;
650 newstate.
curDim = curDim_;
653 std::invalid_argument, errstr );
654 if (newstate.
U != U_) {
659 std::invalid_argument, errstr );
660 if (newstate.
C != C_) {
666 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
667 "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
677 template <
class ScalarType,
class MV,
class OP>
683 if (initialized_ ==
false) {
686 const bool debug =
false;
689 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
690 Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
691 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
692 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
695 std::cout <<
" Iterate Warning: begin from nonzero iter_ ?" << std::endl;
698 std::vector<int> prevInd;
699 Teuchos::RCP<const MV> Uprev;
700 Teuchos::RCP<const MV> Cprev;
701 Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
704 prevInd.resize( prevUdim_ );
705 for(
int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
706 CZ.reshape( prevUdim_ , 1 );
712 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
716 "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
720 "Belos::CGIter::iterate(): mistake in initialization !" );
723 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
724 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
727 std::vector<int> curind(1);
731 curind[0] = curDim_ - 1;
738 std::cout <<
" Input CZ post ortho " << std::endl;
739 CZ.print( std::cout );
741 if( curDim_ == savedBlocks_ ){
742 std::vector<int> zero_index(1);
755 while (stest_->checkStatus(
this) !=
Passed ) {
756 Teuchos::RCP<const MV> P;
760 curind[0] = curDim_ - 1;
763 std::cout << iter_ <<
" " << curDim_ <<
" " << rnorm[0] << std::endl;
765 if( prevUdim_ + iter_ < savedBlocks_ ){
768 lp_->applyOp( *P, *AP );
771 if( prevUdim_ + iter_ == savedBlocks_ ){
773 lp_->applyOp( *P_, *AP );
776 lp_->applyOp( *P_, *AP_ );
781 if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
782 (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
786 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
789 alpha(0,0) = rHz(0,0) / pAp(0,0);
793 "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
796 if( curDim_ < savedBlocks_ ){
797 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
799 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
805 rHz_old(0,0) = rHz(0,0);
809 if( prevUdim_ + iter_ <= savedBlocks_ ){
818 if ( lp_->getLeftPrec() != Teuchos::null ) {
819 lp_->applyLeftPrec( *R_, *Z_ );
826 beta(0,0) = rHz(0,0) / rHz_old(0,0);
828 if( curDim_ < savedBlocks_ ){
830 curind[0] = curDim_ - 1;
837 std::cout <<
" Check CZ " << std::endl;
839 CZ.print( std::cout );
843 if( curDim_ == savedBlocks_ ){
844 std::vector<int> zero_index(1);
848 Pnext = Teuchos::null;
856 std::cout <<
" Check CZ " << std::endl;
858 CZ.print( std::cout );
867 TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error,
"Loop recurrence violated. Please contact Belos team.");
869 if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_;
PCPGIterateFailure is thrown when the PCPGIter object breaks down.
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
bool isInitialized()
States whether the solver has been initialized or not.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< MV > R
The current residual.
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
Declaration of basic traits for the multivector type.
void setSize(int savedBlocks)
Set the maximum number of saved or recycled blocks used by the iterative solver.
int curDim
The current dimension of the reduction.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< MV > P
The current decent direction std::vector.
Teuchos::RCP< MV > AP
The matrix A applied to current decent direction std::vector.
Structure to contain pointers to PCPGIter state variables.
OperatorTraits< ScalarType, MV, OP > OPT
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem...
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
A pure virtual class for defining the status tests for the Belos iterative solvers.
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Class which defines basic traits for the operator type.
PCPGIterLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
PCPGIterLAPACKFailure(const std::string &what_arg)
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
virtual ~PCPGIter()
Destructor.
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
PCPGIterOrthoFailure(const std::string &what_arg)
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SCT::magnitudeType MagnitudeType
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Teuchos::RCP< MV > U
The recycled subspace.
void resetState()
tell the Iterator to "reset" itself; delete and rebuild the seed space.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
void iterate()
PCPGIter iterates CG until the status test either requests a stop or detects an error. In the latter case, std::exception is thrown.
PCPGIterateFailure(const std::string &what_arg)
Teuchos::RCP< MV > Z
The current preconditioned residual.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Class which defines basic traits for the operator type.
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, Teuchos::ParameterList ¶ms)
PCPGIter constructor with linear problem, solver utilities, and parameter list of solver options...
Parent class to all Belos exceptions.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed...
Belos header file which uses auto-configuration information to include necessary C++ headers...
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
PCPGIterInitFailure is thrown when the PCPGIter object is unable to generate an initial iterate in th...
PCPGIterInitFailure(const std::string &what_arg)
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType > SCT