42 #ifndef BELOS_BLOCK_CG_ITER_HPP 43 #define BELOS_BLOCK_CG_ITER_HPP 60 #include "Teuchos_BLAS.hpp" 61 #include "Teuchos_LAPACK.hpp" 62 #include "Teuchos_SerialDenseMatrix.hpp" 63 #include "Teuchos_SerialDenseVector.hpp" 64 #include "Teuchos_SerialSymDenseMatrix.hpp" 65 #include "Teuchos_SerialSpdDenseSolver.hpp" 66 #include "Teuchos_ScalarTraits.hpp" 67 #include "Teuchos_ParameterList.hpp" 68 #include "Teuchos_TimeMonitor.hpp" 80 template<
class ScalarType,
class MV,
class OP,
81 const bool lapackSupportsScalarType =
87 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
94 Teuchos::ParameterList & )
96 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
102 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
106 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
110 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
114 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
118 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
122 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
125 Teuchos::RCP<const MV>
127 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
131 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
135 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
139 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
143 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
147 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
151 void setStateSize() {
152 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
160 template<
class ScalarType,
class MV,
class OP>
170 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
185 Teuchos::ParameterList ¶ms );
300 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
301 const Teuchos::RCP<OutputManager<ScalarType> > om_;
302 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
303 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
322 bool stateStorageInitialized_;
340 Teuchos::RCP<MV> AP_;
344 template<
class ScalarType,
class MV,
class OP>
350 Teuchos::ParameterList& params) :
357 stateStorageInitialized_(false),
361 int bs = params.get(
"Block Size", 1);
365 template <
class ScalarType,
class MV,
class OP>
368 if (! stateStorageInitialized_) {
370 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
371 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
372 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
373 stateStorageInitialized_ =
false;
382 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
383 TEUCHOS_TEST_FOR_EXCEPTION
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>
404 TEUCHOS_TEST_FOR_EXCEPTION
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_) {
430 TEUCHOS_TEST_FOR_EXCEPTION
431 (! stateStorageInitialized_, std::invalid_argument,
432 prefix <<
"Cannot initialize state storage!");
435 const char errstr[] =
"Specified multivectors must have a consistent " 439 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
442 if (newstate.
R != Teuchos::null) {
444 TEUCHOS_TEST_FOR_EXCEPTION
446 std::invalid_argument, prefix << errstr );
447 TEUCHOS_TEST_FOR_EXCEPTION
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 ) {
462 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
463 lp_->applyRightPrec( *Z_, *tmp );
467 else if ( lp_->getRightPrec() != Teuchos::null ) {
468 lp_->applyRightPrec( *R_, *Z_ );
476 TEUCHOS_TEST_FOR_EXCEPTION
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) {
501 Teuchos::LAPACK<int,ScalarType> lapack;
504 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
505 Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
506 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ),
507 rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
508 Teuchos::SerialSymDenseMatrix<int,ScalarType> pApHerm(Teuchos::View, uplo, pAp.values(), blockSize_, blockSize_);
511 Teuchos::SerialSpdDenseSolver<int,ScalarType> lltSolver;
514 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
517 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
520 TEUCHOS_TEST_FOR_EXCEPTION
522 prefix <<
"Current linear system does not have the right number of vectors!" );
523 int rank = ortho_->normalize( *P_, Teuchos::null );
524 TEUCHOS_TEST_FOR_EXCEPTION
526 prefix <<
"Failed to compute initial block of orthonormal direction vectors.");
531 while (stest_->checkStatus(
this) !=
Passed) {
536 lp_->applyOp( *P_, *AP_ );
547 lltSolver.setMatrix( Teuchos::rcp(&pApHerm,
false) );
548 lltSolver.factorWithEquilibration(
true );
549 info = lltSolver.factor();
550 TEUCHOS_TEST_FOR_EXCEPTION
552 prefix <<
"Failed to compute Cholesky factorization using LAPACK routine POTRF.");
556 lltSolver.setVectors (Teuchos::rcpFromRef (alpha), Teuchos::rcpFromRef (alpha));
557 info = lltSolver.solve();
558 TEUCHOS_TEST_FOR_EXCEPTION
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 ) {
573 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
574 lp_->applyRightPrec( *Z_, *tmp );
578 else if ( lp_->getRightPrec() != Teuchos::null ) {
579 lp_->applyRightPrec( *R_, *Z_ );
593 lltSolver.setVectors( Teuchos::rcp( &beta,
false ), Teuchos::rcp( &beta,
false ) );
594 info = lltSolver.solve();
595 TEUCHOS_TEST_FOR_EXCEPTION
597 prefix <<
"Failed to compute beta using Cholesky factorization (POTRS).");
605 rank = ortho_->normalize( *P_, Teuchos::null );
606 TEUCHOS_TEST_FOR_EXCEPTION
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.
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...
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class for defining the status testing capabilities of Belos.
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.
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...