42 #ifndef BELOS_BLOCK_GCRODR_ITER_HPP 43 #define BELOS_BLOCK_GCRODR_ITER_HPP 60 #include "Teuchos_BLAS.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_ScalarTraits.hpp" 64 #include "Teuchos_ParameterList.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 94 template <
class ScalarType,
class MV>
107 Teuchos::RCP<MV>
U,
C;
114 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H;
118 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B;
166 template<
class ScalarType,
class MV,
class OP>
175 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
177 typedef Teuchos::SerialDenseMatrix<int,ScalarType>
SDM;
178 typedef Teuchos::SerialDenseVector<int,ScalarType>
SDV;
195 Teuchos::ParameterList ¶ms );
349 void setSize(
int recycledBlocks,
int numBlocks ) {
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_;
395 Teuchos::SerialDenseVector<int,MagnitudeType>
cs_;
434 Teuchos::RCP<SDM >
H_;
439 Teuchos::RCP<SDM >
B_;
447 Teuchos::RCP<SDM>
R_;
462 template<
class ScalarType,
class MV,
class OP>
467 Teuchos::ParameterList ¶ms ):lp_(problem),
468 om_(printer), stest_(tester), ortho_(ortho) {
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.");
525 template <
class ScalarType,
class MV,
class OP>
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_;
650 template <
class ScalarType,
class MV,
class OP>
652 if (newstate.
V != Teuchos::null && newstate.
H != Teuchos::null) {
653 curDim_ = newstate.
curDim;
660 R_ = Teuchos::rcp(
new SDM(H_->numRows(), H_->numCols() ));
664 SDM Identity(2*blockSize_, 2*blockSize_);
665 for(
int i=0;i<2*blockSize_; i++){
668 for(
int i=0; i<numBlocks_;i++){
669 House_[i].assign(Identity);
673 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
V == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have V initialized.");
674 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
H == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have H initialized.");
686 template <
class ScalarType,
class MV,
class OP>
687 Teuchos::RCP<const MV>
694 if (static_cast<int> (norms->size()) < blockSize_) {
695 norms->resize( blockSize_ );
697 Teuchos::BLAS<int,ScalarType> blas;
698 for (
int j=0; j<blockSize_; j++) {
699 if(trueRHSIndices_[j]){
700 (*norms)[j] = blas.NRM2( blockSize_, &Z_(curDim_-blockSize_+j, j), 1);
706 return Teuchos::null;
709 return Teuchos::null;
715 template <
class ScalarType,
class MV,
class OP>
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;
786 template<
class ScalarType,
class MV,
class OP>
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;
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > V_
The Krylov basis vectors.
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *
virtual ~BlockGCRODRIter()
Destructor.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors...
void iterate()
This method performs block GCRODR iterations until the status test indicates the need to stop or an e...
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
Teuchos::RCP< SDM > B_
Projected matrix from the recycled subspace.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
A pure virtual class for defining the status tests for the Belos iterative solvers.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Class which defines basic traits for the operator type.
Teuchos::RCP< MV > V
The current Krylov basis.
std::vector< SDM > House_
Traits class which defines basic operations on multivectors.
int curDim
The current dimension of the reduction.
Teuchos::RCP< SDM > H_
Projected matrix from the Krylov factorization.
SCT::magnitudeType MagnitudeType
BlockGCRODRIterOrthoFailure(const std::string &what_arg)
BlockGCRODRIter(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)
BlockGCRODRIter constructor with linear problem, solver utilities, and parameter list of solver optio...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::SerialDenseVector< int, MagnitudeType > cs_
BlockGCRODRIterInitFailure(const std::string &what_arg)
A linear system to solve, and its associated information.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void updateLSQR(int dim=-1)
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< MV > U_
Recycled subspace vectors.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
void resetNumIters(int iter=0)
Reset the iteration count.
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
void initialize()
Initialize the solver to an iterate, providing a complete state.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
BlockGCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
SDM Z_
Q applied to right-hand side of the least squares system.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
BlockGCRODRIterInitFailure is thrown when the BlockGCRODRIter object is unable to generate an initial...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
std::vector< bool > trueRHSIndices_
Class which defines basic traits for the operator type.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Parent class to all Belos exceptions.
int getRecycledBlocks() const
Set the maximum number of recycled blocks used by the iterative solver.
Teuchos::SerialDenseVector< int, ScalarType > SDV
const Teuchos::RCP< OutputManager< ScalarType > > om_
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
OperatorTraits< ScalarType, MV, OP > OPT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::RCP< SDM > R_
Upper triangular reduction of H_ (see above).
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...