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;
247 Teuchos::randomSyncedMatrix(
randvec_ );
255 if(r < 0.5) y=std::sqrt(-2.0 * log(r));
256 else y=std::sqrt(-2.0 * log(1.0 - r));
258 p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
259 q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
261 if(r < 0.5) z = (p / q) - y;
262 else z = y - (p / q);
264 randvec_[i] = Teuchos::as<ScalarType,double>(z);
322 template<
class ScalarType,
class MV,
class OP>
333 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) )
340 template <
class ScalarType,
class MV,
class OP>
347 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
353 int numRHS = MVT::GetNumberVecs(*tmp);
359 R_ = MVT::Clone( *tmp, numRHS_ );
360 Z_ = MVT::Clone( *tmp, numRHS_ );
361 P_ = MVT::Clone( *tmp, numRHS_ );
362 AP_ = MVT::Clone( *tmp, numRHS_ );
363 Y_ = MVT::Clone( *tmp, numRHS_ );
367 randvec_.size( numRHS_ );
371 std::string errstr(
"Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
380 std::invalid_argument, errstr );
382 std::invalid_argument, errstr );
385 if (newstate.
R != R_) {
387 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
394 lp_->applyLeftPrec( *R_, *Z_ );
397 lp_->applyRightPrec( *Z_, *tmp2 );
402 lp_->applyRightPrec( *R_, *Z_ );
407 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
412 "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
422 template <
class ScalarType,
class MV,
class OP>
428 if (initialized_ ==
false) {
434 std::vector<int> index(1);
435 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
446 MVT::MvDot( *R_, *Z_, rHz );
448 if ( assertPositiveDefiniteness_ )
449 for (i=0; i<numRHS_; ++i)
452 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
457 while (stest_->checkStatus(
this) !=
Passed) {
463 lp_->applyOp( *P_, *AP_ );
466 MVT::MvDot( *P_, *AP_, pAp );
470 for (i=0; i<numRHS_; ++i) {
471 if ( assertPositiveDefiniteness_ )
475 "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
477 alpha(i,i) = rHz[i] / pAp[i];
486 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
487 lp_->updateSolution();
490 MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
495 for (i=0; i<numRHS_; ++i) {
501 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
507 lp_->applyLeftPrec( *R_, *Z_ );
510 lp_->applyRightPrec( *Z_, *tmp );
515 lp_->applyRightPrec( *R_, *Z_ );
521 MVT::MvDot( *R_, *Z_, rHz );
522 if ( assertPositiveDefiniteness_ )
523 for (i=0; i<numRHS_; ++i)
526 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
529 for (i=0; i<numRHS_; ++i) {
530 beta(i,i) = rHz[i] / rHz_old[i];
534 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
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.
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
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.
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.
Teuchos::SerialDenseVector< int, ScalarType > & normal()
Wrapper for Normal(0,1) random variables.
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.
bool assertPositiveDefiniteness_
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
Teuchos::SerialDenseVector< int, ScalarType > randvec_
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Structure to contain pointers to CGIteration state variables.