44 #ifndef ANASAZI_SIRTR_HPP
45 #define ANASAZI_SIRTR_HPP
85 template <
class ScalarType,
class MV,
class OP>
151 std::vector<std::string> stopReasons_;
155 const MagnitudeType ZERO;
156 const MagnitudeType ONE;
161 void solveTRSubproblem();
164 MagnitudeType rho_prime_;
167 MagnitudeType normgradf0_;
170 trRetType innerStop_;
173 int innerIters_, totalInnerIters_;
181 template <
class ScalarType,
class MV,
class OP>
190 RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,
"SIRTR",true),
196 stopReasons_.push_back(
"n/a");
197 stopReasons_.push_back(
"maximum iterations");
198 stopReasons_.push_back(
"negative curvature");
199 stopReasons_.push_back(
"exceeded TR");
200 stopReasons_.push_back(
"kappa convergence");
201 stopReasons_.push_back(
"theta convergence");
203 rho_prime_ = params.
get(
"Rho Prime",0.5);
205 "Anasazi::SIRTR::constructor: rho_prime must be in (0,1).");
218 template <
class ScalarType,
class MV,
class OP>
229 using Teuchos::tuple;
231 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
238 innerStop_ = MAXIMUM_ITERATIONS;
240 const int n = MVT::GetGlobalLength(*this->eta_);
241 const int p = this->blockSize_;
242 const int d = n*p - (p*p+p)/2;
256 const double D2 = ONE/rho_prime_ - ONE;
258 std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p);
259 std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p);
260 MagnitudeType r0_norm;
262 MVT::MvInit(*this->eta_ ,0.0);
277 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
278 TimeMonitor lcltimer( *this->timerOrtho_ );
280 this->orthman_->projectGen(
282 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
284 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
285 if (this->leftMost_) {
286 MVT::MvScale(*this->R_,2.0);
289 MVT::MvScale(*this->R_,-2.0);
293 r0_norm = MAT::squareroot( RTRBase<ScalarType,MV,OP>::ginner(*this->R_) );
298 MagnitudeType kconv = r0_norm * this->conv_kappa_;
301 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
302 if (this->om_->isVerbosity(
Debug)) {
303 this->om_->stream(
Debug)
304 <<
" >> |r0| : " << r0_norm << endl
305 <<
" >> kappa conv : " << kconv << endl
306 <<
" >> theta conv : " << tconv << endl;
314 if (this->hasPrec_ && this->olsenPrec_)
316 std::vector<int> ind(this->blockSize_);
317 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
320 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321 TimeMonitor prectimer( *this->timerPrec_ );
323 OPT::Apply(*this->Prec_,*this->BX_,*PBX);
324 this->counterPrec_ += this->blockSize_;
335 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
336 TimeMonitor prectimer( *this->timerPrec_ );
338 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
339 this->counterPrec_ += this->blockSize_;
341 if (this->olsenPrec_) {
342 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
343 TimeMonitor orthtimer( *this->timerOrtho_ );
345 this->orthman_->projectGen(
347 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
349 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
352 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
353 TimeMonitor orthtimer( *this->timerOrtho_ );
355 this->orthman_->projectGen(
357 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
359 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
361 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
365 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
366 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
369 if (this->om_->isVerbosity(
Debug )) {
371 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
373 if (this->hasPrec_) chk.checkZ =
true;
374 this->om_->print(
Debug, this->accuracyCheck(chk,
"after computing gradient") );
378 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
380 if (this->hasPrec_) chk.checkZ =
true;
381 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after computing gradient") );
385 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
387 if (this->om_->isVerbosity(
Debug)) {
392 std::vector<MagnitudeType> eAx(this->blockSize_),
393 d_eAe(this->blockSize_),
394 d_eBe(this->blockSize_),
395 d_mxe(this->blockSize_);
398 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
399 TimeMonitor lcltimer( *this->timerAOp_ );
401 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
402 this->counterAOp_ += this->blockSize_;
404 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx);
407 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
408 TimeMonitor lcltimer( *this->timerAOp_ );
410 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
411 this->counterAOp_ += this->blockSize_;
413 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe);
416 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
417 TimeMonitor lcltimer( *this->timerBOp_ );
419 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
420 this->counterBOp_ += this->blockSize_;
423 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
425 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe);
428 if (this->leftMost_) {
429 for (
int j=0; j<this->blockSize_; ++j) {
430 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
434 for (
int j=0; j<this->blockSize_; ++j) {
435 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
438 this->om_->stream(
Debug)
439 <<
" Debugging checks: SIRTR inner iteration " << innerIters_ << endl
440 <<
" >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
441 for (
int j=0; j<this->blockSize_; ++j) {
442 this->om_->stream(
Debug)
443 <<
" >> m_X(eta_" << j <<
") : " << d_mxe[j] << endl;
450 for (innerIters_=1; innerIters_<=d; ++innerIters_) {
458 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
459 TimeMonitor lcltimer( *this->timerAOp_ );
461 OPT::Apply(*this->AOp_,*this->delta_,*this->Z_);
462 this->counterAOp_ += this->blockSize_;
465 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
466 TimeMonitor lcltimer( *this->timerBOp_ );
468 OPT::Apply(*this->BOp_,*this->delta_,*this->Hdelta_);
469 this->counterBOp_ += this->blockSize_;
472 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Hdelta_);
476 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_ ,*this->Hdelta_,eBd);
477 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,dBd);
480 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
481 MVT::MvScale(*this->Hdelta_,theta_comp);
483 if (this->leftMost_) {
484 MVT::MvAddMv( 2.0,*this->Z_,-2.0,*this->Hdelta_,*this->Hdelta_);
487 MVT::MvAddMv(-2.0,*this->Z_, 2.0,*this->Hdelta_,*this->Hdelta_);
491 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
492 TimeMonitor lcltimer( *this->timerOrtho_ );
494 this->orthman_->projectGen(
496 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
498 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
500 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
504 for (
unsigned int j=0; j<alpha.size(); ++j)
506 alpha[j] = z_r[j]/d_Hd[j];
510 for (
unsigned int j=0; j<alpha.size(); ++j)
512 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
515 if (this->om_->isVerbosity(
Debug)) {
516 for (
unsigned int j=0; j<alpha.size(); j++) {
517 this->om_->stream(
Debug)
518 <<
" >> z_r[" << j <<
"] : " << z_r[j]
519 <<
" d_Hd[" << j <<
"] : " << d_Hd[j] << endl
520 <<
" >> eBe[" << j <<
"] : " << eBe[j]
521 <<
" neweBe[" << j <<
"] : " << new_eBe[j] << endl
522 <<
" >> eBd[" << j <<
"] : " << eBd[j]
523 <<
" dBd[" << j <<
"] : " << dBd[j] << endl;
528 std::vector<int> trncstep;
532 bool atleastonenegcur =
false;
533 for (
unsigned int j=0; j<d_Hd.size(); ++j) {
535 trncstep.push_back(j);
536 atleastonenegcur =
true;
538 else if (new_eBe[j] > D2) {
539 trncstep.push_back(j);
543 if (!trncstep.empty())
546 if (this->om_->isVerbosity(
Debug)) {
547 for (
unsigned int j=0; j<trncstep.size(); ++j) {
548 this->om_->stream(
Debug)
549 <<
" >> alpha[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
552 for (
unsigned int j=0; j<trncstep.size(); ++j) {
553 int jj = trncstep[j];
554 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
556 if (this->om_->isVerbosity(
Debug)) {
557 for (
unsigned int j=0; j<trncstep.size(); ++j) {
558 this->om_->stream(
Debug)
559 <<
" >> tau[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
562 if (atleastonenegcur) {
563 innerStop_ = NEGATIVE_CURVATURE;
566 innerStop_ = EXCEEDED_TR;
575 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
576 MVT::MvScale(*this->delta_,alpha_comp);
577 MVT::MvScale(*this->Hdelta_,alpha_comp);
579 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
586 if (this->om_->isVerbosity(
Debug)) {
591 std::vector<MagnitudeType> eAx(this->blockSize_),
592 d_eAe(this->blockSize_),
593 d_eBe(this->blockSize_),
594 d_mxe(this->blockSize_);
597 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
598 TimeMonitor lcltimer( *this->timerAOp_ );
600 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
601 this->counterAOp_ += this->blockSize_;
603 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx);
606 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
607 TimeMonitor lcltimer( *this->timerAOp_ );
609 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
610 this->counterAOp_ += this->blockSize_;
612 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe);
615 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
616 TimeMonitor lcltimer( *this->timerBOp_ );
618 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
619 this->counterBOp_ += this->blockSize_;
622 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
624 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe);
627 if (this->leftMost_) {
628 for (
int j=0; j<this->blockSize_; ++j) {
629 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
633 for (
int j=0; j<this->blockSize_; ++j) {
634 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
637 this->om_->stream(
Debug)
638 <<
" Debugging checks: SIRTR inner iteration " << innerIters_ << endl
639 <<
" >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
640 for (
int j=0; j<this->blockSize_; ++j) {
641 this->om_->stream(
Debug)
642 <<
" >> m_X(eta_" << j <<
") : " << d_mxe[j] << endl;
648 if (!trncstep.empty()) {
656 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
659 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
660 TimeMonitor lcltimer( *this->timerOrtho_ );
662 this->orthman_->projectGen(
664 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
666 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
671 MagnitudeType r_norm = MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->R_,*this->R_));
679 if (this->om_->isVerbosity(
Debug)) {
680 this->om_->stream(
Debug)
681 <<
" >> |r" << innerIters_ <<
"| : " << r_norm << endl;
683 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
684 if (tconv <= kconv) {
685 innerStop_ = THETA_CONVERGENCE;
688 innerStop_ = KAPPA_CONVERGENCE;
700 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
701 TimeMonitor prectimer( *this->timerPrec_ );
703 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
704 this->counterPrec_ += this->blockSize_;
706 if (this->olsenPrec_) {
707 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
708 TimeMonitor orthtimer( *this->timerOrtho_ );
710 this->orthman_->projectGen(
712 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
714 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
717 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
718 TimeMonitor orthtimer( *this->timerOrtho_ );
720 this->orthman_->projectGen(
722 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
724 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
726 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
730 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
731 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
747 for (
unsigned int j=0; j<beta.size(); ++j) {
748 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
751 std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
752 MVT::MvScale(*this->delta_,beta_comp);
754 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
760 if (innerIters_ > d) innerIters_ = d;
762 this->om_->stream(
Debug)
763 <<
" >> stop reason is " << stopReasons_[innerStop_] << endl
769 #define SIRTR_GET_TEMP_MV(mv,workspace) \
771 TEUCHOS_TEST_FOR_EXCEPTION(workspace.size() == 0,std::logic_error,"SIRTR: Request for workspace could not be honored."); \
772 mv = workspace.back(); \
773 workspace.pop_back(); \
776 #define SIRTR_RELEASE_TEMP_MV(mv,workspace) \
778 workspace.push_back(mv); \
779 mv = Teuchos::null; \
785 template <
class ScalarType,
class MV,
class OP>
790 using Teuchos::tuple;
799 if (this->initialized_ ==
false) {
804 BB(this->blockSize_,this->blockSize_),
805 S(this->blockSize_,this->blockSize_);
812 std::vector< RCP<MV> > workspace;
815 workspace.reserve(7);
819 innerStop_ = UNINITIALIZED;
822 while (this->tester_->checkStatus(
this) !=
Passed) {
825 if (this->om_->isVerbosity(
Debug)) {
826 this->currentStatus( this->om_->stream(
Debug) );
837 totalInnerIters_ += innerIters_;
840 if (this->om_->isVerbosity(
Debug ) ) {
845 this->om_->print(
Debug, this->accuracyCheck(chk,
"in iterate() after solveTRSubproblem()") );
859 SIRTR_RELEASE_TEMP_MV(this->delta_ ,workspace);
860 SIRTR_RELEASE_TEMP_MV(this->Hdelta_,workspace);
861 SIRTR_RELEASE_TEMP_MV(this->R_ ,workspace);
862 SIRTR_RELEASE_TEMP_MV(this->Z_ ,workspace);
868 SIRTR_GET_TEMP_MV(XpEta,workspace);
870 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
871 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
873 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*XpEta);
880 MagnitudeType oldfx = this->fx_;
882 rank = this->blockSize_;
886 SIRTR_GET_TEMP_MV(AXpEta,workspace);
888 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
889 TimeMonitor lcltimer( *this->timerAOp_ );
891 OPT::Apply(*this->AOp_,*XpEta,*AXpEta);
892 this->counterAOp_ += this->blockSize_;
896 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
897 TimeMonitor lcltimer( *this->timerLocalProj_ );
899 MVT::MvTransMv(ONE,*XpEta,*AXpEta,AA);
905 SIRTR_GET_TEMP_MV(BXpEta,workspace);
907 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
908 TimeMonitor lcltimer( *this->timerBOp_ );
910 OPT::Apply(*this->BOp_,*XpEta,*BXpEta);
911 this->counterBOp_ += this->blockSize_;
915 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
916 TimeMonitor lcltimer( *this->timerLocalProj_ );
918 MVT::MvTransMv(ONE,*XpEta,*BXpEta,BB);
923 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
924 TimeMonitor lcltimer( *this->timerLocalProj_ );
926 MVT::MvTransMv(ONE,*XpEta,*XpEta,BB);
928 this->om_->stream(
Debug) <<
"AA: " << std::endl << AA << std::endl;;
929 this->om_->stream(
Debug) <<
"BB: " << std::endl << BB << std::endl;;
932 std::vector<MagnitudeType> oldtheta(this->theta_);
934 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
935 TimeMonitor lcltimer( *this->timerDS_ );
937 ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1);
939 this->om_->stream(
Debug) <<
"S: " << std::endl << S << std::endl;;
940 TEUCHOS_TEST_FOR_EXCEPTION(ret != 0,std::logic_error,
"Anasazi::SIRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret <<
"AA: " << AA << std::endl <<
"BB: " << BB << std::endl);
947 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
948 TimeMonitor lcltimer( *this->timerSort_ );
950 std::vector<int> order(this->blockSize_);
952 this->sm_->sort(this->theta_, Teuchos::rcpFromRef(order), this->blockSize_);
954 Utils::permuteVectors(order,S);
958 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
963 SIRTR_GET_TEMP_MV(AX,workspace);
964 if (this->om_->isVerbosity(
Debug ) ) {
970 MagnitudeType rhonum, rhoden, mxeta;
973 rhonum = oldfx - this->fx_;
986 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
987 TimeMonitor lcltimer( *this->timerAOp_ );
989 OPT::Apply(*this->AOp_,*this->eta_,*AX);
990 this->counterAOp_ += this->blockSize_;
996 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
997 TimeMonitor lcltimer( *this->timerBOp_ );
999 OPT::Apply(*this->BOp_,*this->eta_,*AX);
1000 this->counterBOp_ += this->blockSize_;
1004 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*AX);
1007 std::vector<MagnitudeType> eBe(this->blockSize_);
1011 std::vector<ScalarType> oldtheta_complex(oldtheta.begin(),oldtheta.end());
1012 MVT::MvScale( *AX, oldtheta_complex);
1017 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1018 TimeMonitor lcltimer( *this->timerAOp_ );
1020 OPT::Apply(*this->AOp_,*this->X_,*AX);
1021 this->counterAOp_ += this->blockSize_;
1026 mxeta = oldfx - rhoden;
1027 this->rho_ = rhonum / rhoden;
1028 this->om_->stream(
Debug)
1029 <<
" >> old f(x) is : " << oldfx << endl
1030 <<
" >> new f(x) is : " << this->fx_ << endl
1031 <<
" >> m_x(eta) is : " << mxeta << endl
1032 <<
" >> rhonum is : " << rhonum << endl
1033 <<
" >> rhoden is : " << rhoden << endl
1034 <<
" >> rho is : " << this->rho_ << endl;
1036 for (
int j=0; j<this->blockSize_; ++j) {
1037 this->om_->stream(
Debug)
1038 <<
" >> rho[" << j <<
"] : " << 1.0/(1.0+eBe[j]) << endl;
1045 this->X_ = Teuchos::null;
1046 this->BX_ = Teuchos::null;
1048 std::vector<int> ind(this->blockSize_);
1049 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
1051 X = MVT::CloneViewNonConst(*this->V_,ind);
1052 if (this->hasBOp_) {
1053 BX = MVT::CloneViewNonConst(*this->BV_,ind);
1057 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1058 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
1060 MVT::MvTimesMatAddMv(ONE,* XpEta,S,ZERO,*X);
1061 MVT::MvTimesMatAddMv(ONE,*AXpEta,S,ZERO,*AX);
1062 if (this->hasBOp_) {
1063 MVT::MvTimesMatAddMv(ONE,*BXpEta,S,ZERO,*BX);
1069 this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
1070 this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
1074 SIRTR_RELEASE_TEMP_MV(XpEta,workspace);
1075 SIRTR_RELEASE_TEMP_MV(AXpEta,workspace);
1076 if (this->hasBOp_) {
1077 SIRTR_RELEASE_TEMP_MV(BXpEta,workspace);
1085 SIRTR_GET_TEMP_MV(this->R_,workspace);
1087 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1088 TimeMonitor lcltimer( *this->timerCompRes_ );
1090 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
1092 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
1093 MVT::MvScale( *this->R_, theta_comp );
1095 MVT::MvAddMv( ONE, *AX, -ONE, *this->R_, *this->R_ );
1099 this->Rnorms_current_ =
false;
1100 this->R2norms_current_ =
false;
1104 SIRTR_RELEASE_TEMP_MV(AX,workspace);
1107 SIRTR_GET_TEMP_MV(this->delta_,workspace);
1108 SIRTR_GET_TEMP_MV(this->Hdelta_,workspace);
1109 SIRTR_GET_TEMP_MV(this->Z_,workspace);
1113 if (this->om_->isVerbosity(
Debug ) ) {
1119 this->om_->print(
Debug, this->accuracyCheck(chk,
"after local update") );
1125 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after local update") );
1135 template <
class ScalarType,
class MV,
class OP>
1140 os.setf(std::ios::scientific, std::ios::floatfield);
1143 os <<
"================================================================================" << endl;
1145 os <<
" SIRTR Solver Status" << endl;
1147 os <<
"The solver is "<<(this->initialized_ ?
"initialized." :
"not initialized.") << endl;
1148 os <<
"The number of iterations performed is " << this->iter_ << endl;
1149 os <<
"The current block size is " << this->blockSize_ << endl;
1150 os <<
"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
1151 os <<
"The number of operations A*x is " << this->counterAOp_ << endl;
1152 os <<
"The number of operations B*x is " << this->counterBOp_ << endl;
1153 os <<
"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
1154 os <<
"The number of operations Prec*x is " << this->counterPrec_ << endl;
1155 os <<
"Parameter rho_prime is " << rho_prime_ << endl;
1156 os <<
"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
1157 os <<
"Number of inner iterations was " << innerIters_ << endl;
1158 os <<
"Total number of inner iterations is " << totalInnerIters_ << endl;
1159 os <<
"f(x) is " << this->fx_ << endl;
1161 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1163 if (this->initialized_) {
1165 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1166 os << std::setw(20) <<
"Eigenvalue"
1167 << std::setw(20) <<
"Residual(B)"
1168 << std::setw(20) <<
"Residual(2)"
1170 os <<
"--------------------------------------------------------------------------------"<<endl;
1171 for (
int i=0; i<this->blockSize_; i++) {
1172 os << std::setw(20) << this->theta_[i];
1173 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
1174 else os << std::setw(20) <<
"not current";
1175 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
1176 else os << std::setw(20) <<
"not current";
1180 os <<
"================================================================================" << endl;
1187 #endif // ANASAZI_SIRTR_HPP
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.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
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...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
SIRTR(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)
SIRTR constructor with eigenproblem, solver utilities, and parameter list of solver options...
Types and exceptions used within Anasazi solvers and interfaces.
void currentStatus(std::ostream &os)
Impemements Eigensolver. This method requests that the solver print out its current status to screen...
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...
void iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
virtual ~SIRTR()
SIRTR destructor