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.");
229 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
230 if (static_cast<size_type> (
iter_) >=
diag_.size ()) {
244 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
257 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
258 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
259 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
306 template<
class ScalarType,
class MV,
class OP>
310 Teuchos::ParameterList ¶ms ):
317 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
318 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
326 template <
class ScalarType,
class MV,
class OP>
330 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
331 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
332 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
333 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
336 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
339 int numRHS = MVT::GetNumberVecs(*tmp);
344 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
345 R_ = MVT::Clone( *tmp, numRHS_ );
346 Z_ = MVT::Clone( *tmp, numRHS_ );
347 P_ = MVT::Clone( *tmp, numRHS_ );
348 AP_ = MVT::Clone( *tmp, numRHS_ );
352 if(numEntriesForCondEst_ > 0) {
353 diag_.resize(numEntriesForCondEst_);
354 offdiag_.resize(numEntriesForCondEst_-1);
359 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
362 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
363 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
365 if (!Teuchos::is_null(newstate.
R)) {
367 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
368 std::invalid_argument, errstr );
369 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
370 std::invalid_argument, errstr );
373 if (newstate.
R != R_) {
375 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
381 if ( lp_->getLeftPrec() != Teuchos::null ) {
382 lp_->applyLeftPrec( *R_, *Z_ );
383 if ( lp_->getRightPrec() != Teuchos::null ) {
384 Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
385 lp_->applyRightPrec( *Z_, *tmp1 );
389 else if ( lp_->getRightPrec() != Teuchos::null ) {
390 lp_->applyRightPrec( *R_, *Z_ );
395 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
399 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
400 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
410 template <
class ScalarType,
class MV,
class OP>
416 if (initialized_ ==
false) {
422 std::vector<int> index(1);
423 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ ), beta( numRHS_ );
424 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ );
427 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
428 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
431 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
434 MVT::MvDot( *R_, *Z_, rHz );
436 if ( assertPositiveDefiniteness_ )
437 for (i=0; i<numRHS_; ++i)
438 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
440 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
445 while (stest_->checkStatus(
this) !=
Passed) {
451 lp_->applyOp( *P_, *AP_ );
454 MVT::MvDot( *P_, *AP_, pAp );
456 for (i=0; i<numRHS_; ++i) {
457 if ( assertPositiveDefiniteness_ )
459 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
461 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
463 alpha(i,i) = rHz[i] / pAp[i];
469 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
470 lp_->updateSolution();
474 for (i=0; i<numRHS_; ++i) {
480 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
485 if ( lp_->getLeftPrec() != Teuchos::null ) {
486 lp_->applyLeftPrec( *R_, *Z_ );
487 if ( lp_->getRightPrec() != Teuchos::null ) {
488 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
489 lp_->applyRightPrec( *Z_, *tmp );
493 else if ( lp_->getRightPrec() != Teuchos::null ) {
494 lp_->applyRightPrec( *R_, *Z_ );
500 MVT::MvDot( *R_, *Z_, rHz );
501 if ( assertPositiveDefiniteness_ )
502 for (i=0; i<numRHS_; ++i)
503 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
505 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
508 for (i=0; i<numRHS_; ++i) {
509 beta[i] = rHz[i] / rHz_old[i];
511 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
512 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
513 MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
519 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old_ * beta_old_ * pAp_old_ + pAp[0]) / rHz_old[0]);
520 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old_ * pAp_old_ / (sqrt( rHz_old[0] * rHz_old2_)));
523 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
525 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.
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
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...
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_