4 #ifndef ANASAZI_IRTR_HPP 5 #define ANASAZI_IRTR_HPP 13 #include "Teuchos_ScalarTraits.hpp" 15 #include "Teuchos_LAPACK.hpp" 16 #include "Teuchos_BLAS.hpp" 17 #include "Teuchos_SerialDenseMatrix.hpp" 18 #include "Teuchos_ParameterList.hpp" 19 #include "Teuchos_TimeMonitor.hpp" 45 template <
class ScalarType,
class MV,
class OP>
64 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
68 Teuchos::ParameterList ¶ms
99 typedef Teuchos::ScalarTraits<ScalarType> SCT;
100 typedef typename SCT::magnitudeType MagnitudeType;
101 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
111 std::vector<std::string> stopReasons_;
115 const MagnitudeType ZERO;
116 const MagnitudeType ONE;
121 void solveTRSubproblem();
124 MagnitudeType rho_prime_;
127 MagnitudeType normgradf0_;
130 trRetType innerStop_;
133 int innerIters_, totalInnerIters_;
144 template <
class ScalarType,
class MV,
class OP>
147 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
151 Teuchos::ParameterList ¶ms
153 RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,
"IRTR",false),
159 stopReasons_.push_back(
"n/a");
160 stopReasons_.push_back(
"maximum iterations");
161 stopReasons_.push_back(
"negative curvature");
162 stopReasons_.push_back(
"exceeded TR");
163 stopReasons_.push_back(
"kappa convergence");
164 stopReasons_.push_back(
"theta convergence");
166 rho_prime_ = params.get(
"Rho Prime",0.5);
167 TEUCHOS_TEST_FOR_EXCEPTION(rho_prime_ <= 0 || rho_prime_ >= 1,std::invalid_argument,
168 "Anasazi::IRTR::constructor: rho_prime must be in (0,1).");
170 useSA_ = params.get<
bool>(
"Use SA",
false);
183 template <
class ScalarType,
class MV,
class OP>
194 using Teuchos::tuple;
196 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 197 using Teuchos::TimeMonitor;
200 typedef Teuchos::RCP<const MV> PCMV;
201 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
203 innerStop_ = MAXIMUM_ITERATIONS;
206 const int p = this->blockSize_;
207 const int d = n*p - (p*p+p)/2;
221 const double D2 = ONE/rho_prime_ - ONE;
223 std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p);
224 std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p);
225 MagnitudeType r0_norm;
246 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 247 TimeMonitor lcltimer( *this->timerOrtho_ );
249 this->orthman_->projectGen(
251 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
253 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
254 if (this->leftMost_) {
267 MagnitudeType kconv = r0_norm * this->conv_kappa_;
270 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
271 if (this->om_->isVerbosity(
Debug)) {
272 this->om_->stream(
Debug)
273 <<
" >> |r0| : " << r0_norm << endl
274 <<
" >> kappa conv : " << kconv << endl
275 <<
" >> theta conv : " << tconv << endl;
283 if (this->hasPrec_ && this->olsenPrec_)
285 std::vector<int> ind(this->blockSize_);
286 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
289 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 290 TimeMonitor prectimer( *this->timerPrec_ );
293 this->counterPrec_ += this->blockSize_;
304 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 305 TimeMonitor prectimer( *this->timerPrec_ );
308 this->counterPrec_ += this->blockSize_;
310 if (this->olsenPrec_) {
311 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 312 TimeMonitor orthtimer( *this->timerOrtho_ );
314 this->orthman_->projectGen(
316 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
318 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
321 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 322 TimeMonitor orthtimer( *this->timerOrtho_ );
324 this->orthman_->projectGen(
326 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
328 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
336 if (this->om_->isVerbosity(
Debug )) {
340 if (this->hasPrec_) chk.checkZ =
true;
341 this->om_->print(
Debug, this->accuracyCheck(chk,
"after computing gradient") );
347 if (this->hasPrec_) chk.checkZ =
true;
348 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after computing gradient") );
352 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
354 if (this->om_->isVerbosity(
Debug)) {
355 RCP<MV> Heta =
MVT::Clone(*this->eta_,this->blockSize_);
358 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
360 if (this->leftMost_) {
368 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 369 TimeMonitor lcltimer( *this->timerOrtho_ );
371 this->orthman_->projectGen(
373 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
375 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
378 std::vector<MagnitudeType> eg(this->blockSize_),
379 eHe(this->blockSize_);
382 if (this->leftMost_) {
383 for (
int j=0; j<this->blockSize_; ++j) {
384 eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
388 for (
int j=0; j<this->blockSize_; ++j) {
389 eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
392 this->om_->stream(
Debug)
393 <<
" Debugging checks: IRTR inner iteration " << innerIters_ << endl
394 <<
" >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
395 for (
int j=0; j<this->blockSize_; ++j) {
396 this->om_->stream(
Debug)
397 <<
" >> m_X(eta_" << j <<
") : " << eg[j] << endl;
404 for (innerIters_=1; innerIters_<=d; ++innerIters_) {
412 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 413 TimeMonitor lcltimer( *this->timerAOp_ );
415 OPT::Apply(*this->AOp_,*this->delta_,*this->Adelta_);
416 this->counterAOp_ += this->blockSize_;
419 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 420 TimeMonitor lcltimer( *this->timerBOp_ );
422 OPT::Apply(*this->BOp_,*this->delta_,*this->Bdelta_);
423 this->counterBOp_ += this->blockSize_;
426 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Bdelta_);
433 MVT::MvAddMv(ONE,*this->Bdelta_,ZERO,*this->Bdelta_,*this->Hdelta_);
435 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
438 if (this->leftMost_) {
439 MVT::MvAddMv( 2.0,*this->Adelta_,-2.0,*this->Hdelta_,*this->Hdelta_);
442 MVT::MvAddMv(-2.0,*this->Adelta_, 2.0,*this->Hdelta_,*this->Hdelta_);
446 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 447 TimeMonitor lcltimer( *this->timerOrtho_ );
449 this->orthman_->projectGen(
451 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
453 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
459 for (
unsigned int j=0; j<alpha.size(); ++j)
461 alpha[j] = z_r[j]/d_Hd[j];
465 for (
unsigned int j=0; j<alpha.size(); ++j)
467 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
470 if (this->om_->isVerbosity(
Debug)) {
471 for (
unsigned int j=0; j<alpha.size(); j++) {
472 this->om_->stream(
Debug)
473 <<
" >> z_r[" << j <<
"] : " << z_r[j]
474 <<
" d_Hd[" << j <<
"] : " << d_Hd[j] << endl
475 <<
" >> eBe[" << j <<
"] : " << eBe[j]
476 <<
" neweBe[" << j <<
"] : " << new_eBe[j] << endl
477 <<
" >> eBd[" << j <<
"] : " << eBd[j]
478 <<
" dBd[" << j <<
"] : " << dBd[j] << endl;
483 std::vector<int> trncstep;
487 bool atleastonenegcur =
false;
488 for (
unsigned int j=0; j<d_Hd.size(); ++j) {
490 trncstep.push_back(j);
491 atleastonenegcur =
true;
493 else if (new_eBe[j] > D2) {
494 trncstep.push_back(j);
498 if (!trncstep.empty())
501 if (this->om_->isVerbosity(
Debug)) {
502 for (
unsigned int j=0; j<trncstep.size(); ++j) {
503 this->om_->stream(
Debug)
504 <<
" >> alpha[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
507 for (
unsigned int j=0; j<trncstep.size(); ++j) {
508 int jj = trncstep[j];
509 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
511 if (this->om_->isVerbosity(
Debug)) {
512 for (
unsigned int j=0; j<trncstep.size(); ++j) {
513 this->om_->stream(
Debug)
514 <<
" >> tau[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
517 if (atleastonenegcur) {
518 innerStop_ = NEGATIVE_CURVATURE;
521 innerStop_ = EXCEEDED_TR;
530 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
538 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
539 MVT::MvAddMv(ONE,*this->Adelta_,ONE,*this->Aeta_,*this->Aeta_);
541 MVT::MvAddMv(ONE,*this->Bdelta_,ONE,*this->Beta_,*this->Beta_);
549 if (this->om_->isVerbosity(
Debug)) {
550 RCP<MV> Heta =
MVT::Clone(*this->eta_,this->blockSize_);
554 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
557 if (this->leftMost_) {
565 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 566 TimeMonitor lcltimer( *this->timerOrtho_ );
568 this->orthman_->projectGen(
570 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
572 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
575 std::vector<MagnitudeType> eg(this->blockSize_),
576 eHe(this->blockSize_);
579 if (this->leftMost_) {
580 for (
int j=0; j<this->blockSize_; ++j) {
581 eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
585 for (
int j=0; j<this->blockSize_; ++j) {
586 eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
589 this->om_->stream(
Debug)
590 <<
" Debugging checks: IRTR inner iteration " << innerIters_ << endl
591 <<
" >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
592 for (
int j=0; j<this->blockSize_; ++j) {
593 this->om_->stream(
Debug)
594 <<
" >> m_X(eta_" << j <<
") : " << eg[j] << endl;
600 if (!trncstep.empty()) {
608 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
611 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 612 TimeMonitor lcltimer( *this->timerOrtho_ );
614 this->orthman_->projectGen(
616 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
618 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
631 if (this->om_->isVerbosity(
Debug)) {
632 this->om_->stream(
Debug)
633 <<
" >> |r" << innerIters_ <<
"| : " << r_norm << endl;
635 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
636 if (tconv <= kconv) {
637 innerStop_ = THETA_CONVERGENCE;
640 innerStop_ = KAPPA_CONVERGENCE;
652 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 653 TimeMonitor prectimer( *this->timerPrec_ );
656 this->counterPrec_ += this->blockSize_;
658 if (this->olsenPrec_) {
659 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 660 TimeMonitor orthtimer( *this->timerOrtho_ );
662 this->orthman_->projectGen(
664 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
666 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
669 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 670 TimeMonitor orthtimer( *this->timerOrtho_ );
672 this->orthman_->projectGen(
674 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
676 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
697 for (
unsigned int j=0; j<beta.size(); ++j) {
698 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
701 std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
704 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
710 if (innerIters_ > d) innerIters_ = d;
712 this->om_->stream(
Debug)
713 <<
" >> stop reason is " << stopReasons_[innerStop_] << endl
721 template <
class ScalarType,
class MV,
class OP>
726 using Teuchos::tuple;
727 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 728 using Teuchos::TimeMonitor;
737 if (this->initialized_ ==
false) {
741 Teuchos::SerialDenseMatrix<int,ScalarType> AA, BB, S;
742 if (useSA_ ==
true) {
743 AA.reshape(2*this->blockSize_,2*this->blockSize_);
744 BB.reshape(2*this->blockSize_,2*this->blockSize_);
745 S.reshape(2*this->blockSize_,2*this->blockSize_);
748 AA.reshape(this->blockSize_,this->blockSize_);
749 BB.reshape(this->blockSize_,this->blockSize_);
750 S.reshape(this->blockSize_,this->blockSize_);
755 innerStop_ = UNINITIALIZED;
758 while (this->tester_->checkStatus(
this) !=
Passed) {
761 if (this->om_->isVerbosity(
Debug)) {
773 totalInnerIters_ += innerIters_;
776 if (this->om_->isVerbosity(
Debug ) ) {
781 chk.checkAEta =
true;
782 chk.checkBEta =
true;
783 this->om_->print(
Debug, this->accuracyCheck(chk,
"in iterate() after solveTRSubproblem()") );
784 this->om_->stream(
Debug)
789 if (useSA_ ==
false) {
793 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 794 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
796 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*this->delta_);
797 MVT::MvAddMv(ONE,*this->AX_,ONE,*this->Aeta_,*this->Adelta_);
799 MVT::MvAddMv(ONE,*this->BX_,ONE,*this->Beta_,*this->Bdelta_);
803 if (this->om_->isVerbosity(
Debug ) ) {
807 Teuchos::SerialDenseMatrix<int,ScalarType> XE(this->blockSize_,this->blockSize_),
808 E(this->blockSize_,this->blockSize_);
811 this->om_->stream(
Debug)
812 <<
" >> Error in AX+AEta == A(X+Eta) : " <<
Utils::errorEquality(*this->delta_,*this->Adelta_,this->AOp_) << endl
813 <<
" >> Error in BX+BEta == B(X+Eta) : " <<
Utils::errorEquality(*this->delta_,*this->Bdelta_,this->BOp_) << endl
814 <<
" >> norm( (X+Eta)^T B (X+Eta) ) : " << XE.normFrobenius() << endl
815 <<
" >> norm( Eta^T B Eta ) : " << E.normFrobenius() << endl
824 MagnitudeType oldfx = this->fx_;
825 std::vector<MagnitudeType> oldtheta(this->theta_), newtheta(AA.numRows());
828 if (useSA_ ==
true) {
829 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 830 TimeMonitor lcltimer( *this->timerLocalProj_ );
832 Teuchos::SerialDenseMatrix<int,ScalarType> AA11(Teuchos::View,AA,this->blockSize_,this->blockSize_,0,0),
833 AA12(Teuchos::View,AA,this->blockSize_,this->blockSize_,0,this->blockSize_),
834 AA22(Teuchos::View,AA,this->blockSize_,this->blockSize_,this->blockSize_,this->blockSize_);
835 Teuchos::SerialDenseMatrix<int,ScalarType> BB11(Teuchos::View,BB,this->blockSize_,this->blockSize_,0,0),
836 BB12(Teuchos::View,BB,this->blockSize_,this->blockSize_,0,this->blockSize_),
837 BB22(Teuchos::View,BB,this->blockSize_,this->blockSize_,this->blockSize_,this->blockSize_);
846 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 847 TimeMonitor lcltimer( *this->timerLocalProj_ );
852 this->om_->stream(
Debug) <<
"AA: " << std::endl << AA << std::endl;;
853 this->om_->stream(
Debug) <<
"BB: " << std::endl << BB << std::endl;;
855 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 856 TimeMonitor lcltimer( *this->timerDS_ );
860 this->om_->stream(
Debug) <<
"S: " << std::endl << S << std::endl;;
861 TEUCHOS_TEST_FOR_EXCEPTION(ret != 0,std::logic_error,
"Anasazi::IRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret);
862 TEUCHOS_TEST_FOR_EXCEPTION(rank != AA.numRows(),
RTRRitzFailure,
"Anasazi::IRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank);
868 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 869 TimeMonitor lcltimer( *this->timerSort_ );
871 std::vector<int> order(newtheta.size());
873 this->sm_->sort(newtheta, Teuchos::rcpFromRef(order), -1);
879 std::copy(newtheta.begin(), newtheta.begin()+this->blockSize_, this->theta_.begin());
882 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
886 if (this->om_->isVerbosity(
Debug ) ) {
897 MagnitudeType rhonum, rhoden, mxeta;
898 std::vector<MagnitudeType> eBe(this->blockSize_);
902 rhonum = oldfx - this->fx_;
907 for (
int i=0; i<this->blockSize_; ++i) {
908 rhoden += eBe[i]*oldtheta[i];
910 mxeta = oldfx - rhoden;
911 this->rho_ = rhonum / rhoden;
912 this->om_->stream(
Debug)
913 <<
" >> old f(x) is : " << oldfx << endl
914 <<
" >> new f(x) is : " << this->fx_ << endl
915 <<
" >> m_x(eta) is : " << mxeta << endl
916 <<
" >> rhonum is : " << rhonum << endl
917 <<
" >> rhoden is : " << rhoden << endl
918 <<
" >> rho is : " << this->rho_ << endl;
920 for (
int j=0; j<this->blockSize_; ++j) {
921 this->om_->stream(
Debug)
922 <<
" >> rho[" << j <<
"] : " << 1.0/(1.0+eBe[j]) << endl;
933 this->X_ = Teuchos::null;
934 this->BX_ = Teuchos::null;
935 std::vector<int> ind(this->blockSize_);
936 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
937 Teuchos::RCP<MV> X, BX;
942 if (useSA_ ==
false) {
946 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 947 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
960 Teuchos::SerialDenseMatrix<int,ScalarType> Sx(Teuchos::View,S,this->blockSize_,this->blockSize_,0,0),
961 Se(Teuchos::View,S,this->blockSize_,this->blockSize_,this->blockSize_,0);
962 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 963 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
972 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->AX_);
983 this->X_ =
MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
984 this->BX_ =
MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
990 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 991 TimeMonitor lcltimer( *this->timerCompRes_ );
993 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
994 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
996 MVT::MvAddMv( ONE, *this->AX_, -ONE, *this->R_, *this->R_ );
1000 this->Rnorms_current_ =
false;
1001 this->R2norms_current_ =
false;
1006 if (this->om_->isVerbosity(
Debug ) ) {
1013 this->om_->print(
Debug, this->accuracyCheck(chk,
"after local update") );
1019 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after local update") );
1029 template <
class ScalarType,
class MV,
class OP>
1034 os.setf(std::ios::scientific, std::ios::floatfield);
1037 os <<
"================================================================================" << endl;
1039 os <<
" IRTR Solver Status" << endl;
1041 os <<
"The solver is "<<(this->initialized_ ?
"initialized." :
"not initialized.") << endl;
1042 os <<
"The number of iterations performed is " << this->iter_ << endl;
1043 os <<
"The current block size is " << this->blockSize_ << endl;
1044 os <<
"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
1045 os <<
"The number of operations A*x is " << this->counterAOp_ << endl;
1046 os <<
"The number of operations B*x is " << this->counterBOp_ << endl;
1047 os <<
"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
1048 os <<
"The number of operations Prec*x is " << this->counterPrec_ << endl;
1049 os <<
"Parameter rho_prime is " << rho_prime_ << endl;
1050 os <<
"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
1051 os <<
"Number of inner iterations was " << innerIters_ << endl;
1052 os <<
"Total number of inner iterations is " << totalInnerIters_ << endl;
1053 os <<
"f(x) is " << this->fx_ << endl;
1055 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1057 if (this->initialized_) {
1059 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1060 os << std::setw(20) <<
"Eigenvalue" 1061 << std::setw(20) <<
"Residual(B)" 1062 << std::setw(20) <<
"Residual(2)" 1064 os <<
"--------------------------------------------------------------------------------"<<endl;
1065 for (
int i=0; i<this->blockSize_; i++) {
1066 os << std::setw(20) << this->theta_[i];
1067 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
1068 else os << std::setw(20) <<
"not current";
1069 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
1070 else os << std::setw(20) <<
"not current";
1074 os <<
"================================================================================" << endl;
1081 #endif // ANASAZI_IRTR_HPP static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace 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 .
This class defines the interface required by an eigensolver and status test class to compute solution...
Base class for Implicit Riemannian Trust-Region solvers.
Declaration of basic traits for the multivector type.
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers. The class provides the interfaces shared by the IRTR solvers (e.g., getState() and initialize()) as well as the shared implementations (e.g., inner products).
Virtual base class which defines basic traits for the operator type.
static Teuchos::ScalarTraits< ScalarType >::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP< const OP > M=Teuchos::null)
Return the maximum coefficient of the matrix scaled by the maximum coefficient of MX...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated, static class providing utilities for the solvers.
Output managers remove the need for the eigensolver to know any information about the required output...
static int directSolver(int size, const Teuchos::SerialDenseMatrix< int, ScalarType > &KK, Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > MM, Teuchos::SerialDenseMatrix< int, ScalarType > &EV, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &theta, int &nev, int esType=0)
Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK...
IRTR(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
IRTR constructor with eigenproblem, solver utilities, and parameter list of solver options...
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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< 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 iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
Types and exceptions used within Anasazi solvers and interfaces.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
virtual ~IRTR()
IRTR destructor
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
void currentStatus(std::ostream &os)
Impemements Eigensolver. This method requests that the solver print out its current status to screen...