42#ifndef BELOS_RCG_ITER_HPP
43#define BELOS_RCG_ITER_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;
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>
154 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
170 Teuchos::ParameterList ¶ms );
277 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
278 "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
282 void setSize(
int recycleBlocks,
int numBlocks );
298 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
299 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
300 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
346 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Delta_;
349 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
LUUTAU_;
352 Teuchos::RCP<std::vector<int> >
ipiv_;
355 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
rTz_old_;
360 template<
class ScalarType,
class MV,
class OP>
364 Teuchos::ParameterList ¶ms ):
376 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
377 "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
378 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
380 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
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>
409 if (newstate.
P != Teuchos::null &&
410 newstate.
Ap != Teuchos::null &&
411 newstate.
r != Teuchos::null &&
412 newstate.
z != Teuchos::null &&
413 newstate.
U != Teuchos::null &&
414 newstate.
AU != Teuchos::null &&
415 newstate.
Alpha != Teuchos::null &&
416 newstate.
Beta != Teuchos::null &&
417 newstate.
D != Teuchos::null &&
418 newstate.
Delta != Teuchos::null &&
419 newstate.
LUUTAU != Teuchos::null &&
420 newstate.
ipiv != Teuchos::null &&
421 newstate.
rTz_old != Teuchos::null) {
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;
441 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
P == Teuchos::null,std::invalid_argument,
442 "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
444 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Ap == Teuchos::null,std::invalid_argument,
445 "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
447 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
r == Teuchos::null,std::invalid_argument,
448 "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
450 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
z == Teuchos::null,std::invalid_argument,
451 "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
453 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
U == Teuchos::null,std::invalid_argument,
454 "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
456 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
AU == Teuchos::null,std::invalid_argument,
457 "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
459 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Alpha == Teuchos::null,std::invalid_argument,
460 "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
462 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Beta == Teuchos::null,std::invalid_argument,
463 "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
465 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
D == Teuchos::null,std::invalid_argument,
466 "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
468 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Delta == Teuchos::null,std::invalid_argument,
469 "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
471 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
LUUTAU == Teuchos::null,std::invalid_argument,
472 "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
474 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
ipiv == Teuchos::null,std::invalid_argument,
475 "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
477 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
rTz_old == Teuchos::null,std::invalid_argument,
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." );
496 Teuchos::LAPACK<int,ScalarType> lapack;
499 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
500 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
503 std::vector<int> index(1);
504 Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
507 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
510 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
511 "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
514 int searchDim = numBlocks_+1;
524 Teuchos::RCP<const MV> p_ = Teuchos::null;
525 Teuchos::RCP<MV> pnext_ = Teuchos::null;
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 );
557 if ( lp_->getLeftPrec() != Teuchos::null ) {
558 lp_->applyLeftPrec( *r_, *z_ );
560 else if ( lp_->getRightPrec() != Teuchos::null ) {
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 );
583 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
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_ );
603 pnext_ = Teuchos::null;
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
CGIterationLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Teuchos::ScalarTraits< ScalarType > SCT
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
MultiVecTraits< ScalarType, MV > MVT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
virtual ~RCGIter()
Destructor.
void resetNumIters(int iter=0)
Reset the iteration count.
void setSize(int recycleBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors.
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
void setBlockSize(int blockSize)
Set the blocksize.
void setRecycledBlocks(int recycleBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
SCT::magnitudeType MagnitudeType
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
int getRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
Teuchos::RCP< std::vector< int > > ipiv_
void iterate()
This method performs RCG iterations until the status test indicates the need to stop or an error occu...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to RCGIter state variables.
bool existU
Flag to indicate the recycle space should be used.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< MV > Ap
A times current search vector.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > r
The current residual.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
Teuchos::RCP< MV > P
The current Krylov basis.