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" 81 #endif // BELOS_TEUCHOS_TIME_MONITOR 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_);
215 #endif // BELOS_TEUCHOS_TIME_MONITOR 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_;
289 #endif // BELOS_TEUCHOS_TIME_MONITOR 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";
306 #if defined(_WIN32) && defined(__clang__) 307 static constexpr std::ostream * outputStream_default_ =
308 __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
310 static constexpr std::ostream * outputStream_default_ = &std::cout;
312 static constexpr
bool damp_default_ =
false;
313 static constexpr
bool addRoots_default_ =
true;
316 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
317 Teuchos::RCP<Teuchos::ParameterList> params_;
318 Teuchos::RCP<const OP> LP_, RP_;
321 Teuchos::RCP<OutputManager<ScalarType> > printer_;
322 Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcp(outputStream_default_,
false);
325 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
329 int maxDegree_ = maxDegree_default_;
330 int verbosity_ = verbosity_default_;
331 bool randomRHS_ = randomRHS_default_;
332 std::string label_ = label_default_;
333 std::string polyType_ = polyType_default_;
334 std::string orthoType_ = orthoType_default_;
336 bool damp_ = damp_default_;
337 bool addRoots_ = addRoots_default_;
340 mutable Teuchos::RCP<MV> V_, wL_, wR_;
341 Teuchos::SerialDenseMatrix<OT,ScalarType> H_, y_;
342 Teuchos::SerialDenseVector<OT,ScalarType> r0_;
345 bool autoDeg =
false;
346 Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_;
349 Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_;
353 void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index)
const ;
356 void ComputeAddedRoots();
359 template <
class ScalarType,
class MV,
class OP>
363 if (params_in->isParameter(
"Polynomial Type")) {
364 polyType_ = params_in->get(
"Polynomial Type", polyType_default_);
368 if (params_in->isParameter(
"Polynomial Tolerance")) {
369 if (params_in->isType<MagnitudeType> (
"Polynomial Tolerance")) {
370 polyTol_ = params_in->get (
"Polynomial Tolerance",
379 if (params_in->isParameter(
"Maximum Degree")) {
380 maxDegree_ = params_in->get(
"Maximum Degree", maxDegree_default_);
384 if (params_in->isParameter(
"Random RHS")) {
385 randomRHS_ = params_in->get(
"Random RHS", randomRHS_default_);
389 if (params_in->isParameter(
"Verbosity")) {
390 if (Teuchos::isParameterType<int>(*params_in,
"Verbosity")) {
391 verbosity_ = params_in->get(
"Verbosity", verbosity_default_);
394 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,
"Verbosity");
398 if (params_in->isParameter(
"Orthogonalization")) {
399 orthoType_ = params_in->get(
"Orthogonalization",orthoType_default_);
403 if (params_in->isParameter(
"Timer Label")) {
404 label_ = params_in->get(
"Timer Label", label_default_);
408 if (params_in->isParameter(
"Output Stream")) {
409 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,
"Output Stream");
413 if (params_in->isParameter(
"Damped Poly")) {
414 damp_ = params_in->get(
"Damped Poly", damp_default_);
418 if (params_in->isParameter(
"Add Roots")) {
419 addRoots_ = params_in->get(
"Add Roots", addRoots_default_);
423 template <
class ScalarType,
class MV,
class OP>
426 Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+1 );
429 std::vector<int> index(1,0);
430 Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
432 MVT::MvRandom( *V0 );
434 MVT::Assign( *problem_->getRHS(), *V0 );
436 if ( !LP_.is_null() ) {
437 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
438 problem_->applyLeftPrec( *Vtemp, *V0);
441 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
442 problem_->apply( *Vtemp, *V0);
445 for(
int i=0; i< maxDegree_; i++)
448 Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
450 Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
451 problem_->apply( *Vi, *Vip1);
455 Teuchos::Range1D range( 1, maxDegree_);
456 Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
459 Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_, maxDegree_);
460 MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
463 Teuchos::LAPACK< OT, ScalarType > lapack;
468 Teuchos::SerialDenseMatrix< OT, ScalarType > lhs;
469 while( status && dim_ >= 1)
471 Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_, dim_);
472 lapack.POTRF(
'U', dim_, lhstemp.values(), lhstemp.stride(), &infoInt);
479 std::cout <<
"BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
480 std::cout <<
"Error code: " << infoInt << std::endl;
501 pCoeff_.shape( 1, 1);
502 pCoeff_(0,0) = SCT::one();
503 std::cout <<
"Poly Degree is zero. No preconditioner created." << std::endl;
507 pCoeff_.shape( dim_, 1);
509 Teuchos::Range1D rangeSub( 1, dim_);
510 Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
513 MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
514 lapack.POTRS(
'U', dim_, 1, lhs.values(), lhs.stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
517 std::cout <<
"BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
518 std::cout <<
"Error code: " << infoInt << std::endl;
523 template <
class ScalarType,
class MV,
class OP>
526 std::string polyLabel = label_ +
": GmresPolyOp creation";
529 std::vector<int> idx(1,0);
530 Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
531 Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
532 MVT::MvInit( *newX, SCT::zero() );
534 MVT::MvRandom( *newB );
537 MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
539 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
541 newProblem->setInitResVec( newB );
542 newProblem->setLeftPrec( problem_->getLeftPrec() );
543 newProblem->setRightPrec( problem_->getRightPrec() );
544 newProblem->setLabel(polyLabel);
545 newProblem->setProblem();
546 newProblem->setLSIndex( idx );
549 Teuchos::ParameterList polyList;
552 polyList.set(
"Num Blocks",maxDegree_);
553 polyList.set(
"Block Size",1);
554 polyList.set(
"Keep Hessenberg",
true);
560 if (ortho_.is_null()) {
561 params_->set(
"Orthogonalization", orthoType_);
563 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
565 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
569 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
573 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
578 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
582 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
586 Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
587 if ( !LP_.is_null() )
588 newProblem->applyLeftPrec( *newB, *V_0 );
591 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V_0);
592 newProblem->apply( *Vtemp, *V_0 );
599 int rank = ortho_->normalize( *V_0, Teuchos::rcpFromRef(r0_) );
601 "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
606 newstate.
z = Teuchos::rcpFromRef( r0_);
608 gmres_iter->initializeGmres(newstate);
612 gmres_iter->iterate();
616 gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
618 catch (std::exception& e) {
620 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGmresIter::iterate() at iteration " 621 << gmres_iter->getNumIters() << endl << e.what () << endl;
626 Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
636 if(polyType_ ==
"Arnoldi"){
639 y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.
z, dim_, 1 );
645 Teuchos::BLAS<OT,ScalarType> blas;
646 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
647 Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
648 gmresState.
R->values(), gmresState.
R->stride(),
649 y_.values(), y_.stride() );
655 H_ = Teuchos::SerialDenseMatrix<OT,ScalarType>(Teuchos::Copy, *gmresState.
H, dim_, dim_);
657 for(
int i=0; i <= dim_-3; i++) {
658 for(
int k=i+2; k <= dim_-1; k++) {
659 H_(k,i) = SCT::zero();
663 Teuchos::SerialDenseMatrix<OT,ScalarType> Htemp (Teuchos::Copy, H_, dim_, dim_);
666 ScalarType Hlast = (*gmresState.
H)(dim_,dim_-1);
667 Teuchos::SerialDenseMatrix<OT,ScalarType> HlastCol (Teuchos::View, H_, dim_, 1, 0, dim_-1);
670 Teuchos::SerialDenseMatrix< OT, ScalarType > F(dim_,1), E(dim_,1);
671 E.putScalar(SCT::zero());
672 E(dim_-1,0) = SCT::one();
674 Teuchos::SerialDenseSolver< OT, ScalarType > HSolver;
675 HSolver.setMatrix( Teuchos::rcpFromRef(Htemp));
676 HSolver.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
677 HSolver.setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
678 HSolver.factorWithEquilibration(
true );
682 info = HSolver.factor();
684 std::cout <<
"Hsolver factor: info = " << info << std::endl;
686 info = HSolver.solve();
688 std::cout <<
"Hsolver solve : info = " << info << std::endl;
692 F.scale(Hlast*Hlast);
696 Teuchos::LAPACK< OT, ScalarType > lapack;
697 theta_.shape(dim_,2);
704 std::vector<ScalarType> work(1);
705 std::vector<MagnitudeType> rwork(2*dim_);
708 lapack.GEEV(
'N',
'N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
709 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
710 work.resize( lwork );
712 lapack.GEEV(
'N',
'N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
715 std::cout <<
"GEEV solve : info = " << info << std::endl;
720 const MagnitudeType tol = 10.0 * Teuchos::ScalarTraits<MagnitudeType>::eps();
721 std::vector<int> index(dim_);
722 for(
int i=0; i<dim_; ++i){
725 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.");
727 SortModLeja(theta_,index);
736 template <
class ScalarType,
class MV,
class OP>
741 std::vector<std::complex<MagnitudeType>> cmplxHRitz (dim_);
742 for(
unsigned int i=0; i<cmplxHRitz.size(); ++i){
743 cmplxHRitz[i] = std::complex<MagnitudeType>( theta_(i,0), theta_(i,1) );
747 const MagnitudeType one(1.0);
748 std::vector<MagnitudeType> pof (dim_,one);
749 for(
int j=0; j<dim_; ++j) {
750 for(
int i=0; i<dim_; ++i) {
752 pof[j] = std::abs(pof[j]*(one-(cmplxHRitz[j]/cmplxHRitz[i])));
758 std::vector<int> extra (dim_);
760 for(
int i=0; i<dim_; ++i){
761 if (pof[i] > MCT::zero())
762 extra[i] = ceil((log10(pof[i])-MagnitudeType(4.0))/MagnitudeType(14.0));
766 totalExtra += extra[i];
770 printer_->stream(
Warnings) <<
"Warning: Need to add " << totalExtra <<
" extra roots." << std::endl;}
773 if(addRoots_ && totalExtra>0)
775 theta_.reshape(dim_+totalExtra,2);
777 Teuchos::SerialDenseMatrix<OT,MagnitudeType> thetaPert (Teuchos::Copy, theta_, dim_+totalExtra, 2);
781 for(
int i=0; i<dim_; ++i){
782 for(
int j=0; j< extra[i]; ++j){
783 theta_(count,0) = theta_(i,0);
784 theta_(count,1) = theta_(i,1);
785 thetaPert(count,0) = theta_(i,0)+(j+MCT::one())*MagnitudeType(5e-8);
786 thetaPert(count,1) = theta_(i,1);
794 printer_->stream(
Warnings) <<
"New poly degree is: " << dim_ << std::endl;}
797 std::vector<int> index2(dim_);
798 for(
int i=0; i<dim_; ++i){
801 SortModLeja(thetaPert,index2);
803 for(
int i=0; i<dim_; ++i)
805 thetaPert(i,0) = theta_(index2[i],0);
806 thetaPert(i,1) = theta_(index2[i],1);
815 template <
class ScalarType,
class MV,
class OP>
816 void GmresPolyOp<ScalarType, MV, OP>::SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index)
const 821 int dimN = index.size();
822 std::vector<int> newIndex(dimN);
823 Teuchos::SerialDenseMatrix< OT, MagnitudeType > sorted (thetaN.numRows(), thetaN.numCols());
824 Teuchos::SerialDenseVector< OT, MagnitudeType > absVal (thetaN.numRows());
825 Teuchos::SerialDenseVector< OT, MagnitudeType > prod (thetaN.numRows());
828 for(
int i = 0; i < dimN; i++){
829 absVal(i) = hypot(thetaN(i,0), thetaN(i,1));
831 MagnitudeType * maxPointer = std::max_element(absVal.values(), (absVal.values()+dimN));
832 int maxIndex = int (maxPointer- absVal.values());
835 sorted(0,0) = thetaN(maxIndex,0);
836 sorted(0,1) = thetaN(maxIndex,1);
837 newIndex[0] = index[maxIndex];
841 if(sorted(0,1)!= SCT::zero() && !SCT::isComplex)
843 sorted(1,0) = thetaN(maxIndex,0);
844 sorted(1,1) = -thetaN(maxIndex,1);
845 newIndex[1] = index[maxIndex+1];
858 for(
int i = 0; i < dimN; i++)
860 prod(i) = MCT::one();
861 for(
int k = 0; k < j; k++)
863 a = thetaN(i,0) - sorted(k,0);
864 b = thetaN(i,1) - sorted(k,1);
865 if (a*a + b*b > MCT::zero())
866 prod(i) = prod(i) + log10(hypot(a,b));
868 prod(i) = -std::numeric_limits<MagnitudeType>::infinity();
875 MagnitudeType * maxPointer = std::max_element(prod.values(), (prod.values()+dimN));
876 int maxIndex = int (maxPointer- prod.values());
877 sorted(j,0) = thetaN(maxIndex,0);
878 sorted(j,1) = thetaN(maxIndex,1);
879 newIndex[j] = index[maxIndex];
882 if(sorted(j,1)!= SCT::zero() && !SCT::isComplex)
885 sorted(j,0) = thetaN(maxIndex,0);
886 sorted(j,1) = -thetaN(maxIndex,1);
887 newIndex[j] = index[maxIndex+1];
897 template <
class ScalarType,
class MV,
class OP>
901 if (polyType_ ==
"Arnoldi")
902 ApplyArnoldiPoly(x, y);
903 else if (polyType_ ==
"Gmres")
904 ApplyGmresPoly(x, y);
905 else if (polyType_ ==
"Roots")
906 ApplyRootsPoly(x, y);
910 problem_->applyOp( x, y );
914 template <
class ScalarType,
class MV,
class OP>
917 Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
918 Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
921 if (!LP_.is_null()) {
922 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
923 problem_->applyLeftPrec( *AX, *Xtmp );
928 #ifdef BELOS_TEUCHOS_TIME_MONITOR 929 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
931 MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y);
933 for(
int i=1; i < dim_; i++)
935 Teuchos::RCP<MV> X, Y;
946 problem_->apply(*X, *Y);
948 #ifdef BELOS_TEUCHOS_TIME_MONITOR 949 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
951 MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y);
956 if (!RP_.is_null()) {
957 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
958 problem_->applyRightPrec( *Ytmp, y );
962 template <
class ScalarType,
class MV,
class OP>
965 MVT::MvInit( y, SCT::zero() );
966 Teuchos::RCP<MV> prod = MVT::CloneCopy(x);
967 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
968 Teuchos::RCP<MV> Xtmp2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
971 if (!LP_.is_null()) {
972 problem_->applyLeftPrec( *prod, *Xtmp );
979 if(theta_(i,1)== SCT::zero() || SCT::isComplex)
982 #ifdef BELOS_TEUCHOS_TIME_MONITOR 983 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
985 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y);
987 problem_->apply(*prod, *Xtmp);
989 #ifdef BELOS_TEUCHOS_TIME_MONITOR 990 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
992 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/theta_(i,0), *Xtmp, *prod);
998 MagnitudeType mod = theta_(i,0)*theta_(i,0) + theta_(i,1)*theta_(i,1);
999 problem_->apply(*prod, *Xtmp);
1001 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1002 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1004 MVT::MvAddMv(2*theta_(i,0), *prod, -SCT::one(), *Xtmp, *Xtmp);
1005 MVT::MvAddMv(SCT::one(), y, SCT::one()/mod, *Xtmp, y);
1009 problem_->apply(*Xtmp, *Xtmp2);
1011 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1012 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1014 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/mod, *Xtmp2, *prod);
1020 if(theta_(dim_-1,1)== SCT::zero() || SCT::isComplex)
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1023 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1025 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y);
1029 if (!RP_.is_null()) {
1030 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
1031 problem_->applyRightPrec( *Ytmp, y );
1035 template <
class ScalarType,
class MV,
class OP>
1040 V_ = MVT::Clone( x, dim_ );
1041 if (!LP_.is_null()) {
1042 wL_ = MVT::Clone( y, 1 );
1044 if (!RP_.is_null()) {
1045 wR_ = MVT::Clone( y, 1 );
1051 int n = MVT::GetNumberVecs( x );
1052 std::vector<int> idxi(1), idxi2, idxj(1);
1055 for (
int j=0; j<n; ++j) {
1059 Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
1060 Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
1061 if (!LP_.is_null()) {
1062 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1063 problem_->applyLeftPrec( *x_view, *v_curr );
1065 MVT::SetBlock( *x_view, idxi, *V_ );
1068 for (
int i=0; i<dim_-1; ++i) {
1072 for (
int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
1073 Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
1076 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1078 Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
1084 if (!RP_.is_null()) {
1085 problem_->applyRightPrec( *v_curr, *wR_ );
1090 if (LP_.is_null()) {
1094 problem_->applyOp( *wR_, *wL_ );
1096 if (!LP_.is_null()) {
1097 problem_->applyLeftPrec( *wL_, *v_next );
1101 Teuchos::SerialDenseMatrix<OT,ScalarType> h(Teuchos::View,H_,i+1,1,0,i);
1103 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1104 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1106 MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
1110 MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
1114 if (!RP_.is_null()) {
1116 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1117 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1119 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1121 problem_->applyRightPrec( *wR_, *y_view );
1124 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1125 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1127 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.
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
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...
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.
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...
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.
void ApplyGmresPoly(const MV &x, MV &y) const
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
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()
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...
Declaration of basic traits for the multivector type.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
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...
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.
GmresPolyMv * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
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...
void MvRandom()
Fill all the vectors in *this with random numbers.
Structure to contain pointers to GmresIteration state variables.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< const MV > getConstMV() const
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
ETrans
Whether to apply the (conjugate) transpose of an operator.
void ApplyRootsPoly(const MV &x, MV &y) const
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
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.
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.
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.
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 curDim
The current dimension of the reduction.
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
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< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
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...
A linear system to solve, and its associated information.
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.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
void ApplyArnoldiPoly(const MV &x, MV &y) const
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
Belos concrete class for performing the block GMRES iteration.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
NormType
The type of vector norm to compute.
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
Alternative run-time polymorphic interface for operators.
virtual ~GmresPolyOp()
Destructor.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
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.
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.
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms_in)
Process the passed in parameters.
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
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.