42 #ifndef BELOS_RCG_ITER_HPP 43 #define BELOS_RCG_ITER_HPP 91 template <
class ScalarType,
class MV>
195 template<
class ScalarType,
class MV,
class OP>
329 "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
333 void setSize(
int recycleBlocks,
int numBlocks );
411 template<
class ScalarType,
class MV,
class OP>
428 "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
429 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
432 "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
433 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
441 template <
class ScalarType,
class MV,
class OP>
445 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
446 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
447 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument,
"Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
449 numBlocks_ = numBlocks;
450 recycleBlocks_ = recycleBlocks;
456 template <
class ScalarType,
class MV,
class OP>
460 if (newstate.
P != Teuchos::null &&
461 newstate.
Ap != Teuchos::null &&
462 newstate.
r != Teuchos::null &&
463 newstate.
z != Teuchos::null &&
464 newstate.
U != Teuchos::null &&
465 newstate.
AU != Teuchos::null &&
466 newstate.
Alpha != Teuchos::null &&
467 newstate.
Beta != Teuchos::null &&
468 newstate.
D != Teuchos::null &&
469 newstate.
Delta != Teuchos::null &&
470 newstate.
LUUTAU != Teuchos::null &&
471 newstate.
ipiv != Teuchos::null &&
472 newstate.
rTz_old != Teuchos::null) {
474 curDim_ = newstate.
curDim;
479 existU_ = newstate.
existU;
482 Alpha_ = newstate.
Alpha;
483 Beta_ = newstate.
Beta;
485 Delta_ = newstate.
Delta;
486 LUUTAU_ = newstate.
LUUTAU;
487 ipiv_ = newstate.
ipiv;
493 "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
496 "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
499 "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
502 "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
505 "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
508 "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
511 "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
514 "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
517 "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
520 "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
523 "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
526 "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
529 "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
540 template <
class ScalarType,
class MV,
class OP>
544 "Belos::RCGIter::iterate(): RCGIter class not initialized." );
554 std::vector<int> index(1);
562 "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
565 int searchDim = numBlocks_+1;
577 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= searchDim) {
582 p_ = MVT::CloneView( *P_, index );
583 lp_->applyOp( *p_, *Ap_ );
586 MVT::MvTransMv( one, *p_, *Ap_, pAp );
587 (*D_)(i_,0) = pAp(0,0);
590 (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
596 MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
597 lp_->updateSolution();
600 MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
602 std::vector<MagnitudeType> norm(1);
603 MVT::MvNorm( *r_, norm );
607 if ( lp_->getLeftPrec() != Teuchos::null ) {
608 lp_->applyLeftPrec( *r_, *z_ );
610 else if ( lp_->getRightPrec() != Teuchos::null ) {
611 lp_->applyRightPrec( *r_, *z_ );
618 MVT::MvTransMv( one, *r_, *z_, rTz );
621 (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
624 (*rTz_old_)(0,0) = rTz(0,0);
629 pnext_ = MVT::CloneViewNonConst( *P_, index );
634 MVT::MvTransMv( one, *AU_, *z_, mu );
637 lapack.
GETRS(
TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.
values(), mu.
stride(), &info );
639 "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
642 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
644 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
648 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
653 pnext_ = Teuchos::null;
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
ScalarType * values() const
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
RCGIterOrthoFailure(const std::string &what_arg)
RCGIterFailure is thrown when the RCGIter object is unable to compute the next iterate in the RCGIter...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
RCGIterOrthoFailure is thrown when the RCGIter object is unable to compute independent direction vect...
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
Declaration of basic traits for the multivector type.
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
OperatorTraits< ScalarType, MV, OP > OPT
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
MultiVecTraits< ScalarType, MV > MVT
virtual ~RCGIter()
Destructor.
RCGIterLAPACKFailure(const std::string &what_arg)
int getNumIters() const
Get the current iteration count.
Traits class which defines basic operations on multivectors.
void setSize(int recycleBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Structure to contain pointers to RCGIter state variables.
bool existU
Flag to indicate the recycle space should be used.
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
RCGIterInitFailure(const std::string &what_arg)
RCGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList ¶ms)
RCGIter constructor with linear problem, solver utilities, and parameter list of solver options...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
RCGIterFailure(const std::string &what_arg)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
RCGIterInitFailure is thrown when the RCGIter object is unable to generate an initial iterate in the ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
int curDim
The current dimension of the reduction.
void setRecycledBlocks(int recycleBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
bool isParameter(const std::string &name) const
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
Class which defines basic traits for the operator type.
OrdinalType stride() const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
void iterate()
This method performs RCG iterations until the status test indicates the need to stop or an error occu...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Parent class to all Belos exceptions.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< std::vector< int > > ipiv_
RCGIterLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
void setBlockSize(int blockSize)
Set the blocksize.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::ScalarTraits< ScalarType > SCT
void GETRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
SCT::magnitudeType MagnitudeType