42#ifndef BELOS_BLOCK_GMRES_ITER_HPP
43#define BELOS_BLOCK_GMRES_ITER_HPP
60#include "Teuchos_BLAS.hpp"
61#include "Teuchos_LAPACK.hpp"
62#include "Teuchos_SerialDenseMatrix.hpp"
63#include "Teuchos_SerialDenseVector.hpp"
64#include "Teuchos_ScalarTraits.hpp"
65#include "Teuchos_ParameterList.hpp"
66#include "Teuchos_TimeMonitor.hpp"
83template<
class ScalarType,
class MV,
class OP>
93 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
112 Teuchos::ParameterList ¶ms );
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;
289 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
290 const Teuchos::RCP<OutputManager<ScalarType> > om_;
291 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
292 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
304 Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
305 Teuchos::SerialDenseVector<int,MagnitudeType> cs;
318 bool stateStorageInitialized_;
323 bool keepHessenberg_;
327 bool initHessenberg_;
340 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
345 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
346 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;
351 template<
class ScalarType,
class MV,
class OP>
356 Teuchos::ParameterList ¶ms ):
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);
380 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
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_) {
425 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
426 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
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) {
453 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
454 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
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) {
461 Teuchos::RCP<const MV> tmp = V_;
462 V_ = MVT::Clone( *tmp, newsd );
467 if (R_ == Teuchos::null) {
468 R_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
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) {
482 H_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
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) {
500 z_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
502 if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
503 z_->shapeUninitialized( newsd, blockSize_ );
507 stateStorageInitialized_ =
true;
514 template <
class ScalarType,
class MV,
class OP>
521 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
523 return currentUpdate;
525 const ScalarType one = SCT::one();
526 const ScalarType zero = SCT::zero();
527 Teuchos::BLAS<int,ScalarType> blas;
528 currentUpdate = MVT::Clone( *V_, blockSize_ );
532 Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
536 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
537 Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
538 R_->values(), R_->stride(), y.values(), y.stride() );
542 std::vector<int> index(curDim_);
543 for (
int i=0; i<curDim_; i++ ) {
546 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
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_ );
566 Teuchos::BLAS<int,ScalarType> blas;
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_)
585 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
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) {
600 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
V) != MVT::GetGlobalLength(*V_),
601 std::invalid_argument, errstr );
602 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
V) < blockSize_,
603 std::invalid_argument, errstr );
604 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim > blockSize_*(numBlocks_+1),
605 std::invalid_argument, errstr );
607 curDim_ = newstate.
curDim;
608 int lclDim = MVT::GetNumberVecs(*newstate.
V);
611 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
z->numRows() < curDim_ || newstate.
z->numCols() < blockSize_, std::invalid_argument, errstr);
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;
624 Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.
V, nevind );
625 Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_, nevind );
626 MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
629 lclV = Teuchos::null;
633 if (newstate.
z != z_) {
635 Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.
z,curDim_+blockSize_,blockSize_);
636 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclZ;
637 lclZ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
641 lclZ = Teuchos::null;
647 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
V == Teuchos::null,std::invalid_argument,
648 "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
650 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
z == Teuchos::null,std::invalid_argument,
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; }
698 Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
702 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
703 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
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);
714 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
717 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
718 subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
719 ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
720 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
721 AsubH.append( subH );
724 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
725 subH2 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
726 ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
728 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
732 if (keepHessenberg_) {
734 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
735 subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
736 ( Teuchos::View,*R_,lclDim,blockSize_,0,curDim_ ) );
740 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
741 subR2 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
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()) {
794 Teuchos::BLAS<int, ScalarType> blas;
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; }
904 Teuchos::RCP<const MV> lclV,lclF;
905 Teuchos::RCP<MV> lclAV;
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);
931 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
932 Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
933 MVT::MvTimesMatAddMv( -one, *lclV, subH, one, *lclAV );
936 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
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;
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
Teuchos::ScalarTraits< ScalarType > SCT
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SCT::magnitudeType MagnitudeType
virtual ~BlockGmresIter()
Destructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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...
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...
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
MultiVecTraits< ScalarType, MV > MVT
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
bool isInitialized()
States whether the solver has been initialized or not.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
void resetNumIters(int iter=0)
Reset the iteration count.
void setBlockSize(int blockSize)
Set the blocksize.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
void initializeGmres(GmresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
OperatorTraits< ScalarType, MV, OP > OPT
int getNumIters() const
Get the current iteration count.
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to GmresIteration state variables.
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.