4 #ifndef ANASAZI_SIRTR_HPP 5 #define ANASAZI_SIRTR_HPP 45 template <
class ScalarType,
class MV,
class OP>
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_;
141 template <
class ScalarType,
class MV,
class OP>
150 RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,
"SIRTR",true),
156 stopReasons_.push_back(
"n/a");
157 stopReasons_.push_back(
"maximum iterations");
158 stopReasons_.push_back(
"negative curvature");
159 stopReasons_.push_back(
"exceeded TR");
160 stopReasons_.push_back(
"kappa convergence");
161 stopReasons_.push_back(
"theta convergence");
163 rho_prime_ = params.
get(
"Rho Prime",0.5);
165 "Anasazi::SIRTR::constructor: rho_prime must be in (0,1).");
178 template <
class ScalarType,
class MV,
class OP>
189 using Teuchos::tuple;
191 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 198 innerStop_ = MAXIMUM_ITERATIONS;
200 const int n = MVT::GetGlobalLength(*this->eta_);
201 const int p = this->blockSize_;
202 const int d = n*p - (p*p+p)/2;
216 const double D2 = ONE/rho_prime_ - ONE;
218 std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p);
219 std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p);
220 MagnitudeType r0_norm;
222 MVT::MvInit(*this->eta_ ,0.0);
237 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 238 TimeMonitor lcltimer( *this->timerOrtho_ );
240 this->orthman_->projectGen(
242 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
244 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
245 if (this->leftMost_) {
246 MVT::MvScale(*this->R_,2.0);
249 MVT::MvScale(*this->R_,-2.0);
253 r0_norm = MAT::squareroot( RTRBase<ScalarType,MV,OP>::ginner(*this->R_) );
258 MagnitudeType kconv = r0_norm * this->conv_kappa_;
261 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
262 if (this->om_->isVerbosity(
Debug)) {
263 this->om_->stream(
Debug)
264 <<
" >> |r0| : " << r0_norm << endl
265 <<
" >> kappa conv : " << kconv << endl
266 <<
" >> theta conv : " << tconv << endl;
274 if (this->hasPrec_ && this->olsenPrec_)
276 std::vector<int> ind(this->blockSize_);
277 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
280 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 281 TimeMonitor prectimer( *this->timerPrec_ );
283 OPT::Apply(*this->Prec_,*this->BX_,*PBX);
284 this->counterPrec_ += this->blockSize_;
295 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 296 TimeMonitor prectimer( *this->timerPrec_ );
298 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
299 this->counterPrec_ += this->blockSize_;
301 if (this->olsenPrec_) {
302 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 303 TimeMonitor orthtimer( *this->timerOrtho_ );
305 this->orthman_->projectGen(
307 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
309 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
312 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 313 TimeMonitor orthtimer( *this->timerOrtho_ );
315 this->orthman_->projectGen(
317 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
319 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
321 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
325 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
326 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
329 if (this->om_->isVerbosity(
Debug )) {
331 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
333 if (this->hasPrec_) chk.checkZ =
true;
334 this->om_->print(
Debug, this->accuracyCheck(chk,
"after computing gradient") );
338 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
340 if (this->hasPrec_) chk.checkZ =
true;
341 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after computing gradient") );
345 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
347 if (this->om_->isVerbosity(
Debug)) {
352 std::vector<MagnitudeType> eAx(this->blockSize_),
353 d_eAe(this->blockSize_),
354 d_eBe(this->blockSize_),
355 d_mxe(this->blockSize_);
358 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 359 TimeMonitor lcltimer( *this->timerAOp_ );
361 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
362 this->counterAOp_ += this->blockSize_;
364 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx);
367 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 368 TimeMonitor lcltimer( *this->timerAOp_ );
370 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
371 this->counterAOp_ += this->blockSize_;
373 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe);
376 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 377 TimeMonitor lcltimer( *this->timerBOp_ );
379 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
380 this->counterBOp_ += this->blockSize_;
383 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
385 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe);
388 if (this->leftMost_) {
389 for (
int j=0; j<this->blockSize_; ++j) {
390 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
394 for (
int j=0; j<this->blockSize_; ++j) {
395 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
398 this->om_->stream(
Debug)
399 <<
" Debugging checks: SIRTR inner iteration " << innerIters_ << endl
400 <<
" >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
401 for (
int j=0; j<this->blockSize_; ++j) {
402 this->om_->stream(
Debug)
403 <<
" >> m_X(eta_" << j <<
") : " << d_mxe[j] << endl;
410 for (innerIters_=1; innerIters_<=d; ++innerIters_) {
418 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 419 TimeMonitor lcltimer( *this->timerAOp_ );
421 OPT::Apply(*this->AOp_,*this->delta_,*this->Z_);
422 this->counterAOp_ += this->blockSize_;
425 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 426 TimeMonitor lcltimer( *this->timerBOp_ );
428 OPT::Apply(*this->BOp_,*this->delta_,*this->Hdelta_);
429 this->counterBOp_ += this->blockSize_;
432 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Hdelta_);
436 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_ ,*this->Hdelta_,eBd);
437 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,dBd);
440 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
441 MVT::MvScale(*this->Hdelta_,theta_comp);
443 if (this->leftMost_) {
444 MVT::MvAddMv( 2.0,*this->Z_,-2.0,*this->Hdelta_,*this->Hdelta_);
447 MVT::MvAddMv(-2.0,*this->Z_, 2.0,*this->Hdelta_,*this->Hdelta_);
451 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 452 TimeMonitor lcltimer( *this->timerOrtho_ );
454 this->orthman_->projectGen(
456 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
458 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
460 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
464 for (
unsigned int j=0; j<alpha.size(); ++j)
466 alpha[j] = z_r[j]/d_Hd[j];
470 for (
unsigned int j=0; j<alpha.size(); ++j)
472 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
475 if (this->om_->isVerbosity(
Debug)) {
476 for (
unsigned int j=0; j<alpha.size(); j++) {
477 this->om_->stream(
Debug)
478 <<
" >> z_r[" << j <<
"] : " << z_r[j]
479 <<
" d_Hd[" << j <<
"] : " << d_Hd[j] << endl
480 <<
" >> eBe[" << j <<
"] : " << eBe[j]
481 <<
" neweBe[" << j <<
"] : " << new_eBe[j] << endl
482 <<
" >> eBd[" << j <<
"] : " << eBd[j]
483 <<
" dBd[" << j <<
"] : " << dBd[j] << endl;
488 std::vector<int> trncstep;
492 bool atleastonenegcur =
false;
493 for (
unsigned int j=0; j<d_Hd.size(); ++j) {
495 trncstep.push_back(j);
496 atleastonenegcur =
true;
498 else if (new_eBe[j] > D2) {
499 trncstep.push_back(j);
503 if (!trncstep.empty())
506 if (this->om_->isVerbosity(
Debug)) {
507 for (
unsigned int j=0; j<trncstep.size(); ++j) {
508 this->om_->stream(
Debug)
509 <<
" >> alpha[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
512 for (
unsigned int j=0; j<trncstep.size(); ++j) {
513 int jj = trncstep[j];
514 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
516 if (this->om_->isVerbosity(
Debug)) {
517 for (
unsigned int j=0; j<trncstep.size(); ++j) {
518 this->om_->stream(
Debug)
519 <<
" >> tau[" << trncstep[j] <<
"] : " << alpha[trncstep[j]] << endl;
522 if (atleastonenegcur) {
523 innerStop_ = NEGATIVE_CURVATURE;
526 innerStop_ = EXCEEDED_TR;
535 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
536 MVT::MvScale(*this->delta_,alpha_comp);
537 MVT::MvScale(*this->Hdelta_,alpha_comp);
539 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
546 if (this->om_->isVerbosity(
Debug)) {
551 std::vector<MagnitudeType> eAx(this->blockSize_),
552 d_eAe(this->blockSize_),
553 d_eBe(this->blockSize_),
554 d_mxe(this->blockSize_);
557 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 558 TimeMonitor lcltimer( *this->timerAOp_ );
560 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
561 this->counterAOp_ += this->blockSize_;
563 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx);
566 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 567 TimeMonitor lcltimer( *this->timerAOp_ );
569 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
570 this->counterAOp_ += this->blockSize_;
572 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe);
575 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 576 TimeMonitor lcltimer( *this->timerBOp_ );
578 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
579 this->counterBOp_ += this->blockSize_;
582 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
584 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe);
587 if (this->leftMost_) {
588 for (
int j=0; j<this->blockSize_; ++j) {
589 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
593 for (
int j=0; j<this->blockSize_; ++j) {
594 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
597 this->om_->stream(
Debug)
598 <<
" Debugging checks: SIRTR inner iteration " << innerIters_ << endl
599 <<
" >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
600 for (
int j=0; j<this->blockSize_; ++j) {
601 this->om_->stream(
Debug)
602 <<
" >> m_X(eta_" << j <<
") : " << d_mxe[j] << endl;
608 if (!trncstep.empty()) {
616 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
619 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 620 TimeMonitor lcltimer( *this->timerOrtho_ );
622 this->orthman_->projectGen(
624 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
626 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
631 MagnitudeType r_norm = MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->R_,*this->R_));
639 if (this->om_->isVerbosity(
Debug)) {
640 this->om_->stream(
Debug)
641 <<
" >> |r" << innerIters_ <<
"| : " << r_norm << endl;
643 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
644 if (tconv <= kconv) {
645 innerStop_ = THETA_CONVERGENCE;
648 innerStop_ = KAPPA_CONVERGENCE;
660 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 661 TimeMonitor prectimer( *this->timerPrec_ );
663 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
664 this->counterPrec_ += this->blockSize_;
666 if (this->olsenPrec_) {
667 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 668 TimeMonitor orthtimer( *this->timerOrtho_ );
670 this->orthman_->projectGen(
672 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),
false,
674 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
677 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 678 TimeMonitor orthtimer( *this->timerOrtho_ );
680 this->orthman_->projectGen(
682 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),
false,
684 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
686 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
690 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
691 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
707 for (
unsigned int j=0; j<beta.size(); ++j) {
708 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
711 std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
712 MVT::MvScale(*this->delta_,beta_comp);
714 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
720 if (innerIters_ > d) innerIters_ = d;
722 this->om_->stream(
Debug)
723 <<
" >> stop reason is " << stopReasons_[innerStop_] << endl
729 #define SIRTR_GET_TEMP_MV(mv,workspace) \ 731 TEUCHOS_TEST_FOR_EXCEPTION(workspace.size() == 0,std::logic_error,"SIRTR: Request for workspace could not be honored."); \ 732 mv = workspace.back(); \ 733 workspace.pop_back(); \ 736 #define SIRTR_RELEASE_TEMP_MV(mv,workspace) \ 738 workspace.push_back(mv); \ 739 mv = Teuchos::null; \ 745 template <
class ScalarType,
class MV,
class OP>
750 using Teuchos::tuple;
759 if (this->initialized_ ==
false) {
764 BB(this->blockSize_,this->blockSize_),
765 S(this->blockSize_,this->blockSize_);
772 std::vector< RCP<MV> > workspace;
775 workspace.reserve(7);
779 innerStop_ = UNINITIALIZED;
782 while (this->tester_->checkStatus(
this) !=
Passed) {
785 if (this->om_->isVerbosity(
Debug)) {
786 this->currentStatus( this->om_->stream(
Debug) );
797 totalInnerIters_ += innerIters_;
800 if (this->om_->isVerbosity(
Debug ) ) {
805 this->om_->print(
Debug, this->accuracyCheck(chk,
"in iterate() after solveTRSubproblem()") );
819 SIRTR_RELEASE_TEMP_MV(this->delta_ ,workspace);
820 SIRTR_RELEASE_TEMP_MV(this->Hdelta_,workspace);
821 SIRTR_RELEASE_TEMP_MV(this->R_ ,workspace);
822 SIRTR_RELEASE_TEMP_MV(this->Z_ ,workspace);
828 SIRTR_GET_TEMP_MV(XpEta,workspace);
830 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 831 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
833 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*XpEta);
840 MagnitudeType oldfx = this->fx_;
842 rank = this->blockSize_;
846 SIRTR_GET_TEMP_MV(AXpEta,workspace);
848 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 849 TimeMonitor lcltimer( *this->timerAOp_ );
851 OPT::Apply(*this->AOp_,*XpEta,*AXpEta);
852 this->counterAOp_ += this->blockSize_;
856 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 857 TimeMonitor lcltimer( *this->timerLocalProj_ );
859 MVT::MvTransMv(ONE,*XpEta,*AXpEta,AA);
865 SIRTR_GET_TEMP_MV(BXpEta,workspace);
867 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 868 TimeMonitor lcltimer( *this->timerBOp_ );
870 OPT::Apply(*this->BOp_,*XpEta,*BXpEta);
871 this->counterBOp_ += this->blockSize_;
875 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 876 TimeMonitor lcltimer( *this->timerLocalProj_ );
878 MVT::MvTransMv(ONE,*XpEta,*BXpEta,BB);
883 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 884 TimeMonitor lcltimer( *this->timerLocalProj_ );
886 MVT::MvTransMv(ONE,*XpEta,*XpEta,BB);
888 this->om_->stream(
Debug) <<
"AA: " << std::endl << AA << std::endl;;
889 this->om_->stream(
Debug) <<
"BB: " << std::endl << BB << std::endl;;
892 std::vector<MagnitudeType> oldtheta(this->theta_);
894 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 895 TimeMonitor lcltimer( *this->timerDS_ );
897 ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1);
899 this->om_->stream(
Debug) <<
"S: " << std::endl << S << std::endl;;
900 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);
907 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 908 TimeMonitor lcltimer( *this->timerSort_ );
910 std::vector<int> order(this->blockSize_);
912 this->sm_->sort(this->theta_, Teuchos::rcpFromRef(order), this->blockSize_);
914 Utils::permuteVectors(order,S);
918 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
923 SIRTR_GET_TEMP_MV(AX,workspace);
924 if (this->om_->isVerbosity(
Debug ) ) {
930 MagnitudeType rhonum, rhoden, mxeta;
933 rhonum = oldfx - this->fx_;
946 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 947 TimeMonitor lcltimer( *this->timerAOp_ );
949 OPT::Apply(*this->AOp_,*this->eta_,*AX);
950 this->counterAOp_ += this->blockSize_;
956 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 957 TimeMonitor lcltimer( *this->timerBOp_ );
959 OPT::Apply(*this->BOp_,*this->eta_,*AX);
960 this->counterBOp_ += this->blockSize_;
964 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*AX);
967 std::vector<MagnitudeType> eBe(this->blockSize_);
971 std::vector<ScalarType> oldtheta_complex(oldtheta.begin(),oldtheta.end());
972 MVT::MvScale( *AX, oldtheta_complex);
977 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 978 TimeMonitor lcltimer( *this->timerAOp_ );
980 OPT::Apply(*this->AOp_,*this->X_,*AX);
981 this->counterAOp_ += this->blockSize_;
986 mxeta = oldfx - rhoden;
987 this->rho_ = rhonum / rhoden;
988 this->om_->stream(
Debug)
989 <<
" >> old f(x) is : " << oldfx << endl
990 <<
" >> new f(x) is : " << this->fx_ << endl
991 <<
" >> m_x(eta) is : " << mxeta << endl
992 <<
" >> rhonum is : " << rhonum << endl
993 <<
" >> rhoden is : " << rhoden << endl
994 <<
" >> rho is : " << this->rho_ << endl;
996 for (
int j=0; j<this->blockSize_; ++j) {
997 this->om_->stream(
Debug)
998 <<
" >> rho[" << j <<
"] : " << 1.0/(1.0+eBe[j]) << endl;
1005 this->X_ = Teuchos::null;
1006 this->BX_ = Teuchos::null;
1008 std::vector<int> ind(this->blockSize_);
1009 for (
int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
1011 X = MVT::CloneViewNonConst(*this->V_,ind);
1012 if (this->hasBOp_) {
1013 BX = MVT::CloneViewNonConst(*this->BV_,ind);
1017 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1018 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
1020 MVT::MvTimesMatAddMv(ONE,* XpEta,S,ZERO,*X);
1021 MVT::MvTimesMatAddMv(ONE,*AXpEta,S,ZERO,*AX);
1022 if (this->hasBOp_) {
1023 MVT::MvTimesMatAddMv(ONE,*BXpEta,S,ZERO,*BX);
1029 this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
1030 this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
1034 SIRTR_RELEASE_TEMP_MV(XpEta,workspace);
1035 SIRTR_RELEASE_TEMP_MV(AXpEta,workspace);
1036 if (this->hasBOp_) {
1037 SIRTR_RELEASE_TEMP_MV(BXpEta,workspace);
1045 SIRTR_GET_TEMP_MV(this->R_,workspace);
1047 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1048 TimeMonitor lcltimer( *this->timerCompRes_ );
1050 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
1052 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
1053 MVT::MvScale( *this->R_, theta_comp );
1055 MVT::MvAddMv( ONE, *AX, -ONE, *this->R_, *this->R_ );
1059 this->Rnorms_current_ =
false;
1060 this->R2norms_current_ =
false;
1064 SIRTR_RELEASE_TEMP_MV(AX,workspace);
1067 SIRTR_GET_TEMP_MV(this->delta_,workspace);
1068 SIRTR_GET_TEMP_MV(this->Hdelta_,workspace);
1069 SIRTR_GET_TEMP_MV(this->Z_,workspace);
1073 if (this->om_->isVerbosity(
Debug ) ) {
1079 this->om_->print(
Debug, this->accuracyCheck(chk,
"after local update") );
1085 this->om_->print(
OrthoDetails, this->accuracyCheck(chk,
"after local update") );
1095 template <
class ScalarType,
class MV,
class OP>
1100 os.setf(std::ios::scientific, std::ios::floatfield);
1103 os <<
"================================================================================" << endl;
1105 os <<
" SIRTR Solver Status" << endl;
1107 os <<
"The solver is "<<(this->initialized_ ?
"initialized." :
"not initialized.") << endl;
1108 os <<
"The number of iterations performed is " << this->iter_ << endl;
1109 os <<
"The current block size is " << this->blockSize_ << endl;
1110 os <<
"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
1111 os <<
"The number of operations A*x is " << this->counterAOp_ << endl;
1112 os <<
"The number of operations B*x is " << this->counterBOp_ << endl;
1113 os <<
"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
1114 os <<
"The number of operations Prec*x is " << this->counterPrec_ << endl;
1115 os <<
"Parameter rho_prime is " << rho_prime_ << endl;
1116 os <<
"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
1117 os <<
"Number of inner iterations was " << innerIters_ << endl;
1118 os <<
"Total number of inner iterations is " << totalInnerIters_ << endl;
1119 os <<
"f(x) is " << this->fx_ << endl;
1121 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1123 if (this->initialized_) {
1125 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1126 os << std::setw(20) <<
"Eigenvalue" 1127 << std::setw(20) <<
"Residual(B)" 1128 << std::setw(20) <<
"Residual(2)" 1130 os <<
"--------------------------------------------------------------------------------"<<endl;
1131 for (
int i=0; i<this->blockSize_; i++) {
1132 os << std::setw(20) << this->theta_[i];
1133 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
1134 else os << std::setw(20) <<
"not current";
1135 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
1136 else os << std::setw(20) <<
"not current";
1140 os <<
"================================================================================" << endl;
1147 #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.
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...
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