175 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
177 typedef Teuchos::SerialDenseMatrix<int,ScalarType>
SDM;
178 typedef Teuchos::SerialDenseVector<int,ScalarType>
SDV;
195 Teuchos::ParameterList ¶ms );
321 if (!initialized_)
return 0;
349 void setSize(
int recycledBlocks,
int numBlocks ) {
351 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
352 recycledBlocks_ = recycledBlocks;
353 numBlocks_ = numBlocks;
354 cs_.sizeUninitialized( numBlocks_ );
355 sn_.sizeUninitialized( numBlocks_ );
356 Z_.shapeUninitialized( numBlocks_*blockSize_, blockSize_ );
372 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
373 const Teuchos::RCP<OutputManager<ScalarType> > om_;
374 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
375 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
383 int numBlocks_, blockSize_;
388 std::vector<bool> trueRHSIndices_;
395 Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
402 std::vector< SDM >House_;
414 int curDim_, iter_, lclIter_;
425 Teuchos::RCP<MV> U_, C_;
434 Teuchos::RCP<SDM > H_;
439 Teuchos::RCP<SDM > B_;
447 Teuchos::RCP<SDM> R_;
467 Teuchos::ParameterList ¶ms ):lp_(problem),
468 om_(printer), stest_(tester), ortho_(ortho) {
472 initialized_ =
false;
483 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
484 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
486 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
487 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
489 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument,
"Belos::BlockGCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
490 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument,
"Belos::BlockGCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
493 int bs = Teuchos::getParameter<int>(params,
"Block Size");
495 TEUCHOS_TEST_FOR_EXCEPTION(bs <= 0, std::invalid_argument,
"Belos::BlockGCRODRIter() the block size was passed a non-postitive argument.");
499 recycledBlocks_ = rb;
504 trueRHSIndices_.resize(blockSize_);
506 for(i=0; i<blockSize_; i++){
507 trueRHSIndices_[i] =
true;
512 cs_.sizeUninitialized( numBlocks_+1 );
513 sn_.sizeUninitialized( numBlocks_+1 );
514 Z_.shapeUninitialized( (numBlocks_+1)*blockSize_,blockSize_ );
516 House_.resize(numBlocks_);
518 for(i=0; i<numBlocks_;i++){
519 House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
527 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
BlockGCRODRIterInitFailure,
"Belos::BlockGCRODRIter::iterate(): GCRODRIter class not initialized." );
533 setSize( recycledBlocks_, numBlocks_ );
535 Teuchos::RCP<MV> Vnext;
536 Teuchos::RCP<const MV> Vprev;
537 std::vector<int> curind(blockSize_);
543 for(
int i = 0; i<blockSize_; i++){curind[i] = i;};
548 Vnext = MVT::CloneViewNonConst(*V_,curind);
551 Teuchos::RCP<SDM > Z0 =
552 Teuchos::rcp(
new SDM(blockSize_,blockSize_) );
553 int rank = ortho_->normalize(*Vnext,Z0);
558 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,
BlockGCRODRIterOrthoFailure,
"Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank at the initial step.");
560 Teuchos::RCP<SDM > Z_block = Teuchos::rcp(
new SDM(Teuchos::View, Z_, blockSize_,blockSize_) );
561 Z_block->assign(*Z0);
563 std::vector<int> prevind(blockSize_*(numBlocks_ + 1));
570 while( (stest_->checkStatus(
this) !=
Passed) && (curDim_+blockSize_-1) < (numBlocks_*blockSize_)) {
576 int HFirstCol = curDim_-blockSize_;
577 int HLastCol = HFirstCol + blockSize_-1 ;
578 int HLastOrthRow = HLastCol;
579 int HFirstNormRow = HLastOrthRow + 1;
583 for(
int i = 0; i< blockSize_; i++){
584 curind[i] = curDim_ + i;
586 Vnext = MVT::CloneViewNonConst(*V_,curind);
591 for(
int i = 0; i< blockSize_; i++){
592 curind[blockSize_ - 1 - i] = curDim_ - i - 1;
594 Vprev = MVT::CloneView(*V_,curind);
596 lp_->apply(*Vprev,*Vnext);
597 Vprev = Teuchos::null;
603 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
606 subB = Teuchos::rcp(
new SDM ( Teuchos::View,*B_,recycledBlocks_,blockSize_,0, HFirstCol ) );
608 Teuchos::Array<Teuchos::RCP<SDM > > AsubB;
609 AsubB.append( subB );
611 ortho_->project( *Vnext, AsubB, C );
615 prevind.resize(curDim_);
616 for (
int i=0; i<curDim_; i++) { prevind[i] = i; }
617 Vprev = MVT::CloneView(*V_,prevind);
618 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
621 Teuchos::RCP<SDM> subH = Teuchos::rcp(
new SDM ( Teuchos::View,*H_,curDim_,blockSize_,0,HFirstCol ) );
622 Teuchos::Array<Teuchos::RCP<SDM > > AsubH;
623 AsubH.append( subH );
625 Teuchos::RCP<SDM > subR = Teuchos::rcp(
new SDM ( Teuchos::View,*H_,blockSize_,blockSize_,HFirstNormRow,HFirstCol ) );
627 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
630 SDM subR2( Teuchos::View,*R_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
631 SDM subH2( Teuchos::View,*H_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
634 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,
BlockGCRODRIterOrthoFailure,
"Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank.");
640 curDim_ = curDim_ + blockSize_;
721 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
723 if(curDim_<=blockSize_) {
724 return currentUpdate;
727 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
728 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
729 Teuchos::BLAS<int,ScalarType> blas;
730 currentUpdate = MVT::Clone( *V_, blockSize_ );
734 SDM Y( Teuchos::Copy, Z_, curDim_-blockSize_, blockSize_ );
735 Teuchos::RCP<SDM> Rtmp = Teuchos::rcp(
new SDM(Teuchos::View, *R_, curDim_, curDim_-blockSize_));
752 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
753 Teuchos::NON_UNIT_DIAG, curDim_-blockSize_, blockSize_, one,
754 Rtmp->values(), Rtmp->stride(), Y.values(), Y.stride() );
761 std::vector<int> index(curDim_-blockSize_);
762 for (
int i=0; i<curDim_-blockSize_; i++ ) index[i] = i;
763 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
764 MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
771 if (U_ != Teuchos::null) {
772 SDM z(recycledBlocks_,blockSize_);
773 SDM subB( Teuchos::View, *B_, recycledBlocks_, curDim_-blockSize_ );
774 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, Y, zero );
777 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
783 return currentUpdate;
790 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
791 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
798 int curDim = curDim_;
799 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
803 Teuchos::BLAS<int, ScalarType> blas;
812 for (i=0; i<curDim-1; i++) {
816 blas.ROT( 1, &(*R_)(i,curDim-1), 1, &(*R_)(i+1, curDim-1), 1, &cs_[i], &sn_[i] );
821 blas.ROTG( &(*R_)(curDim-1,curDim-1), &(*R_)(curDim,curDim-1), &cs_[curDim-1], &sn_[curDim-1] );
822 (*R_)(curDim,curDim-1) = zero;
826 blas.ROT( 1, &Z_(curDim-1,0), 1, &Z_(curDim,0), 1, &cs_[curDim-1], &sn_[curDim-1] );
841 Teuchos::RCP< SDM > workmatrix = Teuchos::null;
842 Teuchos::RCP< SDV > workvec = Teuchos::null;
843 Teuchos::RCP<SDV> v_refl = Teuchos::null;
844 int R_colStart = curDim_-blockSize_;
845 Teuchos::RCP< SDM >Rblock = Teuchos::null;
850 for(i=0; i<lclIter_-1; i++){
851 int R_rowStart = i*blockSize_;
853 Teuchos::RCP< SDM > RblockCopy = rcp(
new SDM (Teuchos::Copy, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
854 Teuchos::RCP< SDM > RblockView = rcp(
new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
855 blas.GEMM(Teuchos::NO_TRANS,Teuchos::NO_TRANS, 2*blockSize_,blockSize_,2*blockSize_,one,House_[i].values(),House_[i].stride(), RblockCopy->values(),RblockCopy -> stride(), zero, RblockView->values(),RblockView -> stride());
862 Rblock = rcp(
new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, curDim_-blockSize_, curDim_-blockSize_));
865 for(i=0; i<blockSize_; i++){
869 int curcol = (lclIter_ - 1)*blockSize_ + i;
871 ScalarType signDiag = (*R_)(curcol,curcol) / Teuchos::ScalarTraits<ScalarType>::magnitude((*R_)(curcol,curcol));
875 ScalarType nvs = blas.NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
876 ScalarType alpha = -signDiag*nvs;
894 v_refl = rcp(
new SDV(Teuchos::Copy, &((*R_)(curcol,curcol)), blockSize_ + 1 ));
895 (*v_refl)[0] -= alpha;
896 nvs = blas.NRM2(blockSize_+1,v_refl -> values(),1);
912 if(i < blockSize_ - 1){
913 workvec = Teuchos::rcp(
new SDV(blockSize_ - i -1));
915 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, *Rblock, blockSize_+1, blockSize_ - i -1, lclCurcol, lclCurcol +1 ) );
916 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix->values(), workmatrix->stride(), v_refl->values(), 1, zero, workvec->values(), 1);
917 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl->values(),1,workvec->values(),1,workmatrix->values(),workmatrix->stride());
924 workvec = Teuchos::rcp(
new SDV(2*blockSize_));
925 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, House_[lclIter_ -1], blockSize_+1, 2*blockSize_, i, 0 ) );
926 blas.GEMV(Teuchos::TRANS,workmatrix->numRows(),workmatrix->numCols(),one,workmatrix->values(),workmatrix->stride(), v_refl->values(), 1,zero,workvec->values(),1);
927 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl -> values(),1,workvec->values(),1,workmatrix->values(),(*workmatrix).stride());
932 workvec = Teuchos::rcp(
new SDV(blockSize_));
933 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, Z_, blockSize_+1, blockSize_, curcol, 0 ) );
934 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix-> values(), workmatrix->stride(), v_refl -> values(), 1, zero, workvec->values(), 1);
935 blas.GER((*workmatrix).numRows(),(*workmatrix).numCols(), -2.0*one/nvs,v_refl -> values(), 1,&((*workvec)[0]),1,(*workmatrix)[0],(*workmatrix).stride());
940 (*R_)[curcol][curcol] = alpha;
941 for(
int ii=1; ii<= blockSize_; ii++){
942 (*R_)[curcol][curcol+ii] = 0;