42 #ifndef BELOS_BLOCK_CG_ITER_HPP 43 #define BELOS_BLOCK_CG_ITER_HPP 80 template<
class ScalarType,
class MV,
class OP,
81 const bool lapackSupportsScalarType =
151 void setStateSize() {
160 template<
class ScalarType,
class MV,
class OP>
322 bool stateStorageInitialized_;
344 template<
class ScalarType,
class MV,
class OP>
357 stateStorageInitialized_(false),
361 int bs = params.
get(
"Block Size", 1);
365 template <
class ScalarType,
class MV,
class OP>
368 if (! stateStorageInitialized_) {
372 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
373 stateStorageInitialized_ =
false;
384 (tmp == Teuchos::null,std:: invalid_argument,
385 "Belos::BlockCGIter::setStateSize: LinearProblem lacks " 386 "multivectors from which to clone.");
394 stateStorageInitialized_ =
true;
399 template <
class ScalarType,
class MV,
class OP>
405 (blockSize <= 0, std::invalid_argument,
"Belos::BlockGmresIter::" 406 "setBlockSize: blockSize = " << blockSize <<
" <= 0.");
407 if (blockSize == blockSize_) {
410 if (blockSize!=blockSize_) {
411 stateStorageInitialized_ =
false;
413 blockSize_ = blockSize;
414 initialized_ =
false;
419 template <
class ScalarType,
class MV,
class OP>
423 const char prefix[] =
"Belos::BlockCGIter::initialize: ";
426 if (! stateStorageInitialized_) {
431 (! stateStorageInitialized_, std::invalid_argument,
432 prefix <<
"Cannot initialize state storage!");
435 const char errstr[] =
"Specified multivectors must have a consistent " 442 if (newstate.
R != Teuchos::null) {
446 std::invalid_argument, prefix << errstr );
449 std::invalid_argument, prefix << errstr );
452 if (newstate.
R != R_) {
459 if ( lp_->getLeftPrec() != Teuchos::null ) {
460 lp_->applyLeftPrec( *R_, *Z_ );
461 if ( lp_->getRightPrec() != Teuchos::null ) {
463 lp_->applyRightPrec( *Z_, *tmp );
467 else if ( lp_->getRightPrec() != Teuchos::null ) {
468 lp_->applyRightPrec( *R_, *Z_ );
477 (newstate.
R == Teuchos::null, std::invalid_argument,
478 prefix <<
"BlockCGStateIterState does not have initial residual.");
485 template <
class ScalarType,
class MV,
class OP>
488 const char prefix[] =
"Belos::BlockCGIter::iterate: ";
493 if (initialized_ ==
false) {
507 rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
522 prefix <<
"Current linear system does not have the right number of vectors!" );
523 int rank = ortho_->normalize( *P_, Teuchos::null );
526 prefix <<
"Failed to compute initial block of orthonormal direction vectors.");
531 while (stest_->checkStatus(
this) !=
Passed) {
536 lp_->applyOp( *P_, *AP_ );
548 lltSolver.factorWithEquilibration(
true );
549 info = lltSolver.factor();
552 prefix <<
"Failed to compute Cholesky factorization using LAPACK routine POTRF.");
556 lltSolver.setVectors (Teuchos::rcpFromRef (alpha), Teuchos::rcpFromRef (alpha));
557 info = lltSolver.solve();
560 prefix <<
"Failed to compute alpha using Cholesky factorization (POTRS).");
564 lp_->updateSolution();
570 if ( lp_->getLeftPrec() != Teuchos::null ) {
571 lp_->applyLeftPrec( *R_, *Z_ );
572 if ( lp_->getRightPrec() != Teuchos::null ) {
574 lp_->applyRightPrec( *Z_, *tmp );
578 else if ( lp_->getRightPrec() != Teuchos::null ) {
579 lp_->applyRightPrec( *R_, *Z_ );
594 info = lltSolver.solve();
597 prefix <<
"Failed to compute beta using Cholesky factorization (POTRS).");
605 rank = ortho_->normalize( *P_, Teuchos::null );
608 prefix <<
"Failed to compute block of orthonormal direction vectors.");
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< ScalarType > SCT
int getNumIters() const
Get the current iteration count.
ScalarType * values() const
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
BlockCGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &, const Teuchos::RCP< OutputManager< ScalarType > > &, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &, Teuchos::ParameterList &)
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
T & get(const std::string &name, T def_value)
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MultiVecTraits< ScalarType, MV > MVT
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Declaration of basic traits for the multivector type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
int getBlockSize() const
Get the block size to be used by the iterative solver in solving this linear problem.
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.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
OperatorTraits< ScalarType, MV, OP > OPT
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void initializeCG(CGIterationState< ScalarType, MV > &)
Initialize the solver to an iterate, providing a complete state.
CGIterationLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void iterate()
This method performs linear solver iterations until the status test indicates the need to stop or an ...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
SCT::magnitudeType MagnitudeType
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
SCT::magnitudeType MagnitudeType
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
bool isInitialized()
States whether the solver has been initialized or not.
bool isInitialized()
States whether the solver has been initialized or not.
void resetNumIters(int iter=0)
Reset the iteration count to iter.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Teuchos::RCP< const MV > Z
The current preconditioned residual.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
OperatorTraits< ScalarType, MV, OP > OPT
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
CGIterationOrthoFailure is thrown when the CGIteration object is unable to compute independent direct...
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
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.
virtual ~BlockCGIter()
Destructor.
int getNumIters() const
Get the current iteration count.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > P
The current decent direction vector.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...