42#ifndef BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
43#define BELOS_PSEUDO_BLOCK_GMRES_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"
89 template <
class ScalarType,
class MV>
92 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
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>
148 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
166 Teuchos::ParameterList ¶ms );
246 for (
int i=0; i<
numRHS_; ++i) {
250 state.
sn[i] =
sn_[i];
251 state.
cs[i] =
cs_[i];
324 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
325 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
344 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
345 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
346 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
347 const Teuchos::RCP<OrthoManager<ScalarType,MV> >
ortho_;
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>
402 Teuchos::ParameterList ¶ms ):
414 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
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>
446 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
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) {
455 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
456 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
457 Teuchos::BLAS<int,ScalarType> blas;
459 for (
int i=0; i<numRHS_; ++i) {
461 Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
465 Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
469 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
470 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
471 H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
473 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
474 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
477 return currentUpdate;
484 template <
class ScalarType,
class MV,
class OP>
485 Teuchos::RCP<const MV>
489 typedef typename Teuchos::ScalarTraits<ScalarType> STS;
495 if (
static_cast<int> (norms->size()) < numRHS_)
496 norms->resize (numRHS_);
498 Teuchos::BLAS<int, ScalarType> blas;
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.");
529 TEUCHOS_TEST_FOR_EXCEPTION((
int)newstate.
V.size()==0 || (
int)newstate.
Z.size()==0,
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;
552 TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
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 "
562 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
curDim > numBlocks_+1,
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 );
607 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
608 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
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) {
622 Z_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>() );
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]) {
636 Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.
Z[i]->values(),curDim_+1);
637 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
638 lclZ = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
642 lclZ = Teuchos::null;
649 for (
int i=0; i<numRHS_; ++i) {
651 if (H_[i] == Teuchos::null) {
652 H_[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
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]) {
668 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.
H[i],curDim_+1, curDim_);
669 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
670 lclH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
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])
689 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.
cs[i]) );
690 if (sn_[i] != newstate.
sn[i])
691 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.
sn[i]) );
696 for (
int i=0; i<numRHS_; ++i) {
697 if (cs_[i] == Teuchos::null)
698 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
700 cs_[i]->resize(numBlocks_+1);
701 if (sn_[i] == Teuchos::null)
702 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
704 sn_[i]->resize(numBlocks_+1);
715 template <
class ScalarType,
class MV,
class OP>
721 if (initialized_ ==
false) {
725 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
726 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
729 int searchDim = numBlocks_;
734 std::vector<int> index(1);
735 std::vector<int> index2(1);
737 Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
740 Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
742 for (
int i=0; i<numRHS_; ++i) {
744 Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
745 Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
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) {
777 Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
778 Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
783 Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
787 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new
788 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
789 ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
790 Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
792 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
793 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
794 ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
799 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
804 index2[0] = curDim_+1;
805 Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
806 MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
812 Teuchos::RCP<MV> tmp_AU_vec = U_vec;
832 template<
class ScalarType,
class MV,
class OP>
838 int curDim = curDim_;
839 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
844 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
846 Teuchos::BLAS<int, ScalarType> blas;
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] );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which describes the basic interface to the 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.
Parent class to all Belos exceptions.
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...
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
PseudoBlockGmresIterOrthoFailure(const std::string &what_arg)
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
const Teuchos::RCP< OutputManager< ScalarType > > om_
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< MV > cur_block_rhs_
bool isInitialized()
States whether the solver has been initialized or not.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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 ...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the "native" residual vectors.
void setBlockSize(int blockSize)
Set the blocksize.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::RCP< MV > AU_vec_
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.
void resetNumIters(int iter=0)
Reset the iteration count.
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > H_
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, MagnitudeType > > > cs_
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > Z_
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
std::vector< Teuchos::RCP< MV > > V_
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
Teuchos::RCP< MV > U_vec_
Teuchos::RCP< MV > cur_block_sol_
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > sn_
int getNumIters() const
Get the current iteration count.
virtual ~PseudoBlockGmresIter()
Destructor.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > R_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > R
The current upper-triangular matrix from the QR reduction of H.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
PseudoBlockGmresIterState()
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
Teuchos::ScalarTraits< ScalarType > SCT
int curDim
The current dimension of the reduction.
SCT::magnitudeType MagnitudeType