42 #ifndef BELOS_PCPG_ITER_HPP
43 #define BELOS_PCPG_ITER_HPP
86 template <
class ScalarType,
class MV>
123 R(Teuchos::null),
Z(Teuchos::null),
124 P(Teuchos::null),
AP(Teuchos::null),
125 U(Teuchos::null),
C(Teuchos::null),
132 template<
class ScalarType,
class MV,
class OP>
263 if (!initialized_)
return 0;
269 if (!initialized_)
return 0;
293 "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
297 void setSize(
int savedBlocks );
339 bool stateStorageInitialized_;
384 template<
class ScalarType,
class MV,
class OP>
396 stateStorageInitialized_(false),
397 keepDiagonal_(false),
398 initDiagonal_(false),
406 "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
407 int rb = Teuchos::getParameter<int>(params,
"Saved Blocks");
410 keepDiagonal_ = params.get(
"Keep Diagonal",
false);
413 initDiagonal_ = params.get(
"Initialize Diagonal",
false);
421 template <
class ScalarType,
class MV,
class OP>
427 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument,
"Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
429 if ( savedBlocks_ != savedBlocks) {
430 stateStorageInitialized_ =
false;
431 savedBlocks_ = savedBlocks;
432 initialized_ =
false;
441 template <
class ScalarType,
class MV,
class OP>
444 stateStorageInitialized_ =
false;
445 initialized_ =
false;
453 template <
class ScalarType,
class MV,
class OP>
456 if (!stateStorageInitialized_) {
461 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
468 int newsd = savedBlocks_ ;
473 if (Z_ == Teuchos::null) {
475 Z_ = MVT::Clone( *tmp, 1 );
477 if (P_ == Teuchos::null) {
479 P_ = MVT::Clone( *tmp, 1 );
481 if (AP_ == Teuchos::null) {
483 AP_ = MVT::Clone( *tmp, 1 );
486 if (C_ == Teuchos::null) {
491 "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
493 "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
494 C_ = MVT::Clone( *tmp, savedBlocks_ );
498 if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
500 C_ = MVT::Clone( *tmp, savedBlocks_ );
503 if (U_ == Teuchos::null) {
506 "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
507 U_ = MVT::Clone( *tmp, savedBlocks_ );
511 if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
513 U_ = MVT::Clone( *tmp, savedBlocks_ );
517 if (D_ == Teuchos::null) {
521 D_->shape( newsd, newsd );
524 if (D_->numRows() < newsd || D_->numCols() < newsd) {
525 D_->shapeUninitialized( newsd, newsd );
530 stateStorageInitialized_ =
true;
537 template <
class ScalarType,
class MV,
class OP>
542 "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
546 std::string errstr(
"Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
548 if (newstate.
R != Teuchos::null){
551 if (newstate.
U == Teuchos::null){
557 prevUdim_ = newstate.
curDim;
558 if (newstate.
C == Teuchos::null){
559 std::vector<int> index(prevUdim_);
560 for (
int i=0; i< prevUdim_; ++i)
563 newstate.
C = MVT::Clone( *newstate.
U, prevUdim_ );
565 lp_->apply( *Ukeff, *Ckeff );
567 curDim_ = prevUdim_ ;
571 if (!stateStorageInitialized_)
578 newstate.
curDim = curDim_;
582 std::vector<int> zero_index(1);
584 if ( lp_->getLeftPrec() != Teuchos::null ) {
585 lp_->applyLeftPrec( *R_, *Z_ );
586 MVT::SetBlock( *Z_, zero_index , *P_ );
589 MVT::SetBlock( *R_, zero_index, *P_ );
592 std::vector<int> nextind(1);
593 nextind[0] = curDim_;
595 MVT::SetBlock( *P_, nextind, *newstate.
U );
598 newstate.
curDim = curDim_;
601 std::invalid_argument, errstr );
602 if (newstate.
U != U_) {
607 std::invalid_argument, errstr );
608 if (newstate.
C != C_) {
615 "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
625 template <
class ScalarType,
class MV,
class OP>
631 if (initialized_ ==
false) {
634 const bool debug =
false;
643 std::cout <<
" Iterate Warning: begin from nonzero iter_ ?" << std::endl;
646 std::vector<int> prevInd;
652 prevInd.resize( prevUdim_ );
653 for(
int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
655 Uprev = MVT::CloneView(*U_, prevInd);
656 Cprev = MVT::CloneView(*C_, prevInd);
664 "Belos::PCPGIter::iterate(): current linear system has more than one std::vector!" );
668 "Belos::PCPGIter::iterate(): mistake in initialization !" );
675 std::vector<int> curind(1);
676 std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
679 curind[0] = curDim_ - 1;
680 P = MVT::CloneViewNonConst(*U_,curind);
681 MVT::MvTransMv( one, *Cprev, *P, CZ );
682 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );
685 MVT::MvTransMv( one, *Cprev, *P, CZ );
686 std::cout <<
" Input CZ post ortho " << std::endl;
687 CZ.
print( std::cout );
689 if( curDim_ == savedBlocks_ ){
690 std::vector<int> zero_index(1);
692 MVT::SetBlock( *P, zero_index, *P_ );
698 MVT::MvTransMv( one, *R_, *Z_, rHz );
703 while (stest_->checkStatus(
this) !=
Passed ) {
708 curind[0] = curDim_ - 1;
710 MVT::MvNorm(*R_, rnorm);
711 std::cout << iter_ <<
" " << curDim_ <<
" " << rnorm[0] << std::endl;
713 if( prevUdim_ + iter_ < savedBlocks_ ){
714 P = MVT::CloneView(*U_,curind);
715 AP = MVT::CloneViewNonConst(*C_,curind);
716 lp_->applyOp( *P, *AP );
717 MVT::MvTransMv( one, *P, *AP, pAp );
719 if( prevUdim_ + iter_ == savedBlocks_ ){
720 AP = MVT::CloneViewNonConst(*C_,curind);
721 lp_->applyOp( *P_, *AP );
722 MVT::MvTransMv( one, *P_, *AP, pAp );
724 lp_->applyOp( *P_, *AP_ );
725 MVT::MvTransMv( one, *P_, *AP_, pAp );
729 if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
730 (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
734 "Belos::PCPGIter::iterate(): non-positive value for p^H*A*p encountered!" );
737 alpha(0,0) = rHz(0,0) / pAp(0,0);
741 "Belos::PCPGIter::iterate(): non-positive value for alpha encountered!" );
744 if( curDim_ < savedBlocks_ ){
745 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
747 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
753 rHz_old(0,0) = rHz(0,0);
757 if( prevUdim_ + iter_ <= savedBlocks_ ){
758 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
761 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
766 if ( lp_->getLeftPrec() != Teuchos::null ) {
767 lp_->applyLeftPrec( *R_, *Z_ );
772 MVT::MvTransMv( one, *R_, *Z_, rHz );
774 beta(0,0) = rHz(0,0) / rHz_old(0,0);
776 if( curDim_ < savedBlocks_ ){
778 curind[0] = curDim_ - 1;
780 MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
782 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
783 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext );
785 std::cout <<
" Check CZ " << std::endl;
786 MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
787 CZ.
print( std::cout );
791 if( curDim_ == savedBlocks_ ){
792 std::vector<int> zero_index(1);
794 MVT::SetBlock( *Pnext, zero_index, *P_ );
796 Pnext = Teuchos::null;
798 MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
800 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
801 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );
804 std::cout <<
" Check CZ " << std::endl;
805 MVT::MvTransMv( one, *Cprev, *P_, CZ );
806 CZ.
print( std::cout );
815 TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error,
"Loop recurrence violated. Please contact Belos team.");
817 if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_;
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 which augments the basic interface for a conjugate gradient linear solver ite...
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
Declaration of basic traits for the multivector type.
virtual std::ostream & print(std::ostream &os) const
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...
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.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero...
bool isParameter(const std::string &name) const
virtual ~PCPGIter()
Destructor.
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int getNumIters() const
Get the current iteration count.
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.
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
Teuchos::RCP< MV > U
The recycled subspace.
void resetState()
tell the Iterator to "reset" itself; delete and rebuild the seed space.
int reshape(OrdinalType numRows, OrdinalType numCols)
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.
Teuchos::RCP< MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Class which defines basic traits for the operator type.
CGIterationInitFailure is thrown when the CGIteration object is unable to generate an initial iterate...
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...
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?.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType > SCT