42 #ifndef BELOS_RCG_ITER_HPP
43 #define BELOS_RCG_ITER_HPP
91 template <
class ScalarType,
class MV>
136 U(Teuchos::null),
AU(Teuchos::null),
137 Alpha(Teuchos::null),
Beta(Teuchos::null),
D(Teuchos::null),
rTz_old(Teuchos::null),
144 template<
class ScalarType,
class MV,
class OP>
278 "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
282 void setSize(
int recycleBlocks,
int numBlocks );
360 template<
class ScalarType,
class MV,
class OP>
377 "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
378 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
381 "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
382 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
390 template <
class ScalarType,
class MV,
class OP>
394 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
395 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
396 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument,
"Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
398 numBlocks_ = numBlocks;
399 recycleBlocks_ = recycleBlocks;
405 template <
class ScalarType,
class MV,
class OP>
423 curDim_ = newstate.
curDim;
428 existU_ = newstate.
existU;
431 Alpha_ = newstate.
Alpha;
432 Beta_ = newstate.
Beta;
434 Delta_ = newstate.
Delta;
435 LUUTAU_ = newstate.
LUUTAU;
436 ipiv_ = newstate.
ipiv;
442 "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
445 "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
448 "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
451 "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
454 "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
457 "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
460 "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
463 "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
466 "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
469 "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
472 "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
475 "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
478 "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
489 template <
class ScalarType,
class MV,
class OP>
493 "Belos::RCGIter::iterate(): RCGIter class not initialized." );
503 std::vector<int> index(1);
511 "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
514 int searchDim = numBlocks_+1;
526 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= searchDim) {
531 p_ = MVT::CloneView( *P_, index );
532 lp_->applyOp( *p_, *Ap_ );
535 MVT::MvTransMv( one, *p_, *Ap_, pAp );
536 (*D_)(i_,0) = pAp(0,0);
539 (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
543 "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
546 MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
547 lp_->updateSolution();
550 MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
552 std::vector<MagnitudeType> norm(1);
553 MVT::MvNorm( *r_, norm );
558 lp_->applyLeftPrec( *r_, *z_ );
561 lp_->applyRightPrec( *r_, *z_ );
568 MVT::MvTransMv( one, *r_, *z_, rTz );
571 (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
574 (*rTz_old_)(0,0) = rTz(0,0);
579 pnext_ = MVT::CloneViewNonConst( *P_, index );
584 MVT::MvTransMv( one, *AU_, *z_, mu );
587 lapack.
GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.
values(), mu.
stride(), &info );
589 "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
592 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
594 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
598 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
const Teuchos::RCP< OutputManager< ScalarType > > om_
ScalarType * values() const
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.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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...
int getRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
OperatorTraits< ScalarType, MV, OP > OPT
int getNumIters() const
Get the current iteration count.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
MultiVecTraits< ScalarType, MV > MVT
virtual ~RCGIter()
Destructor.
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...
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero...
bool isParameter(const std::string &name) const
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Structure to contain pointers to RCGIter state variables.
CGIterationLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
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.
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
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
int curDim
The current dimension of the reduction.
void setRecycledBlocks(int recycleBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Class which defines basic traits for the operator type.
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...
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< std::vector< int > > ipiv_
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::ScalarTraits< ScalarType > SCT
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
OrdinalType stride() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
SCT::magnitudeType MagnitudeType