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"
79template<
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");
155 void setStateSize() {
156 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Stub");
164template<
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_;
340 bool stateStorageInitialized_;
358 Teuchos::RCP<MV> AP_;
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;
398 if (R_ == Teuchos::null || MVT::GetNumberVecs(*R_)!=blockSize_) {
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.");
405 R_ = MVT::Clone (*tmp, blockSize_);
406 Z_ = MVT::Clone (*tmp, blockSize_);
407 P_ = MVT::Clone (*tmp, blockSize_);
408 AP_ = MVT::Clone (*tmp, blockSize_);
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
462 (MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
463 std::invalid_argument, prefix << errstr );
464 TEUCHOS_TEST_FOR_EXCEPTION
465 (MVT::GetNumberVecs(*newstate.
R) != blockSize_,
466 std::invalid_argument, prefix << errstr );
469 if (newstate.
R != R_) {
471 MVT::Assign( *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_ );
490 MVT::Assign( *Z_, *P_ );
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_ );
560 MVT::MvTransMv( one, *P_, *R_, alpha );
561 MVT::MvTransMv( one, *P_, *AP_, pAp );
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).");
580 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
581 lp_->updateSolution();
584 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
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_ );
608 MVT::MvTransMv( -one, *AP_, *Z_, beta );
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).");
617 Teuchos::RCP<MV> Pnew = MVT::CloneCopy( *Z_ );
618 MVT::MvTimesMatAddMv(one, *P_, beta, one, *Pnew);
622 rank = ortho_->normalize( *P_, Teuchos::null );
623 TEUCHOS_TEST_FOR_EXCEPTION
625 prefix <<
"Failed to compute block of orthonormal direction vectors.");
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.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::ScalarTraits< ScalarType > SCT
void setDoCondEst(bool)
Sets whether or not to store the diagonal for condition estimation.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
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.
MultiVecTraits< ScalarType, MV > MVT
void resetNumIters(int iter=0)
Reset the iteration count.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
virtual ~BlockCGIter()
Destructor.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
int getBlockSize() const
Get the block size to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
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 &)
SCT::magnitudeType MagnitudeType
void resetNumIters(int iter=0)
Reset the iteration count to iter.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
void iterate()
This method performs linear solver iterations until the status test indicates the need to stop or an ...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void initializeCG(CGIterationState< ScalarType, MV > &)
Initialize the solver to an iterate, providing a complete state.
Teuchos::ScalarTraits< ScalarType > SCT
bool isInitialized()
States whether the solver has been initialized or not.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
MultiVecTraits< ScalarType, MV > MVT
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.
CGIterationOrthoFailure is thrown when the CGIteration object is unable to compute independent direct...
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.