316 if (!stateStorageInitialized_) {
319 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
320 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
321 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
322 stateStorageInitialized_ =
false;
329 if (R_ == Teuchos::null) {
331 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
332 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
333 "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
336 W_ = MVT::Clone( *tmp, 3 );
337 std::vector<int> index2(2,0);
338 std::vector<int> index(1,0);
343 S_ = MVT::CloneViewNonConst( *W_, index2 );
348 U_ = MVT::CloneViewNonConst( *W_, index2 );
351 R_ = MVT::CloneViewNonConst( *W_, index );
353 AZ_ = MVT::CloneViewNonConst( *W_, index );
355 Z_ = MVT::CloneViewNonConst( *W_, index );
360 T_ = MVT::CloneViewNonConst( *W_, index2 );
363 V_ = MVT::Clone( *tmp, 2 );
365 AP_ = MVT::CloneViewNonConst( *V_, index );
367 P_ = MVT::CloneViewNonConst( *V_, index );
372 stateStorageInitialized_ =
true;
384 if (!stateStorageInitialized_)
387 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
388 "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
392 std::string errstr(
"Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
394 if (newstate.
R != Teuchos::null) {
396 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
397 std::invalid_argument, errstr );
398 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
399 std::invalid_argument, errstr );
402 if (newstate.
R != R_) {
404 MVT::Assign( *newstate.
R, *R_ );
410 if ( lp_->getLeftPrec() != Teuchos::null ) {
411 lp_->applyLeftPrec( *R_, *Z_ );
412 if ( lp_->getRightPrec() != Teuchos::null ) {
413 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
414 lp_->applyRightPrec( *Z_, *tmp );
415 MVT::Assign( *tmp, *Z_ );
418 else if ( lp_->getRightPrec() != Teuchos::null ) {
419 lp_->applyRightPrec( *R_, *Z_ );
422 MVT::Assign( *R_, *Z_ );
426 lp_->applyOp( *Z_, *AZ_ );
430 MVT::Assign( *U_, *V_);
434 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
435 "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
467 if (initialized_ ==
false) {
472 Teuchos::SerialDenseMatrix<int,ScalarType> sHz( 2, 1 );
473 Teuchos::SerialDenseMatrix<int,ScalarType> sHt( 2, 2 );
474 ScalarType rHz_old, alpha, beta, delta;
477 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
478 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
481 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
484 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
485 "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
487 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
489 MVT::MvTransMv( one, *S_, *T_, sHt );
495 MVT::MvTransMv( one, *S_, *Z_, sHz );
499 if ((Teuchos::ScalarTraits<ScalarType>::magnitude(delta) < Teuchos::ScalarTraits<ScalarType>::eps()) &&
500 (stest_->checkStatus(
this) ==
Passed))
502 alpha = rHz_ / delta;
506 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
511 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
519 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
523 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
527 if ( lp_->getLeftPrec() != Teuchos::null ) {
528 lp_->applyLeftPrec( *R_, *Z_ );
529 if ( lp_->getRightPrec() != Teuchos::null ) {
530 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
531 lp_->applyRightPrec( *tmp, *Z_ );
534 else if ( lp_->getRightPrec() != Teuchos::null ) {
535 lp_->applyRightPrec( *R_, *Z_ );
538 MVT::Assign( *R_, *Z_ );
542 lp_->applyOp( *Z_, *AZ_ );
545 MVT::MvTransMv( one, *S_, *T_, sHt );
558 if (stest_->checkStatus(
this) ==
Passed) {
562 beta = rHz_ / rHz_old;
563 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
567 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
574 MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
585 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
589 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
596 if (stest_->checkStatus(
this) ==
Passed) {
602 if ( lp_->getLeftPrec() != Teuchos::null ) {
603 lp_->applyLeftPrec( *R_, *Z_ );
604 if ( lp_->getRightPrec() != Teuchos::null ) {
605 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
606 lp_->applyRightPrec( *tmp, *Z_ );
609 else if ( lp_->getRightPrec() != Teuchos::null ) {
610 lp_->applyRightPrec( *R_, *Z_ );
613 MVT::Assign( *R_, *Z_ );
617 lp_->applyOp( *Z_, *AZ_ );
620 MVT::MvTransMv( one, *S_, *Z_, sHz );
627 beta = rHz_ / rHz_old;
628 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
632 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
639 MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
CGSingleRedIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > &convTester, Teuchos::ParameterList ¶ms)
CGSingleRedIter constructor with linear problem, solver utilities, and parameter list of solver optio...