42 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP 43 #define BELOS_PSEUDO_BLOCK_CG_ITER_HPP 60 #include "Teuchos_BLAS.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_ScalarTraits.hpp" 64 #include "Teuchos_ParameterList.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 80 template<
class ScalarType,
class MV,
class OP>
90 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
104 Teuchos::ParameterList ¶ms );
210 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
211 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
227 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
228 if (static_cast<size_type> (
iter_) >=
diag_.size ()) {
242 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
255 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
256 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
257 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
303 template<
class ScalarType,
class MV,
class OP>
307 Teuchos::ParameterList ¶ms ):
314 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
315 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) )
322 template <
class ScalarType,
class MV,
class OP>
326 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
327 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
328 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
329 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
332 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
335 int numRHS = MVT::GetNumberVecs(*tmp);
340 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
341 R_ = MVT::Clone( *tmp, numRHS_ );
342 Z_ = MVT::Clone( *tmp, numRHS_ );
343 P_ = MVT::Clone( *tmp, numRHS_ );
344 AP_ = MVT::Clone( *tmp, numRHS_ );
348 if(numEntriesForCondEst_ > 0) {
349 diag_.resize(numEntriesForCondEst_);
350 offdiag_.resize(numEntriesForCondEst_-1);
355 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
358 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
359 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
361 if (!Teuchos::is_null(newstate.
R)) {
363 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
364 std::invalid_argument, errstr );
365 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
366 std::invalid_argument, errstr );
369 if (newstate.
R != R_) {
371 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
377 if ( lp_->getLeftPrec() != Teuchos::null ) {
378 lp_->applyLeftPrec( *R_, *Z_ );
379 if ( lp_->getRightPrec() != Teuchos::null ) {
380 Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
381 lp_->applyRightPrec( *Z_, *tmp1 );
385 else if ( lp_->getRightPrec() != Teuchos::null ) {
386 lp_->applyRightPrec( *R_, *Z_ );
391 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
395 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
396 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
406 template <
class ScalarType,
class MV,
class OP>
412 if (initialized_ ==
false) {
418 std::vector<int> index(1);
419 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
420 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
423 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
424 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
427 ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
430 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
433 MVT::MvDot( *R_, *Z_, rHz );
435 if ( assertPositiveDefiniteness_ )
436 for (i=0; i<numRHS_; ++i)
437 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
439 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
444 while (stest_->checkStatus(
this) !=
Passed) {
450 lp_->applyOp( *P_, *AP_ );
453 MVT::MvDot( *P_, *AP_, pAp );
455 for (i=0; i<numRHS_; ++i) {
456 if ( assertPositiveDefiniteness_ )
458 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
460 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
462 alpha(i,i) = rHz[i] / pAp[i];
468 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
469 lp_->updateSolution();
473 for (i=0; i<numRHS_; ++i) {
479 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
484 if ( lp_->getLeftPrec() != Teuchos::null ) {
485 lp_->applyLeftPrec( *R_, *Z_ );
486 if ( lp_->getRightPrec() != Teuchos::null ) {
487 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
488 lp_->applyRightPrec( *Z_, *tmp );
492 else if ( lp_->getRightPrec() != Teuchos::null ) {
493 lp_->applyRightPrec( *R_, *Z_ );
499 MVT::MvDot( *R_, *Z_, rHz );
500 if ( assertPositiveDefiniteness_ )
501 for (i=0; i<numRHS_; ++i)
502 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
504 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
507 for (i=0; i<numRHS_; ++i) {
508 beta(i,i) = rHz[i] / rHz_old[i];
510 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
511 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
512 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
518 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
519 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
522 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
524 rHz_old2 = rHz_old[0];
525 beta_old = beta(0,0);
Teuchos::RCP< const MV > R
The current residual.
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.
int numEntriesForCondEst_
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.
virtual ~PseudoBlockCGIter()
Destructor.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::ArrayRCP< MagnitudeType > offdiag_
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize.
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.
OperatorTraits< ScalarType, MV, OP > OPT
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Traits class which defines basic operations on multivectors.
void resetNumIters(int iter=0)
Reset the iteration count.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
SCT::magnitudeType MagnitudeType
int getNumIters() const
Get the current iteration count.
MultiVecTraits< ScalarType, MV > MVT
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
bool assertPositiveDefiniteness_
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
PseudoBlockCGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList ¶ms)
PseudoBlockCGIter constructor with linear problem, solver utilities, and parameter list of solver opt...
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< const MV > P
The current decent direction vector.
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ArrayRCP< MagnitudeType > diag_