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 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
156 void setStateSize() {
157 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
165 template<
class ScalarType,
class MV,
class OP>
175 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
190 Teuchos::ParameterList ¶ms );
297 Teuchos::ArrayView<MagnitudeType> temp;
303 Teuchos::ArrayView<MagnitudeType> temp;
319 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
320 const Teuchos::RCP<OutputManager<ScalarType> > om_;
321 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
322 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
341 bool stateStorageInitialized_;
359 Teuchos::RCP<MV> AP_;
363 template<
class ScalarType,
class MV,
class OP>
369 Teuchos::ParameterList& params) :
376 stateStorageInitialized_(false),
380 int bs = params.get(
"Block Size", 1);
384 template <
class ScalarType,
class MV,
class OP>
387 if (! stateStorageInitialized_) {
389 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
390 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
391 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
392 stateStorageInitialized_ =
false;
401 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
402 TEUCHOS_TEST_FOR_EXCEPTION
403 (tmp == Teuchos::null,std:: invalid_argument,
404 "Belos::BlockCGIter::setStateSize: LinearProblem lacks " 405 "multivectors from which to clone.");
413 stateStorageInitialized_ =
true;
418 template <
class ScalarType,
class MV,
class OP>
423 TEUCHOS_TEST_FOR_EXCEPTION
424 (blockSize <= 0, std::invalid_argument,
"Belos::BlockGmresIter::" 425 "setBlockSize: blockSize = " << blockSize <<
" <= 0.");
426 if (blockSize == blockSize_) {
429 if (blockSize!=blockSize_) {
430 stateStorageInitialized_ =
false;
432 blockSize_ = blockSize;
433 initialized_ =
false;
438 template <
class ScalarType,
class MV,
class OP>
442 const char prefix[] =
"Belos::BlockCGIter::initialize: ";
445 if (! stateStorageInitialized_) {
449 TEUCHOS_TEST_FOR_EXCEPTION
450 (! stateStorageInitialized_, std::invalid_argument,
451 prefix <<
"Cannot initialize state storage!");
454 const char errstr[] =
"Specified multivectors must have a consistent " 460 if (newstate.
R != Teuchos::null) {
462 TEUCHOS_TEST_FOR_EXCEPTION
464 std::invalid_argument, prefix << errstr );
465 TEUCHOS_TEST_FOR_EXCEPTION
467 std::invalid_argument, prefix << errstr );
470 if (newstate.
R != R_) {
477 if ( lp_->getLeftPrec() != Teuchos::null ) {
478 lp_->applyLeftPrec( *R_, *Z_ );
479 if ( lp_->getRightPrec() != Teuchos::null ) {
480 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
481 lp_->applyRightPrec( *Z_, *tmp );
485 else if ( lp_->getRightPrec() != Teuchos::null ) {
486 lp_->applyRightPrec( *R_, *Z_ );
494 TEUCHOS_TEST_FOR_EXCEPTION
495 (newstate.
R == Teuchos::null, std::invalid_argument,
496 prefix <<
"BlockCGStateIterState does not have initial residual.");
503 template <
class ScalarType,
class MV,
class OP>
506 const char prefix[] =
"Belos::BlockCGIter::iterate: ";
511 if (initialized_ ==
false) {
519 Teuchos::LAPACK<int,ScalarType> lapack;
522 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
523 Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
524 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ),
525 rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
526 Teuchos::SerialSymDenseMatrix<int,ScalarType> pApHerm(Teuchos::View, uplo, pAp.values(), blockSize_, blockSize_);
529 Teuchos::SerialSpdDenseSolver<int,ScalarType> lltSolver;
532 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
535 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
538 TEUCHOS_TEST_FOR_EXCEPTION
540 prefix <<
"Current linear system does not have the right number of vectors!" );
541 int rank = ortho_->normalize( *P_, Teuchos::null );
542 TEUCHOS_TEST_FOR_EXCEPTION
544 prefix <<
"Failed to compute initial block of orthonormal direction vectors.");
549 while (stest_->checkStatus(
this) !=
Passed) {
554 lp_->applyOp( *P_, *AP_ );
565 lltSolver.setMatrix( Teuchos::rcp(&pApHerm,
false) );
566 lltSolver.factorWithEquilibration(
true );
567 info = lltSolver.factor();
568 TEUCHOS_TEST_FOR_EXCEPTION
570 prefix <<
"Failed to compute Cholesky factorization using LAPACK routine POTRF.");
574 lltSolver.setVectors (Teuchos::rcpFromRef (alpha), Teuchos::rcpFromRef (alpha));
575 info = lltSolver.solve();
576 TEUCHOS_TEST_FOR_EXCEPTION
578 prefix <<
"Failed to compute alpha using Cholesky factorization (POTRS).");
582 lp_->updateSolution();
588 if ( lp_->getLeftPrec() != Teuchos::null ) {
589 lp_->applyLeftPrec( *R_, *Z_ );
590 if ( lp_->getRightPrec() != Teuchos::null ) {
591 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
592 lp_->applyRightPrec( *Z_, *tmp );
596 else if ( lp_->getRightPrec() != Teuchos::null ) {
597 lp_->applyRightPrec( *R_, *Z_ );
611 lltSolver.setVectors( Teuchos::rcp( &beta,
false ), Teuchos::rcp( &beta,
false ) );
612 info = lltSolver.solve();
613 TEUCHOS_TEST_FOR_EXCEPTION
615 prefix <<
"Failed to compute beta using Cholesky factorization (POTRS).");
623 rank = ortho_->normalize( *P_, Teuchos::null );
624 TEUCHOS_TEST_FOR_EXCEPTION
626 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...
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
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 Assign(const MV &A, MV &mv)
mv := A
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 setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
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.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation (NOT_IMPLEMENTED)
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...