42 #ifndef BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP 43 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_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" 84 template<
class ScalarType,
class MV,
class OP>
94 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
108 Teuchos::ParameterList ¶ms );
218 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
219 "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
230 inline ScalarType normal() {
234 const double p0 = -0.322232431088;
235 const double p1 = -1.0;
236 const double p2 = -0.342242088547;
237 const double p3 = -0.204231210245e-1;
238 const double p4 = -0.453642210148e-4;
239 const double q0 = 0.993484626060e-1;
240 const double q1 = 0.588581570495;
241 const double q2 = 0.531103462366;
242 const double q3 = 0.103537752850;
243 const double q4 = 0.38560700634e-2;
247 r=0.5*(Teuchos::ScalarTraits<double>::random() + 1.0);
250 if(r < 0.5) y=std::sqrt(-2.0 * log(r));
251 else y=std::sqrt(-2.0 * log(1.0 - r));
253 p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
254 q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
256 if(r < 0.5) z = (p / q) - y;
257 else z = y - (p / q);
259 return Teuchos::as<ScalarType,double>(z);
265 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
266 const Teuchos::RCP<OutputManager<ScalarType> > om_;
267 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
287 bool assertPositiveDefiniteness_;
302 Teuchos::RCP<MV> AP_;
312 template<
class ScalarType,
class MV,
class OP>
316 Teuchos::ParameterList ¶ms ):
323 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) )
330 template <
class ScalarType,
class MV,
class OP>
334 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
335 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
336 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
337 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
340 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
343 int numRHS = MVT::GetNumberVecs(*tmp);
348 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
349 R_ = MVT::Clone( *tmp, numRHS_ );
350 Z_ = MVT::Clone( *tmp, numRHS_ );
351 P_ = MVT::Clone( *tmp, numRHS_ );
352 AP_ = MVT::Clone( *tmp, numRHS_ );
353 Y_ = MVT::Clone( *tmp, numRHS_ );
358 std::string errstr(
"Belos::BlockPseudoStochasticCGIter::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> tmp2 = MVT::Clone( *Z_, numRHS_ );
384 lp_->applyRightPrec( *Z_, *tmp2 );
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::StochasticCGIter::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_ );
423 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(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::PseudoBlockStochasticCGIter::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::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
462 alpha(i,i) = rHz[i] / pAp[i];
466 ScalarType z = normal();
467 zeta(i,i) = z / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
473 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
474 lp_->updateSolution();
477 MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
482 for (i=0; i<numRHS_; ++i) {
488 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
493 if ( lp_->getLeftPrec() != Teuchos::null ) {
494 lp_->applyLeftPrec( *R_, *Z_ );
495 if ( lp_->getRightPrec() != Teuchos::null ) {
496 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
497 lp_->applyRightPrec( *Z_, *tmp );
501 else if ( lp_->getRightPrec() != Teuchos::null ) {
502 lp_->applyRightPrec( *R_, *Z_ );
508 MVT::MvDot( *R_, *Z_, rHz );
509 if ( assertPositiveDefiniteness_ )
510 for (i=0; i<numRHS_; ++i)
511 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
513 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
516 for (i=0; i<numRHS_; ++i) {
517 beta(i,i) = rHz[i] / rHz_old[i];
519 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
520 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
521 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
Pure virtual base class which augments the basic interface for a stochastic conjugate gradient linear...
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...
Class which manages the output and verbosity of the Belos solvers.
MultiVecTraits< ScalarType, MV > MVT
SCT::magnitudeType MagnitudeType
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
void setBlockSize(int blockSize)
Set the blocksize.
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
virtual ~PseudoBlockStochasticCGIter()
Destructor.
Traits class which defines basic operations on multivectors.
void resetNumIters(int iter=0)
Reset the iteration count.
OperatorTraits< ScalarType, MV, OP > OPT
void iterate()
This method performs stochastic CG iterations on each linear system until the status test indicates t...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
A linear system to solve, and its associated information.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector.
bool isInitialized()
States whether the solver has been initialized or not.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< const MV > R
The current residual.
PseudoBlockStochasticCGIter(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)
PseudoBlockStochasticCGIter constructor with linear problem, solver utilities, and parameter list of ...
void initializeCG(StochasticCGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
This class implements the stochastic pseudo-block CG iteration, where the basic stochastic CG algorit...
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const MV > Y
The current stochastic recurrence vector.
Teuchos::ScalarTraits< ScalarType > SCT
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Structure to contain pointers to CGIteration state variables.