42#ifndef BELOS_GMRESPOLYOP_HPP
43#define BELOS_GMRESPOLYOP_HPP
70#include "Teuchos_BLAS.hpp"
71#include "Teuchos_LAPACK.hpp"
72#include "Teuchos_as.hpp"
73#include "Teuchos_RCP.hpp"
74#include "Teuchos_SerialDenseMatrix.hpp"
75#include "Teuchos_SerialDenseVector.hpp"
76#include "Teuchos_SerialDenseSolver.hpp"
77#include "Teuchos_ParameterList.hpp"
79#ifdef BELOS_TEUCHOS_TIME_MONITOR
80 #include "Teuchos_TimeMonitor.hpp"
96 template <
class ScalarType,
class MV>
106 mv_ = Teuchos::rcp_const_cast<MV>( mv_in );
108 Teuchos::RCP<MV>
getMV() {
return mv_; }
140 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
const ScalarType beta)
163 void MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec,
NormType type =
TwoNorm )
const
180 Teuchos::RCP<MV> mv_;
194 template <
class ScalarType,
class MV,
class OP>
203 const Teuchos::RCP<Teuchos::ParameterList>& params_in
205 : problem_(problem_in),
207 LP_(problem_in->getLeftPrec()),
208 RP_(problem_in->getRightPrec())
212 polyUpdateLabel_ = label_ +
": Hybrid Gmres: Vector Update";
213#ifdef BELOS_TEUCHOS_TIME_MONITOR
214 timerPolyUpdate_ = Teuchos::TimeMonitor::getNewCounter(polyUpdateLabel_);
217 if (polyType_ ==
"Arnoldi" || polyType_==
"Roots")
219 else if (polyType_ ==
"Gmres")
222 TEUCHOS_TEST_FOR_EXCEPTION(polyType_!=
"Arnoldi"&&polyType_!=
"Gmres"&&polyType_!=
"Roots",std::invalid_argument,
223 "Belos::GmresPolyOp: \"Polynomial Type\" must be either \"Arnoldi\", \"Gmres\", or \"Roots\".");
228 : problem_(problem_in)
242 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList>& params_in );
268 void ApplyPoly (
const MV& x, MV& y )
const;
287#ifdef BELOS_TEUCHOS_TIME_MONITOR
288 Teuchos::RCP<Teuchos::Time> timerPolyUpdate_;
290 std::string polyUpdateLabel_;
294 typedef Teuchos::ScalarTraits<ScalarType> SCT ;
295 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
296 typedef Teuchos::ScalarTraits<MagnitudeType> MCT ;
299 static constexpr int maxDegree_default_ = 25;
301 static constexpr bool randomRHS_default_ =
true;
302 static constexpr const char * label_default_ =
"Belos";
303 static constexpr const char * polyType_default_ =
"Roots";
304 static constexpr const char * orthoType_default_ =
"DGKS";
305 static constexpr bool damp_default_ =
false;
306 static constexpr bool addRoots_default_ =
true;
309 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
310 Teuchos::RCP<Teuchos::ParameterList> params_;
311 Teuchos::RCP<const OP> LP_, RP_;
314 Teuchos::RCP<OutputManager<ScalarType> > printer_;
315 Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcpFromRef(std::cout);
318 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
322 int maxDegree_ = maxDegree_default_;
323 int verbosity_ = verbosity_default_;
324 bool randomRHS_ = randomRHS_default_;
325 std::string label_ = label_default_;
326 std::string polyType_ = polyType_default_;
327 std::string orthoType_ = orthoType_default_;
329 bool damp_ = damp_default_;
330 bool addRoots_ = addRoots_default_;
333 mutable Teuchos::RCP<MV> V_, wL_, wR_;
334 Teuchos::SerialDenseMatrix<OT,ScalarType> H_, y_;
335 Teuchos::SerialDenseVector<OT,ScalarType> r0_;
338 bool autoDeg =
false;
339 Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_;
342 Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_;
346 void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index)
const ;
349 void ComputeAddedRoots();
352 template <
class ScalarType,
class MV,
class OP>
356 if (params_in->isParameter(
"Polynomial Type")) {
357 polyType_ = params_in->get(
"Polynomial Type", polyType_default_);
361 if (params_in->isParameter(
"Polynomial Tolerance")) {
362 if (params_in->isType<MagnitudeType> (
"Polynomial Tolerance")) {
363 polyTol_ = params_in->get (
"Polynomial Tolerance",
372 if (params_in->isParameter(
"Maximum Degree")) {
373 maxDegree_ = params_in->get(
"Maximum Degree", maxDegree_default_);
377 if (params_in->isParameter(
"Random RHS")) {
378 randomRHS_ = params_in->get(
"Random RHS", randomRHS_default_);
382 if (params_in->isParameter(
"Verbosity")) {
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");
391 if (params_in->isParameter(
"Orthogonalization")) {
392 orthoType_ = params_in->get(
"Orthogonalization",orthoType_default_);
396 if (params_in->isParameter(
"Timer Label")) {
397 label_ = params_in->get(
"Timer Label", label_default_);
401 if (params_in->isParameter(
"Output Stream")) {
402 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,
"Output Stream");
406 if (params_in->isParameter(
"Damped Poly")) {
407 damp_ = params_in->get(
"Damped Poly", damp_default_);
411 if (params_in->isParameter(
"Add Roots")) {
412 addRoots_ = params_in->get(
"Add Roots", addRoots_default_);
416 template <
class ScalarType,
class MV,
class OP>
419 Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+1 );
422 std::vector<int> index(1,0);
423 Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
425 MVT::MvRandom( *V0 );
427 MVT::Assign( *problem_->getRHS(), *V0 );
429 if ( !LP_.is_null() ) {
430 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
431 problem_->applyLeftPrec( *Vtemp, *V0);
434 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
435 problem_->apply( *Vtemp, *V0);
438 for(
int i=0; i< maxDegree_; i++)
441 Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
443 Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
444 problem_->apply( *Vi, *Vip1);
448 Teuchos::Range1D range( 1, maxDegree_);
449 Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
452 Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_, maxDegree_);
453 MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
456 Teuchos::LAPACK< OT, ScalarType > lapack;
461 Teuchos::SerialDenseMatrix< OT, ScalarType > lhs;
462 while( status && dim_ >= 1)
464 Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_, dim_);
465 lapack.POTRF(
'U', dim_, lhstemp.values(), lhstemp.stride(), &infoInt);
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);
502 Teuchos::Range1D rangeSub( 1, dim_);
503 Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
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);
523 Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
524 Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
525 MVT::MvInit( *newX, SCT::zero() );
527 MVT::MvRandom( *newB );
530 MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
532 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
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 );
542 Teuchos::ParameterList polyList;
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_);
556 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
558 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
562 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
566 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
571 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
575 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
579 Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
580 if ( !LP_.is_null() )
581 newProblem->applyLeftPrec( *newB, *V_0 );
584 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*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;
619 Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
629 if(polyType_ ==
"Arnoldi"){
632 y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.
z, dim_, 1 );
638 Teuchos::BLAS<OT,ScalarType> blas;
639 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
640 Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
641 gmresState.
R->values(), gmresState.
R->stride(),
642 y_.values(), y_.stride() );
648 H_ = Teuchos::SerialDenseMatrix<OT,ScalarType>(Teuchos::Copy, *gmresState.
H, dim_, dim_);
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();
656 Teuchos::SerialDenseMatrix<OT,ScalarType> Htemp (Teuchos::Copy, H_, dim_, dim_);
659 ScalarType Hlast = (*gmresState.
H)(dim_,dim_-1);
660 Teuchos::SerialDenseMatrix<OT,ScalarType> HlastCol (Teuchos::View, H_, dim_, 1, 0, dim_-1);
663 Teuchos::SerialDenseMatrix< OT, ScalarType > F(dim_,1), E(dim_,1);
664 E.putScalar(SCT::zero());
665 E(dim_-1,0) = SCT::one();
667 Teuchos::SerialDenseSolver< OT, ScalarType > HSolver;
668 HSolver.setMatrix( Teuchos::rcpFromRef(Htemp));
669 HSolver.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
670 HSolver.setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
671 HSolver.factorWithEquilibration(
true );
675 info = HSolver.factor();
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);
689 Teuchos::LAPACK< OT, ScalarType > lapack;
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);
702 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
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;
713 const MagnitudeType tol = 10.0 * Teuchos::ScalarTraits<MagnitudeType>::eps();
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) );
740 const MagnitudeType one(1.0);
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())
755 extra[i] = ceil((log10(pof[i])-MagnitudeType(4.0))/MagnitudeType(14.0));
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);
770 Teuchos::SerialDenseMatrix<OT,MagnitudeType> thetaPert (Teuchos::Copy, theta_, 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>
809 void GmresPolyOp<ScalarType, MV, OP>::SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index)
const
814 int dimN = index.size();
815 std::vector<int> newIndex(dimN);
816 Teuchos::SerialDenseMatrix< OT, MagnitudeType > sorted (thetaN.numRows(), thetaN.numCols());
817 Teuchos::SerialDenseVector< OT, MagnitudeType > absVal (thetaN.numRows());
818 Teuchos::SerialDenseVector< OT, MagnitudeType > prod (thetaN.numRows());
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>
910 Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
911 Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
914 if (!LP_.is_null()) {
915 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
916 problem_->applyLeftPrec( *AX, *Xtmp );
921#ifdef BELOS_TEUCHOS_TIME_MONITOR
922 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
924 MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y);
926 for(
int i=1; i < dim_; i++)
928 Teuchos::RCP<MV> X, Y;
939 problem_->apply(*X, *Y);
941#ifdef BELOS_TEUCHOS_TIME_MONITOR
942 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
944 MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y);
949 if (!RP_.is_null()) {
950 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
951 problem_->applyRightPrec( *Ytmp, y );
955 template <
class ScalarType,
class MV,
class OP>
958 MVT::MvInit( y, SCT::zero() );
959 Teuchos::RCP<MV> prod = MVT::CloneCopy(x);
960 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
961 Teuchos::RCP<MV> Xtmp2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
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
976 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
978 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y);
980 problem_->apply(*prod, *Xtmp);
982#ifdef BELOS_TEUCHOS_TIME_MONITOR
983 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
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
995 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
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
1005 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
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
1016 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1018 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y);
1022 if (!RP_.is_null()) {
1023 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
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) {
1052 Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
1053 Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
1054 if (!LP_.is_null()) {
1055 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
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; }
1066 Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
1069 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1071 Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
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 );
1094 Teuchos::SerialDenseMatrix<OT,ScalarType> h(Teuchos::View,H_,i+1,1,0,i);
1096#ifdef BELOS_TEUCHOS_TIME_MONITOR
1097 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
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
1110 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1112 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1114 problem_->applyRightPrec( *wR_, *y_view );
1117#ifdef BELOS_TEUCHOS_TIME_MONITOR
1118 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1120 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
Belos concrete class for performing the block GMRES iteration.
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.
Declaration of basic traits for the multivector type.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Alternative run-time polymorphic interface for operators.
Class which manages the output and verbosity of the Belos solvers.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
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 MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
void MvAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const ScalarType beta, const MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
Teuchos::RCP< const MV > getConstMV() const
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
void MvRandom()
Fill all the vectors in *this with random numbers.
void SetBlock(const MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
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...
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
GmresPolyMv(const Teuchos::RCP< MV > &mv_in)
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 MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
Teuchos::RCP< MV > getMV()
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).
GmresPolyMv * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep 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).
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...
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
GmresPolyOpOrthoFailure(const std::string &what_arg)
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms_in)
Process the passed in parameters.
void ApplyRootsPoly(const MV &x, MV &y) const
void generateGmresPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in)
Given no ParameterList, constructor creates no polynomial and only applies the given operator.
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 ApplyGmresPoly(const MV &x, MV &y) const
void generateArnoldiPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
void ApplyArnoldiPoly(const MV &x, MV &y) const
virtual ~GmresPolyOp()
Destructor.
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...
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in, const Teuchos::RCP< Teuchos::ParameterList > ¶ms_in)
Basic contstructor.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
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 .
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
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 ,...
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 ...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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).
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 void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Interface for multivectors used by Belos' linear solvers.
Alternative run-time polymorphic interface for operators.
Enumeration of all valid Belos (Mat)OrthoManager classes.
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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
A Belos::StatusTest class for specifying a maximum number of iterations.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
NormType
The type of vector norm to compute.
ETrans
Whether to apply the (conjugate) transpose of an operator.
static const double polyTol
Relative residual tolerance for matrix polynomial construction.
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.