10 #ifndef BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
11 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
52 template<
class ScalarType,
class MV,
class OP>
187 "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
202 const double p0 = -0.322232431088;
203 const double p1 = -1.0;
204 const double p2 = -0.342242088547;
205 const double p3 = -0.204231210245e-1;
206 const double p4 = -0.453642210148e-4;
207 const double q0 = 0.993484626060e-1;
208 const double q1 = 0.588581570495;
209 const double q2 = 0.531103462366;
210 const double q3 = 0.103537752850;
211 const double q4 = 0.38560700634e-2;
215 Teuchos::randomSyncedMatrix( randvec_ );
217 for (
int i=0; i<numRHS_; i++)
223 if(r < 0.5) y=std::sqrt(-2.0 * log(r));
224 else y=std::sqrt(-2.0 * log(1.0 - r));
226 p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
227 q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
229 if(r < 0.5) z = (p / q) - y;
230 else z = y - (p / q);
232 randvec_[i] = Teuchos::as<ScalarType,double>(z);
263 bool assertPositiveDefiniteness_;
290 template<
class ScalarType,
class MV,
class OP>
301 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) )
308 template <
class ScalarType,
class MV,
class OP>
315 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
321 int numRHS = MVT::GetNumberVecs(*tmp);
327 R_ = MVT::Clone( *tmp, numRHS_ );
328 Z_ = MVT::Clone( *tmp, numRHS_ );
329 P_ = MVT::Clone( *tmp, numRHS_ );
330 AP_ = MVT::Clone( *tmp, numRHS_ );
331 Y_ = MVT::Clone( *tmp, numRHS_ );
335 randvec_.size( numRHS_ );
339 std::string errstr(
"Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
348 std::invalid_argument, errstr );
350 std::invalid_argument, errstr );
353 if (newstate.
R != R_) {
355 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
361 if ( lp_->getLeftPrec() != Teuchos::null ) {
362 lp_->applyLeftPrec( *R_, *Z_ );
363 if ( lp_->getRightPrec() != Teuchos::null ) {
365 lp_->applyRightPrec( *Z_, *tmp2 );
369 else if ( lp_->getRightPrec() != Teuchos::null ) {
370 lp_->applyRightPrec( *R_, *Z_ );
375 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
380 "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
390 template <
class ScalarType,
class MV,
class OP>
396 if (initialized_ ==
false) {
402 std::vector<int> index(1);
403 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
414 MVT::MvDot( *R_, *Z_, rHz );
416 if ( assertPositiveDefiniteness_ )
417 for (i=0; i<numRHS_; ++i)
420 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
425 while (stest_->checkStatus(
this) !=
Passed) {
431 lp_->applyOp( *P_, *AP_ );
434 MVT::MvDot( *P_, *AP_, pAp );
438 for (i=0; i<numRHS_; ++i) {
439 if ( assertPositiveDefiniteness_ )
443 "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
445 alpha(i,i) = rHz[i] / pAp[i];
454 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
455 lp_->updateSolution();
458 MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
463 for (i=0; i<numRHS_; ++i) {
469 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
474 if ( lp_->getLeftPrec() != Teuchos::null ) {
475 lp_->applyLeftPrec( *R_, *Z_ );
476 if ( lp_->getRightPrec() != Teuchos::null ) {
478 lp_->applyRightPrec( *Z_, *tmp );
482 else if ( lp_->getRightPrec() != Teuchos::null ) {
483 lp_->applyRightPrec( *R_, *Z_ );
489 MVT::MvDot( *R_, *Z_, rHz );
490 if ( assertPositiveDefiniteness_ )
491 for (i=0; i<numRHS_; ++i)
494 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
497 for (i=0; i<numRHS_; ++i) {
498 beta(i,i) = rHz[i] / rHz_old[i];
502 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
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.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(ScalarTypea)
Declaration of basic traits for the multivector type.
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.
virtual ~PseudoBlockStochasticCGIter()
Destructor.
int getNumIters() const
Get the current iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Traits class which defines basic operations on multivectors.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero...
void resetNumIters(int iter=0)
Reset the iteration count.
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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...
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.