42 #ifndef BELOS_BLOCK_CG_ITER_HPP 43 #define BELOS_BLOCK_CG_ITER_HPP 60 #include "Teuchos_LAPACK.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_SerialSymDenseMatrix.hpp" 64 #include "Teuchos_SerialSpdDenseSolver.hpp" 65 #include "Teuchos_ScalarTraits.hpp" 66 #include "Teuchos_ParameterList.hpp" 67 #include "Teuchos_TimeMonitor.hpp" 79 template<
class ScalarType,
class MV,
class OP,
80 const bool lapackSupportsScalarType =
86 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
93 Teuchos::ParameterList & )
95 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
101 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
105 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
109 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
113 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
117 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
121 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
124 Teuchos::RCP<const MV>
126 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
130 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
134 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
138 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
142 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
146 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
150 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
156 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
164 template<
class ScalarType,
class MV,
class OP>
174 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
189 Teuchos::ParameterList ¶ms );
296 Teuchos::ArrayView<MagnitudeType> temp;
302 Teuchos::ArrayView<MagnitudeType> temp;
318 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
319 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
320 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
321 const Teuchos::RCP<OrthoManager<ScalarType,MV> >
ortho_;
362 template<
class ScalarType,
class MV,
class OP>
368 Teuchos::ParameterList& params) :
375 stateStorageInitialized_(false),
379 int bs = params.get(
"Block Size", 1);
383 template <
class ScalarType,
class MV,
class OP>
386 if (! stateStorageInitialized_) {
388 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
389 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
390 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
391 stateStorageInitialized_ =
false;
400 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
401 TEUCHOS_TEST_FOR_EXCEPTION
402 (tmp == Teuchos::null,std:: invalid_argument,
403 "Belos::BlockCGIter::setStateSize: LinearProblem lacks " 404 "multivectors from which to clone.");
412 stateStorageInitialized_ =
true;
417 template <
class ScalarType,
class MV,
class OP>
422 TEUCHOS_TEST_FOR_EXCEPTION
423 (blockSize <= 0, std::invalid_argument,
"Belos::BlockGmresIter::" 424 "setBlockSize: blockSize = " << blockSize <<
" <= 0.");
425 if (blockSize == blockSize_) {
428 if (blockSize!=blockSize_) {
429 stateStorageInitialized_ =
false;
431 blockSize_ = blockSize;
432 initialized_ =
false;
437 template <
class ScalarType,
class MV,
class OP>
441 const char prefix[] =
"Belos::BlockCGIter::initialize: ";
444 if (! stateStorageInitialized_) {
448 TEUCHOS_TEST_FOR_EXCEPTION
449 (! stateStorageInitialized_, std::invalid_argument,
450 prefix <<
"Cannot initialize state storage!");
453 const char errstr[] =
"Specified multivectors must have a consistent " 459 if (newstate.
R != Teuchos::null) {
461 TEUCHOS_TEST_FOR_EXCEPTION
463 std::invalid_argument, prefix << errstr );
464 TEUCHOS_TEST_FOR_EXCEPTION
466 std::invalid_argument, prefix << errstr );
469 if (newstate.
R != R_) {
476 if ( lp_->getLeftPrec() != Teuchos::null ) {
477 lp_->applyLeftPrec( *R_, *Z_ );
478 if ( lp_->getRightPrec() != Teuchos::null ) {
479 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
480 lp_->applyRightPrec( *Z_, *tmp );
484 else if ( lp_->getRightPrec() != Teuchos::null ) {
485 lp_->applyRightPrec( *R_, *Z_ );
493 TEUCHOS_TEST_FOR_EXCEPTION
494 (newstate.
R == Teuchos::null, std::invalid_argument,
495 prefix <<
"BlockCGStateIterState does not have initial residual.");
502 template <
class ScalarType,
class MV,
class OP>
505 const char prefix[] =
"Belos::BlockCGIter::iterate: ";
510 if (initialized_ ==
false) {
518 Teuchos::LAPACK<int,ScalarType> lapack;
521 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
522 Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
523 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ),
524 rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
525 Teuchos::SerialSymDenseMatrix<int,ScalarType> pApHerm(Teuchos::View, uplo, pAp.values(), blockSize_, blockSize_);
528 Teuchos::SerialSpdDenseSolver<int,ScalarType> lltSolver;
531 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
534 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
537 TEUCHOS_TEST_FOR_EXCEPTION
539 prefix <<
"Current linear system does not have the right number of vectors!" );
540 int rank = ortho_->normalize( *P_, Teuchos::null );
541 TEUCHOS_TEST_FOR_EXCEPTION
543 prefix <<
"Failed to compute initial block of orthonormal direction vectors.");
548 while (stest_->checkStatus(
this) !=
Passed) {
553 lp_->applyOp( *P_, *AP_ );
564 lltSolver.setMatrix( Teuchos::rcp(&pApHerm,
false) );
565 lltSolver.factorWithEquilibration(
true );
566 info = lltSolver.factor();
567 TEUCHOS_TEST_FOR_EXCEPTION
569 prefix <<
"Failed to compute Cholesky factorization using LAPACK routine POTRF.");
573 lltSolver.setVectors (Teuchos::rcpFromRef (alpha), Teuchos::rcpFromRef (alpha));
574 info = lltSolver.solve();
575 TEUCHOS_TEST_FOR_EXCEPTION
577 prefix <<
"Failed to compute alpha using Cholesky factorization (POTRS).");
581 lp_->updateSolution();
587 if ( lp_->getLeftPrec() != Teuchos::null ) {
588 lp_->applyLeftPrec( *R_, *Z_ );
589 if ( lp_->getRightPrec() != Teuchos::null ) {
590 Teuchos::RCP<MV> tmp =
MVT::Clone( *Z_, blockSize_ );
591 lp_->applyRightPrec( *Z_, *tmp );
595 else if ( lp_->getRightPrec() != Teuchos::null ) {
596 lp_->applyRightPrec( *R_, *Z_ );
610 lltSolver.setVectors( Teuchos::rcp( &beta,
false ), Teuchos::rcp( &beta,
false ) );
611 info = lltSolver.solve();
612 TEUCHOS_TEST_FOR_EXCEPTION
614 prefix <<
"Failed to compute beta using Cholesky factorization (POTRS).");
622 rank = ortho_->normalize( *P_, Teuchos::null );
623 TEUCHOS_TEST_FOR_EXCEPTION
625 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.
const Teuchos::RCP< OutputManager< ScalarType > > om_
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.
void setDoCondEst(bool)
Sets whether or not to store the diagonal for condition estimation.
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
bool stateStorageInitialized_
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
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 .
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
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.
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 ...
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...