42 #ifndef BELOS_RCG_ITER_HPP 43 #define BELOS_RCG_ITER_HPP 59 #include "Teuchos_LAPACK.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" 90 template <
class ScalarType,
class MV>
118 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Alpha;
119 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Beta;
120 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
D;
121 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
rTz_old;
124 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Delta;
127 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
LUUTAU;
129 Teuchos::RCP<std::vector<int> >
ipiv;
135 U(Teuchos::null),
AU(Teuchos::null),
136 Alpha(Teuchos::null),
Beta(Teuchos::null),
D(Teuchos::null),
rTz_old(Teuchos::null),
184 template<
class ScalarType,
class MV,
class OP>
194 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
210 Teuchos::ParameterList ¶ms );
284 if (!initialized_)
return 0;
317 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
318 "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
322 void setSize(
int recycleBlocks,
int numBlocks );
338 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
339 const Teuchos::RCP<OutputManager<ScalarType> > om_;
340 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
369 Teuchos::RCP<MV> Ap_;
380 Teuchos::RCP<MV> U_, AU_;
383 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_,Beta_,D_;
386 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
389 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
392 Teuchos::RCP<std::vector<int> > ipiv_;
395 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
400 template<
class ScalarType,
class MV,
class OP>
404 Teuchos::ParameterList ¶ms ):
416 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
417 "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
418 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
420 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
421 "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
422 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
430 template <
class ScalarType,
class MV,
class OP>
434 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
435 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument,
"Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
436 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument,
"Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
438 numBlocks_ = numBlocks;
439 recycleBlocks_ = recycleBlocks;
445 template <
class ScalarType,
class MV,
class OP>
449 if (newstate.
P != Teuchos::null &&
450 newstate.
Ap != Teuchos::null &&
451 newstate.
r != Teuchos::null &&
452 newstate.
z != Teuchos::null &&
453 newstate.
U != Teuchos::null &&
454 newstate.
AU != Teuchos::null &&
455 newstate.
Alpha != Teuchos::null &&
456 newstate.
Beta != Teuchos::null &&
457 newstate.
D != Teuchos::null &&
458 newstate.
Delta != Teuchos::null &&
459 newstate.
LUUTAU != Teuchos::null &&
460 newstate.
ipiv != Teuchos::null &&
461 newstate.
rTz_old != Teuchos::null) {
463 curDim_ = newstate.
curDim;
468 existU_ = newstate.
existU;
471 Alpha_ = newstate.
Alpha;
472 Beta_ = newstate.
Beta;
474 Delta_ = newstate.
Delta;
475 LUUTAU_ = newstate.
LUUTAU;
476 ipiv_ = newstate.
ipiv;
481 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
P == Teuchos::null,std::invalid_argument,
482 "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
484 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Ap == Teuchos::null,std::invalid_argument,
485 "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
487 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
r == Teuchos::null,std::invalid_argument,
488 "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
490 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
z == Teuchos::null,std::invalid_argument,
491 "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
493 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
U == Teuchos::null,std::invalid_argument,
494 "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
496 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
AU == Teuchos::null,std::invalid_argument,
497 "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
499 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Alpha == Teuchos::null,std::invalid_argument,
500 "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
502 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Beta == Teuchos::null,std::invalid_argument,
503 "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
505 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
D == Teuchos::null,std::invalid_argument,
506 "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
508 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Delta == Teuchos::null,std::invalid_argument,
509 "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
511 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
LUUTAU == Teuchos::null,std::invalid_argument,
512 "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
514 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
ipiv == Teuchos::null,std::invalid_argument,
515 "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
517 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
rTz_old == Teuchos::null,std::invalid_argument,
518 "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
529 template <
class ScalarType,
class MV,
class OP>
532 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
RCGIterFailure,
533 "Belos::RCGIter::iterate(): RCGIter class not initialized." );
536 Teuchos::LAPACK<int,ScalarType> lapack;
539 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
540 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
543 std::vector<int> index(1);
544 Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
547 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
550 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
RCGIterFailure,
551 "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
554 int searchDim = numBlocks_+1;
564 Teuchos::RCP<const MV> p_ = Teuchos::null;
565 Teuchos::RCP<MV> pnext_ = Teuchos::null;
566 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= searchDim) {
571 p_ = MVT::CloneView( *P_, index );
572 lp_->applyOp( *p_, *Ap_ );
575 MVT::MvTransMv( one, *p_, *Ap_, pAp );
576 (*D_)(i_,0) = pAp(0,0);
579 (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
582 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero,
RCGIterFailure,
"Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
585 MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
586 lp_->updateSolution();
589 MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
591 std::vector<MagnitudeType> norm(1);
592 MVT::MvNorm( *r_, norm );
596 if ( lp_->getLeftPrec() != Teuchos::null ) {
597 lp_->applyLeftPrec( *r_, *z_ );
599 else if ( lp_->getRightPrec() != Teuchos::null ) {
600 lp_->applyRightPrec( *r_, *z_ );
607 MVT::MvTransMv( one, *r_, *z_, rTz );
610 (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
613 (*rTz_old_)(0,0) = rTz(0,0);
618 pnext_ = MVT::CloneViewNonConst( *P_, index );
622 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
623 MVT::MvTransMv( one, *AU_, *z_, mu );
626 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
628 "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
631 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
633 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
637 MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
642 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.
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.
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.
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.
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...
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...
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 ...
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...
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
SCT::magnitudeType MagnitudeType