42 #ifndef BELOS_RCG_ITER_HPP 43 #define BELOS_RCG_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" 91 template <
class ScalarType,
class MV>
119 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Alpha;
120 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Beta;
121 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
D;
122 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
rTz_old;
125 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Delta;
128 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
LUUTAU;
130 Teuchos::RCP<std::vector<int> >
ipiv;
133 RCGIterState() : curDim(0), P(Teuchos::null), Ap(Teuchos::null), r(Teuchos::null),
136 U(Teuchos::null), AU(Teuchos::null),
137 Alpha(Teuchos::null), Beta(Teuchos::null), D(Teuchos::null), rTz_old(Teuchos::null),
138 Delta(Teuchos::null), LUUTAU(Teuchos::null), ipiv(Teuchos::null)
195 template<
class ScalarType,
class MV,
class OP>
205 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
221 Teuchos::ParameterList ¶ms );
295 if (!initialized_)
return 0;
315 void setNumBlocks(
int numBlocks) { setSize( recycleBlocks_, numBlocks ); };
328 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
329 "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
333 void setSize(
int recycleBlocks,
int numBlocks );
349 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
350 const Teuchos::RCP<OutputManager<ScalarType> > om_;
351 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
380 Teuchos::RCP<MV> Ap_;
391 Teuchos::RCP<MV> U_, AU_;
394 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_,Beta_,D_;
397 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
400 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
403 Teuchos::RCP<std::vector<int> > ipiv_;
406 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
411 template<
class ScalarType,
class MV,
class OP>
415 Teuchos::ParameterList ¶ms ):
427 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
428 "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
429 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
431 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
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;
492 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
P == Teuchos::null,std::invalid_argument,
493 "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
495 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Ap == Teuchos::null,std::invalid_argument,
496 "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
498 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
r == Teuchos::null,std::invalid_argument,
499 "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
501 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
z == Teuchos::null,std::invalid_argument,
502 "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
504 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
U == Teuchos::null,std::invalid_argument,
505 "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
507 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
AU == Teuchos::null,std::invalid_argument,
508 "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
510 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Alpha == Teuchos::null,std::invalid_argument,
511 "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
513 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Beta == Teuchos::null,std::invalid_argument,
514 "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
516 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
D == Teuchos::null,std::invalid_argument,
517 "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
519 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Delta == Teuchos::null,std::invalid_argument,
520 "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
522 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
LUUTAU == Teuchos::null,std::invalid_argument,
523 "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
525 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
ipiv == Teuchos::null,std::invalid_argument,
526 "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
528 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
rTz_old == Teuchos::null,std::invalid_argument,
529 "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
540 template <
class ScalarType,
class MV,
class OP>
543 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
RCGIterFailure,
544 "Belos::RCGIter::iterate(): RCGIter class not initialized." );
547 Teuchos::LAPACK<int,ScalarType> lapack;
550 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
551 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
554 std::vector<int> index(1);
555 Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
558 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
562 "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
565 int searchDim = numBlocks_+1;
575 Teuchos::RCP<const MV> p_ = Teuchos::null;
576 Teuchos::RCP<MV> pnext_ = Teuchos::null;
577 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= searchDim) {
583 lp_->applyOp( *p_, *Ap_ );
587 (*D_)(i_,0) = pAp(0,0);
590 (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
593 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero,
RCGIterFailure,
"Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
596 MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
597 lp_->updateSolution();
602 std::vector<MagnitudeType> norm(1);
607 if ( lp_->getLeftPrec() != Teuchos::null ) {
608 lp_->applyLeftPrec( *r_, *z_ );
610 else if ( lp_->getRightPrec() != Teuchos::null ) {
611 lp_->applyRightPrec( *r_, *z_ );
621 (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
624 (*rTz_old_)(0,0) = rTz(0,0);
633 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
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.");
653 pnext_ = Teuchos::null;
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< 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.
RCGIterOrthoFailure(const std::string &what_arg)
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
RCGIterFailure is thrown when the RCGIter object is unable to compute the next iterate in the RCGIter...
RCGIterOrthoFailure is thrown when the RCGIter object is unable to compute independent direction vect...
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< 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.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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...
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.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
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
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
virtual ~RCGIter()
Destructor.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
RCGIterLAPACKFailure(const std::string &what_arg)
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...
Structure to contain pointers to RCGIter state variables.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
bool existU
Flag to indicate the recycle space should be used.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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...
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.
RCGIterFailure(const std::string &what_arg)
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
RCGIterInitFailure is thrown when the RCGIter object is unable to generate an initial iterate in the ...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
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.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
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...
Parent class to all Belos exceptions.
void resetNumIters(int iter=0)
Reset the iteration count.
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
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