42 #ifndef BELOS_PCPG_ITER_HPP 43 #define BELOS_PCPG_ITER_HPP 59 #include "Teuchos_SerialDenseMatrix.hpp" 60 #include "Teuchos_SerialDenseVector.hpp" 61 #include "Teuchos_ScalarTraits.hpp" 62 #include "Teuchos_ParameterList.hpp" 63 #include "Teuchos_TimeMonitor.hpp" 85 template <
class ScalarType,
class MV>
118 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
D;
122 R(Teuchos::null),
Z(Teuchos::null),
123 P(Teuchos::null),
AP(Teuchos::null),
124 U(Teuchos::null),
C(Teuchos::null),
167 template<
class ScalarType,
class MV,
class OP>
177 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
194 Teuchos::ParameterList ¶ms );
298 if (!initialized_)
return 0;
304 if (!initialized_)
return 0;
327 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
328 "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
332 void setSize(
int savedBlocks );
353 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
354 const Teuchos::RCP<OutputManager<ScalarType> > om_;
355 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
356 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
374 bool stateStorageInitialized_;
404 Teuchos::RCP<MV> AP_;
414 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
419 template<
class ScalarType,
class MV,
class OP>
424 Teuchos::ParameterList ¶ms ):
431 stateStorageInitialized_(false),
432 keepDiagonal_(false),
433 initDiagonal_(false),
440 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Saved Blocks"), std::invalid_argument,
441 "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
442 int rb = Teuchos::getParameter<int>(params,
"Saved Blocks");
445 keepDiagonal_ = params.get(
"Keep Diagonal",
false);
448 initDiagonal_ = params.get(
"Initialize Diagonal",
false);
456 template <
class ScalarType,
class MV,
class OP>
462 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument,
"Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
464 if ( savedBlocks_ != savedBlocks) {
465 stateStorageInitialized_ =
false;
466 savedBlocks_ = savedBlocks;
467 initialized_ =
false;
476 template <
class ScalarType,
class MV,
class OP>
479 stateStorageInitialized_ =
false;
480 initialized_ =
false;
488 template <
class ScalarType,
class MV,
class OP>
491 if (!stateStorageInitialized_) {
494 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
495 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
496 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
503 int newsd = savedBlocks_ ;
508 if (Z_ == Teuchos::null) {
509 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
510 Z_ = MVT::Clone( *tmp, 1 );
512 if (P_ == Teuchos::null) {
513 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
514 P_ = MVT::Clone( *tmp, 1 );
516 if (AP_ == Teuchos::null) {
517 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
518 AP_ = MVT::Clone( *tmp, 1 );
521 if (C_ == Teuchos::null) {
524 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
525 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
526 "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
527 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
528 "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
529 C_ = MVT::Clone( *tmp, savedBlocks_ );
533 if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
534 Teuchos::RCP<const MV> tmp = C_;
535 C_ = MVT::Clone( *tmp, savedBlocks_ );
538 if (U_ == Teuchos::null) {
539 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
540 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
541 "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
542 U_ = MVT::Clone( *tmp, savedBlocks_ );
546 if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
547 Teuchos::RCP<const MV> tmp = U_;
548 U_ = MVT::Clone( *tmp, savedBlocks_ );
552 if (D_ == Teuchos::null) {
553 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
556 D_->shape( newsd, newsd );
559 if (D_->numRows() < newsd || D_->numCols() < newsd) {
560 D_->shapeUninitialized( newsd, newsd );
565 stateStorageInitialized_ =
true;
572 template <
class ScalarType,
class MV,
class OP>
576 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
577 "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
581 std::string errstr(
"Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
583 if (newstate.
R != Teuchos::null){
586 if (newstate.
U == Teuchos::null){
592 prevUdim_ = newstate.
curDim;
593 if (newstate.
C == Teuchos::null){
594 std::vector<int> index(prevUdim_);
595 for (
int i=0; i< prevUdim_; ++i)
597 Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.
U, index );
598 newstate.
C = MVT::Clone( *newstate.
U, prevUdim_ );
599 Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.
C, index );
600 lp_->apply( *Ukeff, *Ckeff );
602 curDim_ = prevUdim_ ;
606 if (!stateStorageInitialized_)
613 newstate.
curDim = curDim_;
617 std::vector<int> zero_index(1);
619 if ( lp_->getLeftPrec() != Teuchos::null ) {
620 lp_->applyLeftPrec( *R_, *Z_ );
621 MVT::SetBlock( *Z_, zero_index , *P_ );
624 MVT::SetBlock( *R_, zero_index, *P_ );
627 std::vector<int> nextind(1);
628 nextind[0] = curDim_;
630 MVT::SetBlock( *P_, nextind, *newstate.
U );
633 newstate.
curDim = curDim_;
635 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
U) != savedBlocks_ ,
636 std::invalid_argument, errstr );
637 if (newstate.
U != U_) {
641 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
C) != savedBlocks_ ,
642 std::invalid_argument, errstr );
643 if (newstate.
C != C_) {
649 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
650 "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
660 template <
class ScalarType,
class MV,
class OP>
666 if (initialized_ ==
false) {
669 const bool debug =
false;
672 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
673 Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
674 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
675 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
678 std::cout <<
" Iterate Warning: begin from nonzero iter_ ?" << std::endl;
681 std::vector<int> prevInd;
682 Teuchos::RCP<const MV> Uprev;
683 Teuchos::RCP<const MV> Cprev;
684 Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
687 prevInd.resize( prevUdim_ );
688 for(
int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
689 CZ.reshape( prevUdim_ , 1 );
690 Uprev = MVT::CloneView(*U_, prevInd);
691 Cprev = MVT::CloneView(*C_, prevInd);
695 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
699 "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
703 "Belos::CGIter::iterate(): mistake in initialization !" );
706 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
707 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
710 std::vector<int> curind(1);
711 std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
714 curind[0] = curDim_ - 1;
715 P = MVT::CloneViewNonConst(*U_,curind);
716 MVT::MvTransMv( one, *Cprev, *P, CZ );
717 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );
720 MVT::MvTransMv( one, *Cprev, *P, CZ );
721 std::cout <<
" Input CZ post ortho " << std::endl;
722 CZ.print( std::cout );
724 if( curDim_ == savedBlocks_ ){
725 std::vector<int> zero_index(1);
727 MVT::SetBlock( *P, zero_index, *P_ );
733 MVT::MvTransMv( one, *R_, *Z_, rHz );
738 while (stest_->checkStatus(
this) !=
Passed ) {
739 Teuchos::RCP<const MV> P;
743 curind[0] = curDim_ - 1;
745 MVT::MvNorm(*R_, rnorm);
746 std::cout << iter_ <<
" " << curDim_ <<
" " << rnorm[0] << std::endl;
748 if( prevUdim_ + iter_ < savedBlocks_ ){
749 P = MVT::CloneView(*U_,curind);
750 AP = MVT::CloneViewNonConst(*C_,curind);
751 lp_->applyOp( *P, *AP );
752 MVT::MvTransMv( one, *P, *AP, pAp );
754 if( prevUdim_ + iter_ == savedBlocks_ ){
755 AP = MVT::CloneViewNonConst(*C_,curind);
756 lp_->applyOp( *P_, *AP );
757 MVT::MvTransMv( one, *P_, *AP, pAp );
759 lp_->applyOp( *P_, *AP_ );
760 MVT::MvTransMv( one, *P_, *AP_, pAp );
764 if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
765 (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
769 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
772 alpha(0,0) = rHz(0,0) / pAp(0,0);
776 "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
779 if( curDim_ < savedBlocks_ ){
780 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
782 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
788 rHz_old(0,0) = rHz(0,0);
792 if( prevUdim_ + iter_ <= savedBlocks_ ){
793 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
796 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
801 if ( lp_->getLeftPrec() != Teuchos::null ) {
802 lp_->applyLeftPrec( *R_, *Z_ );
807 MVT::MvTransMv( one, *R_, *Z_, rHz );
809 beta(0,0) = rHz(0,0) / rHz_old(0,0);
811 if( curDim_ < savedBlocks_ ){
813 curind[0] = curDim_ - 1;
814 Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
815 MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
817 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
818 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext );
820 std::cout <<
" Check CZ " << std::endl;
821 MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
822 CZ.print( std::cout );
826 if( curDim_ == savedBlocks_ ){
827 std::vector<int> zero_index(1);
829 MVT::SetBlock( *Pnext, zero_index, *P_ );
831 Pnext = Teuchos::null;
833 MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
835 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
836 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );
839 std::cout <<
" Check CZ " << std::endl;
840 MVT::MvTransMv( one, *Cprev, *P_, CZ );
841 CZ.print( std::cout );
850 TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error,
"Loop recurrence violated. Please contact Belos team.");
852 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.
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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
A pure virtual class for defining the status tests for the Belos iterative solvers.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
virtual ~PCPGIter()
Destructor.
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
int getNumIters() const
Get the current iteration count.
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)
SCT::magnitudeType MagnitudeType
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > U
The recycled subspace.
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem...
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.
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
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
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