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 );
527 Z_ = MVT::Clone( *tmp, 1 );
529 if (P_ == Teuchos::null) {
530 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
531 P_ = MVT::Clone( *tmp, 1 );
533 if (AP_ == Teuchos::null) {
534 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
535 AP_ = MVT::Clone( *tmp, 1 );
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.");
546 C_ = MVT::Clone( *tmp, savedBlocks_ );
550 if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
551 Teuchos::RCP<const MV> tmp = C_;
552 C_ = MVT::Clone( *tmp, savedBlocks_ );
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.");
559 U_ = MVT::Clone( *tmp, savedBlocks_ );
563 if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
564 Teuchos::RCP<const MV> tmp = U_;
565 U_ = MVT::Clone( *tmp, savedBlocks_ );
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)
614 Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.
U, index );
615 newstate.
C = MVT::Clone( *newstate.
U, prevUdim_ );
616 Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.
C, index );
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_ );
638 MVT::SetBlock( *Z_, zero_index , *P_ );
641 MVT::SetBlock( *R_, zero_index, *P_ );
644 std::vector<int> nextind(1);
645 nextind[0] = curDim_;
647 MVT::SetBlock( *P_, nextind, *newstate.
U );
650 newstate.
curDim = curDim_;
652 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
U) != savedBlocks_ ,
653 std::invalid_argument, errstr );
654 if (newstate.
U != U_) {
658 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
C) != savedBlocks_ ,
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 );
707 Uprev = MVT::CloneView(*U_, prevInd);
708 Cprev = MVT::CloneView(*C_, prevInd);
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);
728 std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
731 curind[0] = curDim_ - 1;
732 P = MVT::CloneViewNonConst(*U_,curind);
733 MVT::MvTransMv( one, *Cprev, *P, CZ );
734 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );
737 MVT::MvTransMv( one, *Cprev, *P, CZ );
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);
744 MVT::SetBlock( *P, zero_index, *P_ );
750 MVT::MvTransMv( one, *R_, *Z_, rHz );
755 while (stest_->checkStatus(
this) !=
Passed ) {
756 Teuchos::RCP<const MV> P;
760 curind[0] = curDim_ - 1;
762 MVT::MvNorm(*R_, rnorm);
763 std::cout << iter_ <<
" " << curDim_ <<
" " << rnorm[0] << std::endl;
765 if( prevUdim_ + iter_ < savedBlocks_ ){
766 P = MVT::CloneView(*U_,curind);
767 AP = MVT::CloneViewNonConst(*C_,curind);
768 lp_->applyOp( *P, *AP );
769 MVT::MvTransMv( one, *P, *AP, pAp );
771 if( prevUdim_ + iter_ == savedBlocks_ ){
772 AP = MVT::CloneViewNonConst(*C_,curind);
773 lp_->applyOp( *P_, *AP );
774 MVT::MvTransMv( one, *P_, *AP, pAp );
776 lp_->applyOp( *P_, *AP_ );
777 MVT::MvTransMv( one, *P_, *AP_, pAp );
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_ ){
810 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
813 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
818 if ( lp_->getLeftPrec() != Teuchos::null ) {
819 lp_->applyLeftPrec( *R_, *Z_ );
824 MVT::MvTransMv( one, *R_, *Z_, rHz );
826 beta(0,0) = rHz(0,0) / rHz_old(0,0);
828 if( curDim_ < savedBlocks_ ){
830 curind[0] = curDim_ - 1;
831 Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
832 MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
834 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
835 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext );
837 std::cout <<
" Check CZ " << std::endl;
838 MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
839 CZ.print( std::cout );
843 if( curDim_ == savedBlocks_ ){
844 std::vector<int> zero_index(1);
846 MVT::SetBlock( *Pnext, zero_index, *P_ );
848 Pnext = Teuchos::null;
850 MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
852 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
853 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );
856 std::cout <<
" Check CZ " << std::endl;
857 MVT::MvTransMv( one, *Cprev, *P_, CZ );
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.
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?.
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.
PCPGIterLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
PCPGIterLAPACKFailure(const std::string &what_arg)
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
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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