42 #ifndef BELOS_BLOCK_GMRES_ITER_HPP
43 #define BELOS_BLOCK_GMRES_ITER_HPP
83 template<
class ScalarType,
class MV,
class OP>
226 if (!initialized_)
return 0;
260 void setSize(
int blockSize,
int numBlocks);
275 CheckList() : checkV(false), checkArn(false) {};
281 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
318 bool stateStorageInitialized_;
323 bool keepHessenberg_;
327 bool initHessenberg_;
351 template<
class ScalarType,
class MV,
class OP>
364 stateStorageInitialized_(false),
365 keepHessenberg_(false),
366 initHessenberg_(false),
371 if ( om_->isVerbosity(
Debug ) )
372 keepHessenberg_ =
true;
374 keepHessenberg_ = params.
get(
"Keep Hessenberg",
false);
377 initHessenberg_ = params.
get(
"Initialize Hessenberg",
false);
381 "Belos::BlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
382 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
385 int bs = params.get(
"Block Size", 1);
391 template <
class ScalarType,
class MV,
class OP>
397 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument,
"Belos::BlockGmresIter::setSize was passed a non-positive argument.");
398 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
403 if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
404 stateStorageInitialized_ =
false;
406 blockSize_ = blockSize;
407 numBlocks_ = numBlocks;
409 initialized_ =
false;
419 template <
class ScalarType,
class MV,
class OP>
422 if (!stateStorageInitialized_) {
427 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
428 stateStorageInitialized_ =
false;
436 int newsd = blockSize_*(numBlocks_+1);
443 beta.resize( newsd );
447 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_*static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
448 "Belos::BlockGmresIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
451 if (V_ == Teuchos::null) {
455 "Belos::BlockGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
456 V_ = MVT::Clone( *tmp, newsd );
460 if (MVT::GetNumberVecs(*V_) < newsd) {
462 V_ = MVT::Clone( *tmp, newsd );
467 if (R_ == Teuchos::null) {
470 if (initHessenberg_) {
471 R_->shape( newsd, newsd-blockSize_ );
474 if (R_->numRows() < newsd || R_->numCols() < newsd-blockSize_) {
475 R_->shapeUninitialized( newsd, newsd-blockSize_ );
480 if (keepHessenberg_) {
481 if (H_ == Teuchos::null) {
484 if (initHessenberg_) {
485 H_->shape( newsd, newsd-blockSize_ );
488 if (H_->numRows() < newsd || H_->numCols() < newsd-blockSize_) {
489 H_->shapeUninitialized( newsd, newsd-blockSize_ );
499 if (z_ == Teuchos::null) {
502 if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
503 z_->shapeUninitialized( newsd, blockSize_ );
507 stateStorageInitialized_ =
true;
514 template <
class ScalarType,
class MV,
class OP>
523 return currentUpdate;
525 const ScalarType one = SCT::one();
526 const ScalarType zero = SCT::zero();
528 currentUpdate = MVT::Clone( *V_, blockSize_ );
542 std::vector<int> index(curDim_);
543 for (
int i=0; i<curDim_; i++ ) {
547 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
549 return currentUpdate;
556 template <
class ScalarType,
class MV,
class OP>
562 if ( norms && (
int)norms->size() < blockSize_ )
563 norms->resize( blockSize_ );
567 for (
int j=0; j<blockSize_; j++) {
568 (*norms)[j] = blas.
NRM2( blockSize_, &(*z_)(curDim_, j), 1);
571 return Teuchos::null;
578 template <
class ScalarType,
class MV,
class OP>
582 if (!stateStorageInitialized_)
586 "Belos::BlockGmresIter::initialize(): Cannot initialize state storage!");
588 const ScalarType zero = SCT::zero(), one = SCT::one();
594 std::string errstr(
"Belos::BlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
596 if (newstate.
V != Teuchos::null && newstate.
z != Teuchos::null) {
601 std::invalid_argument, errstr );
603 std::invalid_argument, errstr );
605 std::invalid_argument, errstr );
607 curDim_ = newstate.
curDim;
608 int lclDim = MVT::GetNumberVecs(*newstate.
V);
615 if (newstate.
V != V_) {
617 if (curDim_ == 0 && lclDim > blockSize_) {
618 om_->stream(
Warnings) <<
"Belos::BlockGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
619 <<
"The block size however is only " << blockSize_ << std::endl
620 <<
"The last " << lclDim - blockSize_ <<
" vectors will be discarded." << std::endl;
622 std::vector<int> nevind(curDim_+blockSize_);
623 for (
int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
626 MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
629 lclV = Teuchos::null;
633 if (newstate.
z != z_) {
641 lclZ = Teuchos::null;
648 "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
651 "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial norms z_0.");
657 if (om_->isVerbosity(
Debug ) ) {
662 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
670 template <
class ScalarType,
class MV,
class OP>
676 if (initialized_ ==
false) {
681 int searchDim = blockSize_*numBlocks_;
688 while (stest_->checkStatus(
this) !=
Passed && curDim_+blockSize_ <= searchDim) {
693 int lclDim = curDim_ + blockSize_;
696 std::vector<int> curind(blockSize_);
697 for (
int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
702 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
706 lp_->apply(*Vprev,*Vnext);
707 Vprev = Teuchos::null;
711 std::vector<int> prevind(lclDim);
712 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
713 Vprev = MVT::CloneView(*V_,prevind);
726 (
Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
728 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
732 if (keepHessenberg_) {
742 (
Teuchos::View,*R_,blockSize_,blockSize_,lclDim,curDim_ ) );
743 subR2->assign(*subH2);
747 "Belos::BlockGmresIter::iterate(): couldn't generate basis of full rank.");
757 Vnext = Teuchos::null;
758 curDim_ += blockSize_;
761 if (om_->isVerbosity(
Debug ) ) {
766 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
771 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );
779 template<
class ScalarType,
class MV,
class OP>
783 ScalarType sigma, mu, vscale, maxelem;
784 const ScalarType zero = SCT::zero();
789 int curDim = curDim_;
790 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
799 if (blockSize_ == 1) {
803 for (i=0; i<curDim; i++) {
807 blas.
ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
812 blas.
ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
813 (*R_)(curDim+1,curDim) = zero;
817 blas.
ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
823 for (j=0; j<blockSize_; j++) {
827 for (i=0; i<curDim+j; i++) {
828 sigma = blas.
DOT( blockSize_, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
829 sigma += (*R_)(i,curDim+j);
830 sigma *= SCT::conjugate(beta[i]);
831 blas.
AXPY(blockSize_, ScalarType(-sigma), &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
832 (*R_)(i,curDim+j) -= sigma;
837 maxidx = blas.
IAMAX( blockSize_+1, &(*R_)(curDim+j,curDim+j), 1 );
838 maxelem = SCT::magnitude((*R_)(curDim+j+maxidx-1,curDim+j));
839 for (i=0; i<blockSize_+1; i++)
840 (*R_)(curDim+j+i,curDim+j) /= maxelem;
841 sigma = blas.
DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j), 1,
842 &(*R_)(curDim+j+1,curDim+j), 1 );
843 MagnitudeType sign_Rjj = -SCT::real((*R_)(curDim+j,curDim+j)) /
844 SCT::magnitude(SCT::real(((*R_)(curDim+j,curDim+j))));
846 beta[curDim + j] = zero;
848 mu = SCT::squareroot(SCT::conjugate((*R_)(curDim+j,curDim+j))*(*R_)(curDim+j,curDim+j)+sigma);
849 vscale = (*R_)(curDim+j,curDim+j) - Teuchos::as<ScalarType>(sign_Rjj)*mu;
850 beta[curDim+j] = -Teuchos::as<ScalarType>(sign_Rjj) * vscale / mu;
851 (*R_)(curDim+j,curDim+j) = Teuchos::as<ScalarType>(sign_Rjj)*maxelem*mu;
852 for (i=0; i<blockSize_; i++)
853 (*R_)(curDim+j+1+i,curDim+j) /= vscale;
858 for (i=0; i<blockSize_; i++) {
859 sigma = blas.
DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j),
860 1, &(*z_)(curDim+j+1,i), 1);
861 sigma += (*z_)(curDim+j,i);
862 sigma *= SCT::conjugate(beta[curDim+j]);
863 blas.
AXPY(blockSize_, ScalarType(-sigma), &(*R_)(curDim+j+1,curDim+j),
864 1, &(*z_)(curDim+j+1,i), 1);
865 (*z_)(curDim+j,i) -= sigma;
871 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
872 curDim_ = dim + blockSize_;
888 template <
class ScalarType,
class MV,
class OP>
891 std::stringstream os;
893 os.setf(std::ios::scientific, std::ios::floatfield);
896 os <<
" Debugging checks: iteration " << iter_ << where << std::endl;
899 std::vector<int> lclind(curDim_);
900 for (
int i=0; i<curDim_; i++) lclind[i] = i;
901 std::vector<int> bsind(blockSize_);
902 for (
int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
907 lclV = MVT::CloneView(*V_,lclind);
908 lclF = MVT::CloneView(*V_,bsind);
912 tmp = ortho_->orthonormError(*lclV);
913 os <<
" >> Error in V^H M V == I : " << tmp << std::endl;
915 tmp = ortho_->orthonormError(*lclF);
916 os <<
" >> Error in F^H M F == I : " << tmp << std::endl;
918 tmp = ortho_->orthogError(*lclV,*lclF);
919 os <<
" >> Error in V^H M F == 0 : " << tmp << std::endl;
927 lclAV = MVT::Clone(*V_,curDim_);
928 lp_->apply(*lclV,*lclAV);
933 MVT::MvTimesMatAddMv( -one, *lclV, subH, one, *lclAV );
937 blockSize_,curDim_, curDim_ );
938 MVT::MvTimesMatAddMv( -one, *lclF, curB, one, *lclAV );
941 std::vector<MagnitudeType> arnNorms( curDim_ );
942 ortho_->norm( *lclAV, arnNorms );
944 for (
int i=0; i<curDim_; i++) {
945 os <<
" >> Error in Krylov factorization (R = AV-VH-FB^H), ||R[" << i <<
"]|| : " << arnNorms[i] << std::endl;
ScalarType * values() const
Collection of types and exceptions used within the Belos solvers.
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
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Array< T > & append(const T &x)
Class which manages the output and verbosity of the Belos solvers.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
T & get(const std::string &name, T def_value)
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Declaration of basic traits for the multivector type.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
OrdinalType IAMAX(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Structure to contain pointers to GmresIteration state variables.
void resetNumIters(int iter=0)
Reset the iteration count.
SCT::magnitudeType MagnitudeType
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this linear prob...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
virtual ~BlockGmresIter()
Destructor.
Traits class which defines basic operations on multivectors.
bool isParameter(const std::string &name) const
bool isInitialized()
States whether the solver has been initialized or not.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int curDim
The current dimension of the reduction.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
ScalarType DOT(const OrdinalType &n, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy) const
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Teuchos::ScalarTraits< ScalarType > SCT
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
OrdinalType numCols() const
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.
void initializeGmres(GmresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
BlockGmresIter(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)
BlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver option...
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
void setBlockSize(int blockSize)
Set the blocksize.
OrdinalType stride() const
OrdinalType numRows() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.