42 #ifndef BELOS_GMRESPOLYOP_HPP
43 #define BELOS_GMRESPOLYOP_HPP
79 #ifdef BELOS_TEUCHOS_TIME_MONITOR
81 #endif // BELOS_TEUCHOS_TIME_MONITOR
96 template <
class ScalarType,
class MV>
106 mv_ = Teuchos::rcp_const_cast<MV>( mv_in );
194 template <
class ScalarType,
class MV,
class OP>
207 LP_(problem_in->getLeftPrec()),
208 RP_(problem_in->getRightPrec())
213 #ifdef BELOS_TEUCHOS_TIME_MONITOR
215 #endif // BELOS_TEUCHOS_TIME_MONITOR
223 "Belos::GmresPolyOp: \"Polynomial Type\" must be either \"Arnoldi\", \"Gmres\", or \"Roots\".");
268 void ApplyPoly (
const MV& x, MV& y )
const;
287 #ifdef BELOS_TEUCHOS_TIME_MONITOR
289 #endif // BELOS_TEUCHOS_TIME_MONITOR
352 template <
class ScalarType,
class MV,
class OP>
357 polyType_ = params_in->
get(
"Polynomial Type", polyType_default_);
361 if (params_in->
isParameter(
"Polynomial Tolerance")) {
363 polyTol_ = params_in->
get (
"Polynomial Tolerance",
373 maxDegree_ = params_in->
get(
"Maximum Degree", maxDegree_default_);
378 randomRHS_ = params_in->
get(
"Random RHS", randomRHS_default_);
383 if (Teuchos::isParameterType<int>(*params_in,
"Verbosity")) {
384 verbosity_ = params_in->
get(
"Verbosity", verbosity_default_);
387 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,
"Verbosity");
392 orthoType_ = params_in->
get(
"Orthogonalization",orthoType_default_);
397 label_ = params_in->
get(
"Timer Label", label_default_);
402 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,
"Output Stream");
407 damp_ = params_in->
get(
"Damped Poly", damp_default_);
412 addRoots_ = params_in->
get(
"Add Roots", addRoots_default_);
416 template <
class ScalarType,
class MV,
class OP>
422 std::vector<int> index(1,0);
425 MVT::MvRandom( *V0 );
427 MVT::Assign( *problem_->getRHS(), *V0 );
429 if ( !LP_.is_null() ) {
431 problem_->applyLeftPrec( *Vtemp, *V0);
435 problem_->apply( *Vtemp, *V0);
438 for(
int i=0; i< maxDegree_; i++)
444 problem_->apply( *Vi, *Vip1);
453 MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
462 while( status && dim_ >= 1)
472 std::cout <<
"BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
473 std::cout <<
"Error code: " << infoInt << std::endl;
494 pCoeff_.
shape( 1, 1);
495 pCoeff_(0,0) = SCT::one();
496 std::cout <<
"Poly Degree is zero. No preconditioner created." << std::endl;
500 pCoeff_.shape( dim_, 1);
506 MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
507 lapack.
POTRS(
'U', dim_, 1, lhs.
values(), lhs.
stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
510 std::cout <<
"BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
511 std::cout <<
"Error code: " << infoInt << std::endl;
516 template <
class ScalarType,
class MV,
class OP>
519 std::string polyLabel = label_ +
": GmresPolyOp creation";
522 std::vector<int> idx(1,0);
525 MVT::MvInit( *newX, SCT::zero() );
527 MVT::MvRandom( *newB );
530 MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
534 newProblem->setInitResVec( newB );
535 newProblem->setLeftPrec( problem_->getLeftPrec() );
536 newProblem->setRightPrec( problem_->getRightPrec() );
537 newProblem->setLabel(polyLabel);
538 newProblem->setProblem();
539 newProblem->setLSIndex( idx );
545 polyList.
set(
"Num Blocks",maxDegree_);
546 polyList.
set(
"Block Size",1);
547 polyList.
set(
"Keep Hessenberg",
true);
553 if (ortho_.is_null()) {
554 params_->set(
"Orthogonalization", orthoType_);
580 if ( !LP_.is_null() )
581 newProblem->applyLeftPrec( *newB, *V_0 );
585 newProblem->apply( *Vtemp, *V_0 );
592 int rank = ortho_->normalize( *V_0, Teuchos::rcpFromRef(r0_) );
594 "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
599 newstate.
z = Teuchos::rcpFromRef( r0_);
601 gmres_iter->initializeGmres(newstate);
605 gmres_iter->iterate();
609 gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
611 catch (std::exception& e) {
613 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGmresIter::iterate() at iteration "
614 << gmres_iter->getNumIters() << endl << e.what () << endl;
629 if(polyType_ ==
"Arnoldi"){
642 y_.values(), y_.stride() );
650 for(
int i=0; i <= dim_-3; i++) {
651 for(
int k=i+2; k <= dim_-1; k++) {
652 H_(k,i) = SCT::zero();
659 ScalarType Hlast = (*gmresState.
H)(dim_,dim_-1);
665 E(dim_-1,0) = SCT::one();
668 HSolver.
setMatrix( Teuchos::rcpFromRef(Htemp));
670 HSolver.
setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
677 std::cout <<
"Hsolver factor: info = " << info << std::endl;
679 info = HSolver.
solve();
681 std::cout <<
"Hsolver solve : info = " << info << std::endl;
685 F.scale(Hlast*Hlast);
690 theta_.shape(dim_,2);
697 std::vector<ScalarType> work(1);
698 std::vector<MagnitudeType> rwork(2*dim_);
701 lapack.
GEEV(
'N',
'N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
703 work.resize( lwork );
705 lapack.
GEEV(
'N',
'N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
708 std::cout <<
"GEEV solve : info = " << info << std::endl;
714 std::vector<int> index(dim_);
715 for(
int i=0; i<dim_; ++i){
718 TEUCHOS_TEST_FOR_EXCEPTION(hypot(theta_(i,0),theta_(i,1)) < tol, std::runtime_error,
"BelosGmresPolyOp Error: One of the computed polynomial roots is approximately zero. This will cause a divide by zero error! Your matrix may be close to singular. Please select a lower polynomial degree or give a shifted matrix.");
720 SortModLeja(theta_,index);
729 template <
class ScalarType,
class MV,
class OP>
734 std::vector<std::complex<MagnitudeType>> cmplxHRitz (dim_);
735 for(
unsigned int i=0; i<cmplxHRitz.size(); ++i){
736 cmplxHRitz[i] = std::complex<MagnitudeType>( theta_(i,0), theta_(i,1) );
741 std::vector<MagnitudeType> pof (dim_,one);
742 for(
int j=0; j<dim_; ++j) {
743 for(
int i=0; i<dim_; ++i) {
745 pof[j] = std::abs(pof[j]*(one-(cmplxHRitz[j]/cmplxHRitz[i])));
751 std::vector<int> extra (dim_);
753 for(
int i=0; i<dim_; ++i){
754 if (pof[i] > MCT::zero())
759 totalExtra += extra[i];
763 printer_->stream(
Warnings) <<
"Warning: Need to add " << totalExtra <<
" extra roots." << std::endl;}
766 if(addRoots_ && totalExtra>0)
768 theta_.reshape(dim_+totalExtra,2);
774 for(
int i=0; i<dim_; ++i){
775 for(
int j=0; j< extra[i]; ++j){
776 theta_(count,0) = theta_(i,0);
777 theta_(count,1) = theta_(i,1);
778 thetaPert(count,0) = theta_(i,0)+(j+MCT::one())*
MagnitudeType(5e-8);
779 thetaPert(count,1) = theta_(i,1);
787 printer_->stream(
Warnings) <<
"New poly degree is: " << dim_ << std::endl;}
790 std::vector<int> index2(dim_);
791 for(
int i=0; i<dim_; ++i){
794 SortModLeja(thetaPert,index2);
796 for(
int i=0; i<dim_; ++i)
798 thetaPert(i,0) = theta_(index2[i],0);
799 thetaPert(i,1) = theta_(index2[i],1);
808 template <
class ScalarType,
class MV,
class OP>
814 int dimN = index.size();
815 std::vector<int> newIndex(dimN);
821 for(
int i = 0; i < dimN; i++){
822 absVal(i) = hypot(thetaN(i,0), thetaN(i,1));
824 MagnitudeType * maxPointer = std::max_element(absVal.values(), (absVal.values()+dimN));
825 int maxIndex = int (maxPointer- absVal.values());
828 sorted(0,0) = thetaN(maxIndex,0);
829 sorted(0,1) = thetaN(maxIndex,1);
830 newIndex[0] = index[maxIndex];
834 if(sorted(0,1)!= SCT::zero() && !SCT::isComplex)
836 sorted(1,0) = thetaN(maxIndex,0);
837 sorted(1,1) = -thetaN(maxIndex,1);
838 newIndex[1] = index[maxIndex+1];
851 for(
int i = 0; i < dimN; i++)
853 prod(i) = MCT::one();
854 for(
int k = 0; k < j; k++)
856 a = thetaN(i,0) - sorted(k,0);
857 b = thetaN(i,1) - sorted(k,1);
858 if (a*a + b*b > MCT::zero())
859 prod(i) = prod(i) + log10(hypot(a,b));
861 prod(i) = -std::numeric_limits<MagnitudeType>::infinity();
868 maxPointer = std::max_element(prod.values(), (prod.values()+dimN));
869 maxIndex = int (maxPointer- prod.values());
870 sorted(j,0) = thetaN(maxIndex,0);
871 sorted(j,1) = thetaN(maxIndex,1);
872 newIndex[j] = index[maxIndex];
875 if(sorted(j,1)!= SCT::zero() && !SCT::isComplex)
878 sorted(j,0) = thetaN(maxIndex,0);
879 sorted(j,1) = -thetaN(maxIndex,1);
880 newIndex[j] = index[maxIndex+1];
890 template <
class ScalarType,
class MV,
class OP>
894 if (polyType_ ==
"Arnoldi")
895 ApplyArnoldiPoly(x, y);
896 else if (polyType_ ==
"Gmres")
897 ApplyGmresPoly(x, y);
898 else if (polyType_ ==
"Roots")
899 ApplyRootsPoly(x, y);
903 problem_->applyOp( x, y );
907 template <
class ScalarType,
class MV,
class OP>
914 if (!LP_.is_null()) {
916 problem_->applyLeftPrec( *AX, *Xtmp );
921 #ifdef BELOS_TEUCHOS_TIME_MONITOR
924 MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y);
926 for(
int i=1; i < dim_; i++)
939 problem_->apply(*X, *Y);
941 #ifdef BELOS_TEUCHOS_TIME_MONITOR
944 MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y);
949 if (!RP_.is_null()) {
951 problem_->applyRightPrec( *Ytmp, y );
955 template <
class ScalarType,
class MV,
class OP>
958 MVT::MvInit( y, SCT::zero() );
964 if (!LP_.is_null()) {
965 problem_->applyLeftPrec( *prod, *Xtmp );
972 if(theta_(i,1)== SCT::zero() || SCT::isComplex)
975 #ifdef BELOS_TEUCHOS_TIME_MONITOR
978 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y);
980 problem_->apply(*prod, *Xtmp);
982 #ifdef BELOS_TEUCHOS_TIME_MONITOR
985 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/theta_(i,0), *Xtmp, *prod);
991 MagnitudeType mod = theta_(i,0)*theta_(i,0) + theta_(i,1)*theta_(i,1);
992 problem_->apply(*prod, *Xtmp);
994 #ifdef BELOS_TEUCHOS_TIME_MONITOR
997 MVT::MvAddMv(2*theta_(i,0), *prod, -SCT::one(), *Xtmp, *Xtmp);
998 MVT::MvAddMv(SCT::one(), y, SCT::one()/mod, *Xtmp, y);
1002 problem_->apply(*Xtmp, *Xtmp2);
1004 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1007 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/mod, *Xtmp2, *prod);
1013 if(theta_(dim_-1,1)== SCT::zero() || SCT::isComplex)
1015 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1018 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y);
1022 if (!RP_.is_null()) {
1024 problem_->applyRightPrec( *Ytmp, y );
1028 template <
class ScalarType,
class MV,
class OP>
1033 V_ = MVT::Clone( x, dim_ );
1034 if (!LP_.is_null()) {
1035 wL_ = MVT::Clone( y, 1 );
1037 if (!RP_.is_null()) {
1038 wR_ = MVT::Clone( y, 1 );
1044 int n = MVT::GetNumberVecs( x );
1045 std::vector<int> idxi(1), idxi2, idxj(1);
1048 for (
int j=0; j<n; ++j) {
1054 if (!LP_.is_null()) {
1056 problem_->applyLeftPrec( *x_view, *v_curr );
1058 MVT::SetBlock( *x_view, idxi, *V_ );
1061 for (
int i=0; i<dim_-1; ++i) {
1065 for (
int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
1077 if (!RP_.is_null()) {
1078 problem_->applyRightPrec( *v_curr, *wR_ );
1083 if (LP_.is_null()) {
1087 problem_->applyOp( *wR_, *wL_ );
1089 if (!LP_.is_null()) {
1090 problem_->applyLeftPrec( *wL_, *v_next );
1096 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1099 MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
1103 MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
1107 if (!RP_.is_null()) {
1109 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1112 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1114 problem_->applyRightPrec( *wR_, *y_view );
1117 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1120 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
ScalarType * values() const
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
Collection of types and exceptions used within the Belos solvers.
const GmresPolyMv * CloneView(const std::vector< int > &index) const
Creates a new Belos::MultiVec that shares the selected contents of *this. The index of the numvecs ve...
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...
static constexpr bool randomRHS_default_
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
Class which manages the output and verbosity of the Belos solvers.
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
void MvAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const ScalarType beta, const MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
Teuchos::RCP< OutputManager< ScalarType > > printer_
static magnitudeType eps()
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
static constexpr const char * label_default_
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< MV > getMV()
T & get(ParameterList &l, const std::string &name)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void ApplyGmresPoly(const MV &x, MV &y) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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).
void generateGmresPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
MultiVecTraits< ScalarType, MV > MVT
Declaration of basic traits for the multivector type.
void POTRS(const char &UPLO, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< const OP > RP_
GmresPolyMv * CloneViewNonConst(const std::vector< int > &index)
Creates a new Belos::MultiVec that shares the selected contents of *this. The index of the numvecs ve...
GmresPolyMv * CloneCopy(const std::vector< int > &index) const
Creates a new Belos::MultiVec and copies the selected contents of *this into the new multivector (dee...
MultiVecTraits< ScalarType, MV > MVT
void MvRandom()
Fill all the vectors in *this with random numbers.
Structure to contain pointers to GmresIteration state variables.
static constexpr int verbosity_default_
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in, const Teuchos::RCP< Teuchos::ParameterList > ¶ms_in)
Basic contstructor.
Belos::StatusTest class for specifying a maximum number of iterations.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static const double polyTol
Relative residual tolerance for matrix polynomial construction.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
Class which defines basic traits for the operator type.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
A factory class for generating StatusTestOutput objects.
static constexpr const char * polyType_default_
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Teuchos::ScalarTraits< MagnitudeType > MCT
ETrans
Whether to apply the (conjugate) transpose of an operator.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
A Belos::StatusTest class for specifying a maximum number of iterations.
static constexpr int maxDegree_default_
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).
void MvTimesMatAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta)
Update *this with alpha * A * B + beta * (*this).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Alternative run-time polymorphic interface for operators.
Teuchos::SerialDenseMatrix< OT, ScalarType > y_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int curDim
The current dimension of the reduction.
void ApplyPoly(const MV &x, MV &y) const
This routine takes the MV x and applies the polynomial operator phi(OP) to it resulting in the MV y...
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
void factorWithEquilibration(bool flag)
void ApplyArnoldiPoly(const MV &x, MV &y) const
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
void solveWithTransposeFlag(Teuchos::ETransp trans)
A linear system to solve, and its associated information.
Teuchos::RCP< const OP > LP_
Class which describes the linear problem to be solved by the iterative solver.
GmresPolyOpOrthoFailure(const std::string &what_arg)
void SetBlock(const MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
void ApplyRootsPoly(const MV &x, MV &y) const
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_
Teuchos::RCP< const MV > getConstMV() const
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_
Belos concrete class for performing the block GMRES iteration.
std::string polyUpdateLabel_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
OrdinalType numCols() const
Teuchos::ScalarTraits< ScalarType > SCT
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector< int > &index) const
NormType
The type of vector norm to compute.
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
void MvDot(const MultiVec< ScalarType > &A, std::vector< ScalarType > &b) const
Compute the dot product of each column of *this with the corresponding column of A.
void MvTransMv(const ScalarType alpha, const MultiVec< ScalarType > &A, Teuchos::SerialDenseMatrix< int, ScalarType > &B) const
Compute a dense matrix B through the matrix-matrix multiply alpha * A^T * (*this).
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Alternative run-time polymorphic interface for operators.
virtual ~GmresPolyOp()
Destructor.
void GEEV(const char &JOBVL, const char &JOBVR, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, MagnitudeType *WR, MagnitudeType *WI, ScalarType *VL, const OrdinalType &ldvl, ScalarType *VR, const OrdinalType &ldvr, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) const
Teuchos::RCP< Teuchos::ParameterList > params_
static constexpr bool damp_default_
GmresPolyMv * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
bool isType(const std::string &name) const
A class for extending the status testing capabilities of Belos via logical combinations.
Interface for multivectors used by Belos' linear solvers.
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in)
Given no ParameterList, constructor creates no polynomial and only applies the given operator...
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
void generateArnoldiPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Parent class to all Belos exceptions.
int shape(OrdinalType numRows, OrdinalType numCols)
static constexpr bool addRoots_default_
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...
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static constexpr const char * orthoType_default_
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
Teuchos::SerialDenseMatrix< OT, ScalarType > H_
OrdinalType stride() const
void Apply(const MultiVec< ScalarType > &x, MultiVec< ScalarType > &y, ETrans=NOTRANS) const
This routine casts the MultiVec to GmresPolyMv to retrieve the MV. Then the above apply method is cal...
Teuchos::RCP< std::ostream > outputStream_
Teuchos::SerialDenseVector< OT, ScalarType > r0_
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms_in)
Process the passed in parameters.
OrdinalType numRows() const
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Interface for multivectors used by Belos' linear solvers.
GmresPolyMv(const Teuchos::RCP< MV > &mv_in)
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.