42 #ifndef BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
43 #define BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
89 template <
class ScalarType,
class MV>
101 std::vector<Teuchos::RCP<const MV> >
V;
107 std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > >
H;
109 std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > >
R;
111 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > >
Z;
113 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > >
sn;
114 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > >
cs;
138 template<
class ScalarType,
class MV,
class OP>
241 state.
V.resize(numRHS_);
242 state.
H.resize(numRHS_);
243 state.
Z.resize(numRHS_);
244 state.
sn.resize(numRHS_);
245 state.
cs.resize(numRHS_);
246 for (
int i=0; i<numRHS_; ++i) {
250 state.
sn[i] = sn_[i];
251 state.
cs[i] = cs_[i];
303 if (!initialized_)
return 0;
325 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
358 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
359 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
381 std::vector<Teuchos::RCP<MV> > V_;
386 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
391 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
392 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
397 template<
class ScalarType,
class MV,
class OP>
415 "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
416 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
423 template <
class ScalarType,
class MV,
class OP>
429 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
431 numBlocks_ = numBlocks;
434 initialized_ =
false;
439 template <
class ScalarType,
class MV,
class OP>
448 return currentUpdate;
450 currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
451 std::vector<int> index(1), index2(curDim_);
452 for (
int i=0; i<curDim_; ++i) {
459 for (
int i=0; i<numRHS_; ++i) {
461 Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
471 H_[i]->values(), H_[i]->stride(), y.
values(), y.
stride() );
474 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
477 return currentUpdate;
484 template <
class ScalarType,
class MV,
class OP>
495 if (static_cast<int> (norms->size()) < numRHS_)
496 norms->resize (numRHS_);
499 for (
int j = 0; j < numRHS_; ++j)
501 const ScalarType curNativeResid = (*Z_[j])(curDim_);
502 (*norms)[j] = STS::magnitude (curNativeResid);
505 return Teuchos::null;
509 template <
class ScalarType,
class MV,
class OP>
518 this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
524 std::string errstr (
"Belos::PseudoBlockGmresIter::initialize(): "
525 "Specified multivectors must have a consistent "
526 "length and width.");
530 std::invalid_argument,
531 "Belos::PseudoBlockGmresIter::initialize(): "
532 "V and/or Z was not specified in the input state; "
533 "the V and/or Z arrays have length zero.");
544 RCP<const MV> lhsMV = lp_->getLHS();
545 RCP<const MV> rhsMV = lp_->getRHS();
549 RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
553 std::invalid_argument,
554 "Belos::PseudoBlockGmresIter::initialize(): "
555 "The linear problem to solve does not specify multi"
556 "vectors from which we can clone basis vectors. The "
557 "right-hand side(s), left-hand side(s), or both should "
563 std::invalid_argument,
565 curDim_ = newstate.
curDim;
571 for (
int i=0; i<numRHS_; ++i) {
575 if (V_[i].
is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
576 V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
580 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
V[i]) != MVT::GetGlobalLength(*V_[i]),
581 std::invalid_argument, errstr );
582 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
V[i]) < newstate.
curDim,
583 std::invalid_argument, errstr );
588 int lclDim = MVT::GetNumberVecs(*newstate.
V[i]);
589 if (newstate.
V[i] != V_[i]) {
591 if (curDim_ == 0 && lclDim > 1) {
593 <<
"Belos::PseudoBlockGmresIter::initialize(): the solver was "
594 <<
"initialized with a kernel of " << lclDim
596 <<
"The block size however is only " << 1
598 <<
"The last " << lclDim - 1 <<
" vectors will be discarded."
601 std::vector<int> nevind (curDim_ + 1);
602 for (
int j = 0; j < curDim_ + 1; ++j)
605 RCP<const MV> newV = MVT::CloneView (*newstate.
V[i], nevind);
606 RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
609 MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
612 lclV = Teuchos::null;
619 for (
int i=0; i<numRHS_; ++i) {
621 if (Z_[i] == Teuchos::null) {
624 if (Z_[i]->length() < numBlocks_+1) {
625 Z_[i]->shapeUninitialized(numBlocks_+1, 1);
629 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
632 if (newstate.
Z[i] != Z_[i]) {
642 lclZ = Teuchos::null;
649 for (
int i=0; i<numRHS_; ++i) {
651 if (H_[i] == Teuchos::null) {
654 if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
655 H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
659 if ((
int)newstate.
H.size() == numRHS_) {
662 TEUCHOS_TEST_FOR_EXCEPTION((newstate.
H[i]->numRows() < curDim_ || newstate.
H[i]->numCols() < curDim_), std::invalid_argument,
663 "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
665 if (newstate.
H[i] != H_[i]) {
674 lclH = Teuchos::null;
686 if ((
int)newstate.
cs.size() == numRHS_ && (int)newstate.
sn.size() == numRHS_) {
687 for (
int i=0; i<numRHS_; ++i) {
688 if (cs_[i] != newstate.
cs[i])
690 if (sn_[i] != newstate.
sn[i])
696 for (
int i=0; i<numRHS_; ++i) {
697 if (cs_[i] == Teuchos::null)
700 cs_[i]->resize(numBlocks_+1);
701 if (sn_[i] == Teuchos::null)
704 sn_[i]->resize(numBlocks_+1);
715 template <
class ScalarType,
class MV,
class OP>
721 if (initialized_ ==
false) {
729 int searchDim = numBlocks_;
734 std::vector<int> index(1);
735 std::vector<int> index2(1);
742 for (
int i=0; i<numRHS_; ++i) {
746 MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
754 while (stest_->checkStatus(
this) !=
Passed && curDim_ < searchDim) {
760 lp_->apply( *U_vec, *AU_vec );
765 int num_prev = curDim_+1;
766 index.resize( num_prev );
767 for (
int i=0; i<num_prev; ++i) {
773 for (
int i=0; i<numRHS_; ++i) {
799 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
804 index2[0] = curDim_+1;
806 MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
832 template<
class ScalarType,
class MV,
class OP>
838 int curDim = curDim_;
839 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
848 for (i=0; i<numRHS_; ++i) {
854 for (j=0; j<curDim; j++) {
858 blas.
ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
863 blas.
ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
864 (*H_[i])(curDim+1,curDim) = zero;
868 blas.
ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
ScalarType * values() const
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
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...
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > R
The current upper-triangular matrix from the QR reduction of H.
Teuchos::ScalarTraits< ScalarType > SCT
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
virtual ~PseudoBlockGmresIter()
Destructor.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
int resize(OrdinalType length_in)
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the "native" residual vectors.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
Structure to contain pointers to PseudoBlockGmresIter state variables.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
bool isInitialized()
States whether the solver has been initialized or not.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
PseudoBlockGmresIter(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)
PseudoBlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver ...
PseudoBlockGmresIterState()
PseudoBlockGmresIterOrthoFailure(const std::string &what_arg)
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void resetNumIters(int iter=0)
Reset the iteration count.
MultiVecTraits< ScalarType, MV > MVT
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
OperatorTraits< ScalarType, MV, OP > OPT
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
Class which defines basic traits for the operator type.
SCT::magnitudeType MagnitudeType
int curDim
The current dimension of the reduction.
Parent class to all Belos exceptions.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
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...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
OrdinalType stride() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...