163 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
181 Teuchos::ParameterList ¶ms );
295 if (!initialized_)
return 0;
328 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
332 void setSize(
int recycledBlocks,
int numBlocks ) {
334 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
335 recycledBlocks_ = recycledBlocks;
336 numBlocks_ = numBlocks;
337 cs_.sizeUninitialized( numBlocks_+1 );
338 sn_.sizeUninitialized( numBlocks_+1 );
339 z_.sizeUninitialized( numBlocks_+1 );
340 R_.shapeUninitialized( numBlocks_+1,numBlocks );
358 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
359 const Teuchos::RCP<OutputManager<ScalarType> > om_;
360 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
361 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
373 Teuchos::SerialDenseVector<int,ScalarType> sn_;
374 Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
394 Teuchos::RCP<MV> U_, C_;
398 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
401 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
406 Teuchos::SerialDenseMatrix<int,ScalarType> R_;
407 Teuchos::SerialDenseVector<int,ScalarType> z_;
417 Teuchos::ParameterList ¶ms ):
418 lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
421 initialized_ =
false;
431 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
"Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
432 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
434 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
435 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
437 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument,
"Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
438 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument,
"Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
441 recycledBlocks_ = rb;
443 cs_.sizeUninitialized( numBlocks_+1 );
444 sn_.sizeUninitialized( numBlocks_+1 );
445 z_.sizeUninitialized( numBlocks_+1 );
446 R_.shapeUninitialized( numBlocks_+1,numBlocks_ );
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.