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;
205 const int n = MVT::GetGlobalLength(*this->eta_);
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;
227 MVT::MvInit(*this->eta_ ,0.0);
228 MVT::MvInit(*this->Aeta_,0.0);
230 MVT::MvInit(*this->Beta_,0.0);
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_) {
255 MVT::MvScale(*this->R_,2.0);
258 MVT::MvScale(*this->R_,-2.0);
262 r0_norm = MAT::squareroot( RTRBase<ScalarType,MV,OP>::ginner(*this->R_) );
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;
287 Teuchos::RCP<MV> PBX = MVT::CloneViewNonConst(*this->PBV_,ind);
289 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 290 TimeMonitor prectimer( *this->timerPrec_ );
292 OPT::Apply(*this->Prec_,*this->BX_,*PBX);
293 this->counterPrec_ += this->blockSize_;
304 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 305 TimeMonitor prectimer( *this->timerPrec_ );
307 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
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_));
330 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
333 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
336 if (this->om_->isVerbosity(
Debug )) {
338 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
340 if (this->hasPrec_) chk.checkZ =
true;
341 this->om_->print(
Debug, this->accuracyCheck(chk,
"after computing gradient") );
345 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
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_);
357 MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
358 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
359 MVT::MvScale(*Heta,theta_comp);
360 if (this->leftMost_) {
361 MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta);
364 MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta);
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_);
380 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->AX_,eg);
381 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*Heta,eHe);
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_);
430 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_ ,*this->Bdelta_,eBd);
431 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Bdelta_,dBd);
433 MVT::MvAddMv(ONE,*this->Bdelta_,ZERO,*this->Bdelta_,*this->Hdelta_);
435 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
436 MVT::MvScale(*this->Hdelta_,theta_comp);
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_));
455 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
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());
531 MVT::MvScale(*this->delta_,alpha_comp);
532 MVT::MvScale(*this->Adelta_,alpha_comp);
534 MVT::MvScale(*this->Bdelta_,alpha_comp);
536 MVT::MvScale(*this->Hdelta_,alpha_comp);
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_);
552 MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
554 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
555 MVT::MvScale(*Heta,theta_comp);
557 if (this->leftMost_) {
558 MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta);
561 MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta);
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_);
577 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->AX_,eg);
578 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*Heta,eHe);
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_));
623 MagnitudeType r_norm = MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->R_,*this->R_));
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_ );
655 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
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_));
678 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
681 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
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());
702 MVT::MvScale(*this->delta_,beta_comp);
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)) {
762 this->currentStatus( this->om_->stream(
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_);
809 MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,XE);
810 MVT::MvTransMv(ONE,*this->eta_,*this->Beta_,E);
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_);
838 MVT::MvTransMv(ONE,*this->X_ ,*this->AX_ ,AA11);
839 MVT::MvTransMv(ONE,*this->X_ ,*this->Aeta_,AA12);
840 MVT::MvTransMv(ONE,*this->eta_,*this->Aeta_,AA22);
841 MVT::MvTransMv(ONE,*this->X_ ,*this->BX_ ,BB11);
842 MVT::MvTransMv(ONE,*this->X_ ,*this->Beta_,BB12);
843 MVT::MvTransMv(ONE,*this->eta_,*this->Beta_,BB22);
846 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 847 TimeMonitor lcltimer( *this->timerLocalProj_ );
849 MVT::MvTransMv(ONE,*this->delta_,*this->Adelta_,AA);
850 MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,BB);
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_ );
858 ret = Utils::directSolver(AA.numRows(),AA,Teuchos::rcpFromRef(BB),S,newtheta,rank,1);
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);
875 Utils::permuteVectors(order,S);
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;
938 X = MVT::CloneViewNonConst(*this->V_,ind);
940 BX = MVT::CloneViewNonConst(*this->BV_,ind);
942 if (useSA_ ==
false) {
946 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 947 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
949 MVT::MvTimesMatAddMv(ONE,*this->delta_,S,ZERO,*X);
950 MVT::MvTimesMatAddMv(ONE,*this->Adelta_,S,ZERO,*this->AX_);
952 MVT::MvTimesMatAddMv(ONE,*this->Bdelta_,S,ZERO,*BX);
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_ );
966 MVT::MvTimesMatAddMv(ONE,* X ,Sx,ZERO,*this->delta_);
967 MVT::MvTimesMatAddMv(ONE,*this->eta_,Se,ONE ,*this->delta_);
968 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*X);
970 MVT::MvTimesMatAddMv(ONE,*this->AX_ ,Sx,ZERO,*this->delta_);
971 MVT::MvTimesMatAddMv(ONE,*this->Aeta_,Se,ONE ,*this->delta_);
972 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->AX_);
975 MVT::MvTimesMatAddMv(ONE,* BX ,Sx,ZERO,*this->delta_);
976 MVT::MvTimesMatAddMv(ONE,*this->Beta_,Se,ONE ,*this->delta_);
977 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*BX);
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());
995 MVT::MvScale( *this->R_, theta_comp );
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 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.
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...
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.
Virtual base class which defines basic traits for the operator type.
void iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
Types and exceptions used within Anasazi solvers and interfaces.
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...
void currentStatus(std::ostream &os)
Impemements Eigensolver. This method requests that the solver print out its current status to screen...