42 #ifndef BELOS_BLOCK_GCRODR_ITER_HPP 43 #define BELOS_BLOCK_GCRODR_ITER_HPP 94 template <
class ScalarType,
class MV>
166 template<
class ScalarType,
class MV,
class OP>
321 if (!initialized_)
return 0;
349 void setSize(
int recycledBlocks,
int numBlocks ) {
351 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
352 recycledBlocks_ = recycledBlocks;
353 numBlocks_ = numBlocks;
383 int numBlocks_, blockSize_;
388 std::vector<bool> trueRHSIndices_;
402 std::vector< SDM >House_;
414 int curDim_, iter_, lclIter_;
462 template<
class ScalarType,
class MV,
class OP>
468 om_(printer), stest_(tester), ortho_(ortho) {
472 initialized_ =
false;
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;
516 House_.resize(numBlocks_);
518 for(i=0; i<numBlocks_;i++){
519 House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
525 template <
class ScalarType,
class MV,
class OP>
533 setSize( recycledBlocks_, numBlocks_ );
537 std::vector<int> curind(blockSize_);
543 for(
int i = 0; i<blockSize_; i++){curind[i] = i;};
548 Vnext = MVT::CloneViewNonConst(*V_,curind);
553 int rank = ortho_->normalize(*Vnext,Z0);
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;
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);
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);
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;
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>
694 if (static_cast<int> (norms->size()) < blockSize_) {
695 norms->resize( blockSize_ );
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>
723 if(curDim_<=blockSize_) {
724 return currentUpdate;
730 currentUpdate = MVT::Clone( *V_, blockSize_ );
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;
764 MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
771 if (U_ != Teuchos::null) {
772 SDM z(recycledBlocks_,blockSize_);
777 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
783 return currentUpdate;
786 template<
class ScalarType,
class MV,
class OP>
798 int curDim = curDim_;
799 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
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] );
844 int R_colStart = curDim_-blockSize_;
850 for(i=0; i<lclIter_-1; i++){
851 int R_rowStart = i*blockSize_;
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;
875 ScalarType nvs = blas.
NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
876 ScalarType alpha = -signDiag*nvs;
895 (*v_refl)[0] -= alpha;
896 nvs = blas.
NRM2(blockSize_+1,v_refl -> values(),1);
912 if(i < blockSize_ - 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());
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());
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...
Array< T > & append(const T &x)
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
ScalarType * values() const
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
void GER(const OrdinalType m, const OrdinalType n, const alpha_type alpha, const x_type *x, const OrdinalType incx, const y_type *y, const OrdinalType incy, ScalarType *A, const OrdinalType lda) const
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 TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb) const
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.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, const x_type *x, const OrdinalType incx, const beta_type beta, ScalarType *y, const OrdinalType incy) const
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
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.
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Class which defines basic traits for the operator type.
Teuchos::RCP< MV > V
The current Krylov basis.
Traits class which defines basic operations on multivectors.
int curDim
The current dimension of the reduction.
SCT::magnitudeType MagnitudeType
void ROT(const OrdinalType n, ScalarType *dx, const OrdinalType incx, ScalarType *dy, const OrdinalType incy, MagnitudeType *c, ScalarType *s) const
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...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
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
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.
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 isParameter(const std::string &name) const
bool isInitialized()
States whether the solver has been initialized or not.
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.
OrdinalType stride() const
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
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.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos::ScalarTraits< ScalarType > SCT
OperatorTraits< ScalarType, MV, OP > OPT
int sizeUninitialized(OrdinalType length_in)
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...