458 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
460 return currentUpdate;
462 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
463 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464 Teuchos::BLAS<int,ScalarType> blas;
465 currentUpdate = MVT::Clone( *V_, 1 );
469 Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, z_, curDim_, 1 );
473 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
474 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
475 R_.values(), R_.stride(), y.values(), y.stride() );
479 std::vector<int> index(curDim_);
480 for (
int i=0; i<curDim_; i++ ) index[i] = i;
481 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
482 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
486 if (U_ != Teuchos::null) {
487 Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
488 Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
489 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
490 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
493 return currentUpdate;
546 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
GCRODRIterInitFailure,
"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
549 setSize( recycledBlocks_, numBlocks_ );
551 Teuchos::RCP<MV> Vnext;
552 Teuchos::RCP<const MV> Vprev;
553 std::vector<int> curind(1);
560 Vnext = MVT::CloneViewNonConst(*V_,curind);
563 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z0 =
564 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(1,1) );
565 int rank = ortho_->normalize( *Vnext, z0 );
566 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
570 std::vector<int> prevind(numBlocks_+1);
577 if (U_ == Teuchos::null) {
578 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= numBlocks_) {
580 int lclDim = curDim_ + 1;
584 Vnext = MVT::CloneViewNonConst(*V_,curind);
588 Vprev = MVT::CloneView(*V_,curind);
591 lp_->apply(*Vprev,*Vnext);
596 prevind.resize(lclDim);
597 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
598 Vprev = MVT::CloneView(*V_,prevind);
599 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
602 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
603 subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
604 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
605 AsubH.append( subH );
608 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
609 subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
612 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
615 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
616 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
619 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
629 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= numBlocks_) {
631 int lclDim = curDim_ + 1;
635 Vnext = MVT::CloneViewNonConst(*V_,curind);
640 Vprev = MVT::CloneView(*V_,curind);
643 lp_->apply(*Vprev,*Vnext);
644 Vprev = Teuchos::null;
647 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
648 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
649 subB = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
650 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
651 AsubB.append( subB );
654 ortho_->project( *Vnext, AsubB, C );
658 prevind.resize(lclDim);
659 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
660 Vprev = MVT::CloneView(*V_,prevind);
661 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
664 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
665 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
666 AsubH.append( subH );
669 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
672 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
675 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
676 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
679 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
697 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
702 int curDim = curDim_;
703 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
706 Teuchos::BLAS<int, ScalarType> blas;
713 for (i=0; i<curDim; i++) {
717 blas.ROT( 1, &R_(i,curDim), 1, &R_(i+1, curDim), 1, &cs_[i], &sn_[i] );
722 blas.ROTG( &R_(curDim,curDim), &R_(curDim+1,curDim), &cs_[curDim], &sn_[curDim] );
723 R_(curDim+1,curDim) = zero;
727 blas.ROT( 1, &z_(curDim), 1, &z_(curDim+1), 1, &cs_[curDim], &sn_[curDim] );
GCRODRIter(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< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
GCRODRIter constructor with linear problem, solver utilities, and parameter list of solver options.