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;
195 template<
class ScalarType,
class MV,
class OP>
205 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
221 Teuchos::ParameterList ¶ms );
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_;
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();
561 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
RCGIterFailure,
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) {
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);
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();
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 );
633 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
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...
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.
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.
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.
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
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