42#ifndef BELOS_PCPG_ITER_HPP
43#define BELOS_PCPG_ITER_HPP
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
86 template <
class ScalarType,
class MV>
119 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
D;
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>
142 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
159 Teuchos::ParameterList ¶ms );
263 if (!initialized_)
return 0;
269 if (!initialized_)
return 0;
292 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
293 "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
297 void setSize(
int savedBlocks );
318 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
319 const Teuchos::RCP<OutputManager<ScalarType> > om_;
320 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
321 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
339 bool stateStorageInitialized_;
369 Teuchos::RCP<MV> AP_;
379 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
384 template<
class ScalarType,
class MV,
class OP>
389 Teuchos::ParameterList ¶ms ):
396 stateStorageInitialized_(false),
397 keepDiagonal_(false),
398 initDiagonal_(false),
405 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Saved Blocks"), std::invalid_argument,
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_) {
459 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
460 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
461 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
468 int newsd = savedBlocks_ ;
473 if (Z_ == Teuchos::null) {
474 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
475 Z_ = MVT::Clone( *tmp, 1 );
477 if (P_ == Teuchos::null) {
478 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
479 P_ = MVT::Clone( *tmp, 1 );
481 if (AP_ == Teuchos::null) {
482 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
483 AP_ = MVT::Clone( *tmp, 1 );
486 if (C_ == Teuchos::null) {
489 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
490 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
491 "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
492 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
493 "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
494 C_ = MVT::Clone( *tmp, savedBlocks_ );
498 if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
499 Teuchos::RCP<const MV> tmp = C_;
500 C_ = MVT::Clone( *tmp, savedBlocks_ );
503 if (U_ == Teuchos::null) {
504 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
505 TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
506 "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
507 U_ = MVT::Clone( *tmp, savedBlocks_ );
511 if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
512 Teuchos::RCP<const MV> tmp = U_;
513 U_ = MVT::Clone( *tmp, savedBlocks_ );
517 if (D_ == Teuchos::null) {
518 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
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>
541 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
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)
562 Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.
U, index );
563 newstate.
C = MVT::Clone( *newstate.
U, prevUdim_ );
564 Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.
C, index );
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_;
600 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
U) != savedBlocks_ ,
601 std::invalid_argument, errstr );
602 if (newstate.
U != U_) {
606 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
C) != savedBlocks_ ,
607 std::invalid_argument, errstr );
608 if (newstate.
C != C_) {
614 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
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;
637 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
638 Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
639 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
640 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
643 std::cout <<
" Iterate Warning: begin from nonzero iter_ ?" << std::endl;
646 std::vector<int> prevInd;
647 Teuchos::RCP<const MV> Uprev;
648 Teuchos::RCP<const MV> Cprev;
649 Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
652 prevInd.resize( prevUdim_ );
653 for(
int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
654 CZ.reshape( prevUdim_ , 1 );
655 Uprev = MVT::CloneView(*U_, prevInd);
656 Cprev = MVT::CloneView(*C_, prevInd);
660 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
664 "Belos::PCPGIter::iterate(): current linear system has more than one std::vector!" );
668 "Belos::PCPGIter::iterate(): mistake in initialization !" );
671 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
672 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
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 ) {
704 Teuchos::RCP<const MV> P;
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;
779 Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
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_;
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
CGIterationInitFailure is thrown when the CGIteration object is unable to generate an initial iterate...
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed....
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
MultiVecTraits< ScalarType, MV > MVT
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
void resetNumIters(int iter=0)
Reset the iteration count.
virtual ~PCPGIter()
Destructor.
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
void resetState()
tell the Iterator to "reset" itself; delete and rebuild the seed space.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
void setSize(int savedBlocks)
Set the maximum number of saved or recycled blocks used by the iterative solver.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
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.
void iterate()
PCPGIter iterates CG until the status test either requests a stop or detects an error....
bool isInitialized()
States whether the solver has been initialized or not.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PCPGIter state variables.
Teuchos::RCP< MV > AP
The matrix A applied to current decent direction std::vector.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Teuchos::RCP< MV > Z
The current preconditioned residual.
Teuchos::RCP< MV > U
The recycled subspace.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > P
The current decent direction std::vector.
Teuchos::RCP< MV > R
The current residual.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.