42 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP 43 #define BELOS_PSEUDO_BLOCK_CG_ITER_HPP 60 #include "Teuchos_SerialDenseMatrix.hpp" 61 #include "Teuchos_SerialDenseVector.hpp" 62 #include "Teuchos_ScalarTraits.hpp" 63 #include "Teuchos_ParameterList.hpp" 64 #include "Teuchos_TimeMonitor.hpp" 79 template<
class ScalarType,
class MV,
class OP>
89 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
103 Teuchos::ParameterList ¶ms );
209 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
210 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
220 if (numEntriesForCondEst_) doCondEst_=val;
228 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
229 if (static_cast<size_type> (iter_) >= diag_.size ()) {
232 return diag_ (0, iter_);
243 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
244 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
247 return offdiag_ (0, iter_);
256 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
257 const Teuchos::RCP<OutputManager<ScalarType> > om_;
258 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
278 bool assertPositiveDefiniteness_;
281 Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
282 ScalarType pAp_old_, beta_old_, rHz_old2_;
283 int numEntriesForCondEst_;
299 Teuchos::RCP<MV> AP_;
305 template<
class ScalarType,
class MV,
class OP>
309 Teuchos::ParameterList ¶ms ):
316 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
317 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
325 template <
class ScalarType,
class MV,
class OP>
329 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
330 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
331 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
332 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
335 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
338 int numRHS = MVT::GetNumberVecs(*tmp);
343 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
344 R_ = MVT::Clone( *tmp, numRHS_ );
345 Z_ = MVT::Clone( *tmp, numRHS_ );
346 P_ = MVT::Clone( *tmp, numRHS_ );
347 AP_ = MVT::Clone( *tmp, numRHS_ );
351 if(numEntriesForCondEst_ > 0) {
352 diag_.resize(numEntriesForCondEst_);
353 offdiag_.resize(numEntriesForCondEst_-1);
358 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
361 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
362 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
364 if (!Teuchos::is_null(newstate.
R)) {
366 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
367 std::invalid_argument, errstr );
368 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
369 std::invalid_argument, errstr );
372 if (newstate.
R != R_) {
374 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
380 if ( lp_->getLeftPrec() != Teuchos::null ) {
381 lp_->applyLeftPrec( *R_, *Z_ );
382 if ( lp_->getRightPrec() != Teuchos::null ) {
383 Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
384 lp_->applyRightPrec( *Z_, *tmp1 );
388 else if ( lp_->getRightPrec() != Teuchos::null ) {
389 lp_->applyRightPrec( *R_, *Z_ );
394 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
398 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
399 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
409 template <
class ScalarType,
class MV,
class OP>
415 if (initialized_ ==
false) {
421 std::vector<int> index(1);
422 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ ), beta( numRHS_ );
423 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ );
426 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
427 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
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] = 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], *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];
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.
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.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
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.
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 ...
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
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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...
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 ...