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 );
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)) {
367 std::invalid_argument, errstr );
369 std::invalid_argument, errstr );
372 if (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_ );
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();
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_ );
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]);
474 lp_->updateSolution();
482 for (i=0; i<numRHS_; ++i) {
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_ );
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];
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
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
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.
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
virtual ~PseudoBlockStochasticCGIter()
Destructor.
int getNumIters() const
Get the current iteration count.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Traits class which defines basic operations on multivectors.
void resetNumIters(int iter=0)
Reset the iteration count.
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the 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.
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.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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.