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_ );
249 for (
int i=0; i<numRHS_; i++)
252 r=0.5*SCT::real(randvec_[i]) + 1.0;
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_;
295 bool assertPositiveDefiniteness_;
310 Teuchos::RCP<MV> AP_;
316 Teuchos::SerialDenseVector<int,ScalarType> randvec_;
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.");
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 );
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 ...