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;
154 template<
class ScalarType,
class MV,
class OP>
164 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
182 Teuchos::ParameterList ¶ms );
257 state.
V.resize(numRHS_);
258 state.
H.resize(numRHS_);
259 state.
Z.resize(numRHS_);
260 state.
sn.resize(numRHS_);
261 state.
cs.resize(numRHS_);
262 for (
int i=0; i<numRHS_; ++i) {
266 state.
sn[i] = sn_[i];
267 state.
cs[i] = cs_[i];
301 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms )
const;
309 Teuchos::RCP<MV> getCurrentUpdate()
const;
315 void updateLSQR(
int dim = -1 );
319 if (!initialized_)
return 0;
340 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
341 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
348 void setNumBlocks(
int numBlocks);
360 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
361 const Teuchos::RCP<OutputManager<ScalarType> > om_;
362 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
363 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
374 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
375 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
378 Teuchos::RCP<MV> U_vec_, AU_vec_;
381 Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
397 std::vector<Teuchos::RCP<MV> > V_;
402 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
407 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
408 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
413 template<
class ScalarType,
class MV,
class OP>
418 Teuchos::ParameterList ¶ms ):
430 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
431 "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
432 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
439 template <
class ScalarType,
class MV,
class OP>
445 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
447 numBlocks_ = numBlocks;
450 initialized_ =
false;
455 template <
class ScalarType,
class MV,
class OP>
462 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
464 return currentUpdate;
466 currentUpdate =
MVT::Clone(*(V_[0]), numRHS_);
467 std::vector<int> index(1), index2(curDim_);
468 for (
int i=0; i<curDim_; ++i) {
471 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
472 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
473 Teuchos::BLAS<int,ScalarType> blas;
475 for (
int i=0; i<numRHS_; ++i) {
481 Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
485 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
486 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
487 H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
493 return currentUpdate;
500 template <
class ScalarType,
class MV,
class OP>
501 Teuchos::RCP<const MV>
505 typedef typename Teuchos::ScalarTraits<ScalarType> STS;
511 if (static_cast<int> (norms->size()) < numRHS_)
512 norms->resize (numRHS_);
514 Teuchos::BLAS<int, ScalarType> blas;
515 for (
int j = 0; j < numRHS_; ++j)
517 const ScalarType curNativeResid = (*Z_[j])(curDim_);
518 (*norms)[j] = STS::magnitude (curNativeResid);
521 return Teuchos::null;
525 template <
class ScalarType,
class MV,
class OP>
540 std::string errstr (
"Belos::PseudoBlockGmresIter::initialize(): " 541 "Specified multivectors must have a consistent " 542 "length and width.");
545 TEUCHOS_TEST_FOR_EXCEPTION((
int)newstate.
V.size()==0 || (int)newstate.
Z.size()==0,
546 std::invalid_argument,
547 "Belos::PseudoBlockGmresIter::initialize(): " 548 "V and/or Z was not specified in the input state; " 549 "the V and/or Z arrays have length zero.");
560 RCP<const MV> lhsMV = lp_->getLHS();
561 RCP<const MV> rhsMV = lp_->getRHS();
565 RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
568 TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
569 std::invalid_argument,
570 "Belos::PseudoBlockGmresIter::initialize(): " 571 "The linear problem to solve does not specify multi" 572 "vectors from which we can clone basis vectors. The " 573 "right-hand side(s), left-hand side(s), or both should " 578 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
curDim > numBlocks_+1,
579 std::invalid_argument,
581 curDim_ = newstate.
curDim;
587 for (
int i=0; i<numRHS_; ++i) {
592 V_[i] =
MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
597 std::invalid_argument, errstr );
599 std::invalid_argument, errstr );
605 if (newstate.
V[i] != V_[i]) {
607 if (curDim_ == 0 && lclDim > 1) {
609 <<
"Belos::PseudoBlockGmresIter::initialize(): the solver was " 610 <<
"initialized with a kernel of " << lclDim
612 <<
"The block size however is only " << 1
614 <<
"The last " << lclDim - 1 <<
" vectors will be discarded." 617 std::vector<int> nevind (curDim_ + 1);
618 for (
int j = 0; j < curDim_ + 1; ++j)
623 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
624 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
628 lclV = Teuchos::null;
635 for (
int i=0; i<numRHS_; ++i) {
637 if (Z_[i] == Teuchos::null) {
638 Z_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>() );
640 if (Z_[i]->length() < numBlocks_+1) {
641 Z_[i]->shapeUninitialized(numBlocks_+1, 1);
645 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
648 if (newstate.
Z[i] != Z_[i]) {
652 Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.
Z[i]->values(),curDim_+1);
653 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
654 lclZ = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
658 lclZ = Teuchos::null;
665 for (
int i=0; i<numRHS_; ++i) {
667 if (H_[i] == Teuchos::null) {
668 H_[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
670 if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
671 H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
675 if ((
int)newstate.
H.size() == numRHS_) {
678 TEUCHOS_TEST_FOR_EXCEPTION((newstate.
H[i]->numRows() < curDim_ || newstate.
H[i]->numCols() < curDim_), std::invalid_argument,
679 "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
681 if (newstate.
H[i] != H_[i]) {
684 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.
H[i],curDim_+1, curDim_);
685 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
686 lclH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
690 lclH = Teuchos::null;
702 if ((
int)newstate.
cs.size() == numRHS_ && (int)newstate.
sn.size() == numRHS_) {
703 for (
int i=0; i<numRHS_; ++i) {
704 if (cs_[i] != newstate.
cs[i])
705 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.
cs[i]) );
706 if (sn_[i] != newstate.
sn[i])
707 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.
sn[i]) );
712 for (
int i=0; i<numRHS_; ++i) {
713 if (cs_[i] == Teuchos::null)
714 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
716 cs_[i]->resize(numBlocks_+1);
717 if (sn_[i] == Teuchos::null)
718 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
720 sn_[i]->resize(numBlocks_+1);
742 template <
class ScalarType,
class MV,
class OP>
748 if (initialized_ ==
false) {
752 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
753 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
756 int searchDim = numBlocks_;
761 std::vector<int> index(1);
762 std::vector<int> index2(1);
764 Teuchos::RCP<MV> U_vec =
MVT::Clone( *V_[0], numRHS_ );
767 Teuchos::RCP<MV> AU_vec =
MVT::Clone( *V_[0], numRHS_ );
769 for (
int i=0; i<numRHS_; ++i) {
773 MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
781 while (stest_->checkStatus(
this) !=
Passed && curDim_ < searchDim) {
787 lp_->apply( *U_vec, *AU_vec );
792 int num_prev = curDim_+1;
793 index.resize( num_prev );
794 for (
int i=0; i<num_prev; ++i) {
800 for (
int i=0; i<numRHS_; ++i) {
805 Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
814 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new
815 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
816 ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
817 Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
819 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
820 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
821 ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
826 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
831 index2[0] = curDim_+1;
839 Teuchos::RCP<MV> tmp_AU_vec = U_vec;
875 template<
class ScalarType,
class MV,
class OP>
881 int curDim = curDim_;
887 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
889 Teuchos::BLAS<int, ScalarType> blas;
891 for (i=0; i<numRHS_; ++i) {
897 for (j=0; j<curDim; j++) {
901 blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
906 blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
907 (*H_[i])(curDim+1,curDim) = zero;
911 blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
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.
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.
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
Pure virtual base class for defining the status testing capabilities of Belos.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
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.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
PseudoBlockGmresIterInitFailure(const std::string &what_arg)
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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 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 ...
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
PseudoBlockGmresIterInitFailure is thrown when the PseudoBlockGmresIter object is unable to generate ...