42#ifndef BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
43#define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_SerialDenseHelpers.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 Teuchos::SerialDenseVector<int,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 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);
273 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
274 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
275 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
316 Teuchos::SerialDenseVector<int,ScalarType>
randvec_;
322 template<
class ScalarType,
class MV,
class OP>
326 Teuchos::ParameterList ¶ms ):
333 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) )
340 template <
class ScalarType,
class MV,
class OP>
344 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
345 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
346 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
347 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
350 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
353 int numRHS = MVT::GetNumberVecs(*tmp);
358 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
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.");
374 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
375 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
377 if (!Teuchos::is_null(newstate.
R)) {
379 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
380 std::invalid_argument, errstr );
381 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
382 std::invalid_argument, errstr );
385 if (newstate.
R != R_) {
387 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
393 if ( lp_->getLeftPrec() != Teuchos::null ) {
394 lp_->applyLeftPrec( *R_, *Z_ );
395 if ( lp_->getRightPrec() != Teuchos::null ) {
396 Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
397 lp_->applyRightPrec( *Z_, *tmp2 );
401 else if ( lp_->getRightPrec() != Teuchos::null ) {
402 lp_->applyRightPrec( *R_, *Z_ );
407 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
411 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
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_ );
436 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(numRHS_,numRHS_);
439 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
443 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
446 MVT::MvDot( *R_, *Z_, rHz );
448 if ( assertPositiveDefiniteness_ )
449 for (i=0; i<numRHS_; ++i)
450 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
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 );
468 Teuchos::SerialDenseVector<int,ScalarType>& z = normal();
470 for (i=0; i<numRHS_; ++i) {
471 if ( assertPositiveDefiniteness_ )
473 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
475 "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
477 alpha(i,i) = rHz[i] / pAp[i];
480 zeta(i,i) = z[i] / Teuchos::ScalarTraits<ScalarType>::squareroot(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_ );
506 if ( lp_->getLeftPrec() != Teuchos::null ) {
507 lp_->applyLeftPrec( *R_, *Z_ );
508 if ( lp_->getRightPrec() != Teuchos::null ) {
509 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
510 lp_->applyRightPrec( *Z_, *tmp );
514 else if ( lp_->getRightPrec() != Teuchos::null ) {
515 lp_->applyRightPrec( *R_, *Z_ );
521 MVT::MvDot( *R_, *Z_, rHz );
522 if ( assertPositiveDefiniteness_ )
523 for (i=0; i<numRHS_; ++i)
524 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
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];
532 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
533 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
534 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
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.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the stochastic pseudo-block CG iteration, where the basic stochastic CG algorit...
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::SerialDenseVector< int, ScalarType > & normal()
Wrapper for Normal(0,1) random variables.
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 ...
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::SerialDenseVector< int, ScalarType > randvec_
int getNumIters() const
Get the current iteration count.
Teuchos::ScalarTraits< ScalarType > SCT
SCT::magnitudeType MagnitudeType
void initializeCG(StochasticCGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
MultiVecTraits< ScalarType, MV > MVT
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector
void resetNumIters(int iter=0)
Reset the iteration count.
void iterate()
This method performs stochastic CG iterations on each linear system until the status test indicates t...
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
const Teuchos::RCP< OutputManager< ScalarType > > om_
bool assertPositiveDefiniteness_
void setBlockSize(int blockSize)
Set the blocksize.
virtual ~PseudoBlockStochasticCGIter()
Destructor.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > Y
The current stochastic recurrence vector.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > R
The current residual.