42 #ifndef BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP 43 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP 84 template<
class ScalarType,
class MV,
class OP>
219 "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
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;
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);
312 template<
class ScalarType,
class MV,
class OP>
323 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) )
330 template <
class ScalarType,
class MV,
class OP>
337 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
343 int numRHS = MVT::GetNumberVecs(*tmp);
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.");
367 std::invalid_argument, errstr );
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 ) {
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_ );
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_ );
433 MVT::MvDot( *R_, *Z_, rHz );
435 if ( assertPositiveDefiniteness_ )
436 for (i=0; i<numRHS_; ++i)
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_ )
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();
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 ) {
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)
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];
521 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
const Teuchos::RCP< OutputManager< ScalarType > > om_
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
void setBlockSize(int blockSize)
Set the blocksize.
bool is_null(const std::shared_ptr< T > &p)
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
ScalarType normal()
Wrapper for Normal(0,1) random variables.
bool assertPositiveDefiniteness_
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.