00001
00002
00003
00004 #ifndef ANASAZI_IRTR_HPP
00005 #define ANASAZI_IRTR_HPP
00006
00007 #include "AnasaziTypes.hpp"
00008 #include "AnasaziRTRBase.hpp"
00009
00010 #include "AnasaziEigensolver.hpp"
00011 #include "AnasaziMultiVecTraits.hpp"
00012 #include "AnasaziOperatorTraits.hpp"
00013 #include "Teuchos_ScalarTraits.hpp"
00014
00015 #include "Teuchos_LAPACK.hpp"
00016 #include "Teuchos_BLAS.hpp"
00017 #include "Teuchos_SerialDenseMatrix.hpp"
00018 #include "Teuchos_ParameterList.hpp"
00019 #include "Teuchos_TimeMonitor.hpp"
00020
00040
00041
00042
00043 namespace Anasazi {
00044
00045 template <class ScalarType, class MV, class OP>
00046 class IRTR : public RTRBase<ScalarType,MV,OP> {
00047 public:
00048
00050
00051
00063 IRTR( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00064 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00065 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00066 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00067 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00068 Teuchos::ParameterList ¶ms
00069 );
00070
00072 virtual ~IRTR() {};
00073
00075
00077
00078
00080 void iterate();
00081
00083
00085
00086
00088 void currentStatus(std::ostream &os);
00089
00091
00092 private:
00093
00094
00095
00096 typedef SolverUtils<ScalarType,MV,OP> Utils;
00097 typedef MultiVecTraits<ScalarType,MV> MVT;
00098 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00099 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00100 typedef typename SCT::magnitudeType MagnitudeType;
00101 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
00102 enum trRetType {
00103 UNINITIALIZED = 0,
00104 MAXIMUM_ITERATIONS,
00105 NEGATIVE_CURVATURE,
00106 EXCEEDED_TR,
00107 KAPPA_CONVERGENCE,
00108 THETA_CONVERGENCE
00109 };
00110
00111 std::vector<std::string> stopReasons_;
00112
00113
00114
00115 const MagnitudeType ZERO;
00116 const MagnitudeType ONE;
00117
00118
00119
00121 void solveTRSubproblem();
00122
00123
00124 MagnitudeType rho_prime_;
00125
00126
00127 MagnitudeType normgradf0_;
00128
00129
00130 trRetType innerStop_;
00131
00132
00133 int innerIters_;
00134 };
00135
00136
00137
00138
00140
00141 template <class ScalarType, class MV, class OP>
00142 IRTR<ScalarType,MV,OP>::IRTR(
00143 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00144 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00145 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00146 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00147 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00148 Teuchos::ParameterList ¶ms
00149 ) :
00150 RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,"IRTR",false),
00151 ZERO(MAT::zero()),
00152 ONE(MAT::one())
00153 {
00154
00155 stopReasons_.push_back("n/a");
00156 stopReasons_.push_back("maximum iterations");
00157 stopReasons_.push_back("negative curvature");
00158 stopReasons_.push_back("exceeded TR");
00159 stopReasons_.push_back("kappa convergence");
00160 stopReasons_.push_back("theta convergence");
00161
00162 rho_prime_ = params.get("Rho Prime",0.5);
00163 TEST_FOR_EXCEPTION(rho_prime_ <= 0 || rho_prime_ >= 1,std::invalid_argument,
00164 "Anasazi::IRTR::constructor: rho_prime must be in (0,1).");
00165 }
00166
00167
00169
00170
00171
00172
00173
00174
00175
00176
00177 template <class ScalarType, class MV, class OP>
00178 void IRTR<ScalarType,MV,OP>::solveTRSubproblem() {
00179
00180
00181
00182
00183
00184
00185
00186
00187 using Teuchos::RCP;
00188 using Teuchos::tuple;
00189 using Teuchos::null;
00190 using Teuchos::TimeMonitor;
00191 using std::endl;
00192 typedef Teuchos::RCP<const MV> PCMV;
00193 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
00194
00195 innerStop_ = MAXIMUM_ITERATIONS;
00196
00197 const int n = MVT::GetVecLength(*this->eta_);
00198 const int p = this->blockSize_;
00199 const int d = n*p - (p*p+p)/2;
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 const double D2 = ONE/rho_prime_ - ONE;
00214
00215 std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p);
00216 std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p);
00217 MagnitudeType r0_norm;
00218
00219 MVT::MvInit(*this->eta_ ,0.0);
00220 MVT::MvInit(*this->Aeta_,0.0);
00221 if (this->hasBOp_) {
00222 MVT::MvInit(*this->Beta_,0.0);
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 {
00238 TimeMonitor lcltimer( *this->timerOrtho_ );
00239 this->orthman_->projectGen(
00240 *this->R_,
00241 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00242 tuple<PSDM>(null),
00243 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00244 if (this->leftMost_) {
00245 MVT::MvScale(*this->R_,2.0);
00246 }
00247 else {
00248 MVT::MvScale(*this->R_,-2.0);
00249 }
00250 }
00251
00252 r0_norm = MAT::squareroot( ginner(*this->R_) );
00253
00254
00255
00256
00257 MagnitudeType kconv = r0_norm * this->conv_kappa_;
00258
00259
00260 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
00261 if (this->om_->isVerbosity(Debug)) {
00262 this->om_->stream(Debug)
00263 << " >> |r0| : " << r0_norm << endl
00264 << " >> kappa conv : " << kconv << endl
00265 << " >> theta conv : " << tconv << endl;
00266 }
00267
00268
00269
00270
00271
00272 if (this->hasPrec_)
00273 {
00274 std::vector<int> ind(this->blockSize_);
00275 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
00276 Teuchos::RCP<MV> PBX = MVT::CloneView(*this->PBV_,ind);
00277 {
00278 TimeMonitor prectimer( *this->timerPrec_ );
00279 OPT::Apply(*this->Prec_,*this->BX_,*PBX);
00280 this->counterPrec_ += this->blockSize_;
00281 }
00282 PBX = Teuchos::null;
00283 }
00284
00285
00286
00287 if (this->hasPrec_)
00288 {
00289 TimeMonitor prectimer( *this->timerPrec_ );
00290 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
00291 this->counterPrec_ += this->blockSize_;
00292
00293 {
00294 TimeMonitor orthtimer( *this->timerOrtho_ );
00295 this->orthman_->projectGen(
00296 *this->Z_,
00297 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false,
00298 tuple<PSDM>(null),
00299 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00300 }
00301 ginnersep(*this->R_,*this->Z_,z_r);
00302 }
00303 else {
00304 ginnersep(*this->R_,z_r);
00305 }
00306
00307 if (this->om_->isVerbosity( Debug )) {
00308
00309 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00310 chk.checkBR = true;
00311 if (this->hasPrec_) chk.checkZ = true;
00312 this->om_->print( Debug, this->accuracyCheck(chk, "after computing gradient") );
00313 }
00314 else if (this->om_->isVerbosity( OrthoDetails )) {
00315
00316 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00317 chk.checkBR = true;
00318 if (this->hasPrec_) chk.checkZ = true;
00319 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after computing gradient") );
00320 }
00321
00322
00323 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
00324
00325 if (this->om_->isVerbosity(Debug)) {
00326 RCP<MV> Heta = MVT::Clone(*this->eta_,this->blockSize_);
00327
00328 MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
00329 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00330 MVT::MvScale(*Heta,theta_comp);
00331 if (this->leftMost_) {
00332 MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta);
00333 }
00334 else {
00335 MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta);
00336 }
00337
00338 {
00339 TimeMonitor lcltimer( *this->timerOrtho_ );
00340 this->orthman_->projectGen(
00341 *Heta,
00342 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00343 tuple<PSDM>(null),
00344 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00345 }
00346
00347 std::vector<MagnitudeType> eg(this->blockSize_),
00348 eHe(this->blockSize_);
00349 ginnersep(*this->eta_,*this->AX_,eg);
00350 ginnersep(*this->eta_,*Heta,eHe);
00351 if (this->leftMost_) {
00352 for (int j=0; j<this->blockSize_; ++j) {
00353 eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
00354 }
00355 }
00356 else {
00357 for (int j=0; j<this->blockSize_; ++j) {
00358 eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
00359 }
00360 }
00361 this->om_->stream(Debug)
00362 << " Debugging checks: IRTR inner iteration " << innerIters_ << endl
00363 << " >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
00364 for (int j=0; j<this->blockSize_; ++j) {
00365 this->om_->stream(Debug)
00366 << " >> m_X(eta_" << j << ") : " << eg[j] << endl;
00367 }
00368 }
00369
00372
00373 for (innerIters_=1; innerIters_<=d; ++innerIters_) {
00374
00375
00376
00377
00378
00379
00380 {
00381 TimeMonitor lcltimer( *this->timerAOp_ );
00382 OPT::Apply(*this->AOp_,*this->delta_,*this->Adelta_);
00383 this->counterAOp_ += this->blockSize_;
00384 }
00385 if (this->hasBOp_) {
00386 TimeMonitor lcltimer( *this->timerBOp_ );
00387 OPT::Apply(*this->BOp_,*this->delta_,*this->Bdelta_);
00388 this->counterBOp_ += this->blockSize_;
00389 }
00390 else {
00391 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Bdelta_);
00392 }
00393
00394
00395 ginnersep(*this->eta_ ,*this->Bdelta_,eBd);
00396 ginnersep(*this->delta_,*this->Bdelta_,dBd);
00397
00398 MVT::MvAddMv(ONE,*this->Bdelta_,ZERO,*this->Bdelta_,*this->Hdelta_);
00399 {
00400 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00401 MVT::MvScale(*this->Hdelta_,theta_comp);
00402 }
00403 if (this->leftMost_) {
00404 MVT::MvAddMv( 2.0,*this->Adelta_,-2.0,*this->Hdelta_,*this->Hdelta_);
00405 }
00406 else {
00407 MVT::MvAddMv(-2.0,*this->Adelta_, 2.0,*this->Hdelta_,*this->Hdelta_);
00408 }
00409
00410 {
00411 TimeMonitor lcltimer( *this->timerOrtho_ );
00412 this->orthman_->projectGen(
00413 *this->Hdelta_,
00414 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00415 tuple<PSDM>(null),
00416 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00417 }
00418 ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
00419
00420
00421
00422 for (unsigned int j=0; j<alpha.size(); ++j)
00423 {
00424 alpha[j] = z_r[j]/d_Hd[j];
00425 }
00426
00427
00428 for (unsigned int j=0; j<alpha.size(); ++j)
00429 {
00430 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
00431 }
00432
00433 if (this->om_->isVerbosity(Debug)) {
00434 for (unsigned int j=0; j<alpha.size(); j++) {
00435 this->om_->stream(Debug)
00436 << " >> z_r[" << j << "] : " << z_r[j]
00437 << " d_Hd[" << j << "] : " << d_Hd[j] << endl
00438 << " >> eBe[" << j << "] : " << eBe[j]
00439 << " neweBe[" << j << "] : " << new_eBe[j] << endl
00440 << " >> eBd[" << j << "] : " << eBd[j]
00441 << " dBd[" << j << "] : " << dBd[j] << endl;
00442 }
00443 }
00444
00445
00446 std::vector<int> trncstep;
00447 trncstep.reserve(p);
00448
00449
00450 bool atleastonenegcur = false;
00451 for (unsigned int j=0; j<d_Hd.size(); ++j) {
00452 if (d_Hd[j] <= 0) {
00453 trncstep.push_back(j);
00454 atleastonenegcur = true;
00455 }
00456 else if (new_eBe[j] > D2) {
00457 trncstep.push_back(j);
00458 }
00459 }
00460
00461 if (!trncstep.empty())
00462 {
00463
00464 if (this->om_->isVerbosity(Debug)) {
00465 for (unsigned int j=0; j<trncstep.size(); ++j) {
00466 this->om_->stream(Debug)
00467 << " >> alpha[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
00468 }
00469 }
00470 for (unsigned int j=0; j<trncstep.size(); ++j) {
00471 int jj = trncstep[j];
00472 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
00473 }
00474 if (this->om_->isVerbosity(Debug)) {
00475 for (unsigned int j=0; j<trncstep.size(); ++j) {
00476 this->om_->stream(Debug)
00477 << " >> tau[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
00478 }
00479 }
00480 if (atleastonenegcur) {
00481 innerStop_ = NEGATIVE_CURVATURE;
00482 }
00483 else {
00484 innerStop_ = EXCEEDED_TR;
00485 }
00486 }
00487
00488
00489
00490
00491
00492 {
00493 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
00494 MVT::MvScale(*this->delta_,alpha_comp);
00495 MVT::MvScale(*this->Adelta_,alpha_comp);
00496 if (this->hasBOp_) {
00497 MVT::MvScale(*this->Bdelta_,alpha_comp);
00498 }
00499 MVT::MvScale(*this->Hdelta_,alpha_comp);
00500 }
00501 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
00502 MVT::MvAddMv(ONE,*this->Adelta_,ONE,*this->Aeta_,*this->Aeta_);
00503 if (this->hasBOp_) {
00504 MVT::MvAddMv(ONE,*this->Bdelta_,ONE,*this->Beta_,*this->Beta_);
00505 }
00506
00507
00508 eBe = new_eBe;
00509
00510
00511
00512 if (this->om_->isVerbosity(Debug)) {
00513 RCP<MV> Heta = MVT::Clone(*this->eta_,this->blockSize_);
00514
00515 MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
00516 {
00517 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00518 MVT::MvScale(*Heta,theta_comp);
00519 }
00520 if (this->leftMost_) {
00521 MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta);
00522 }
00523 else {
00524 MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta);
00525 }
00526
00527 {
00528 TimeMonitor lcltimer( *this->timerOrtho_ );
00529 this->orthman_->projectGen(
00530 *Heta,
00531 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00532 tuple<PSDM>(null),
00533 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00534 }
00535
00536 std::vector<MagnitudeType> eg(this->blockSize_),
00537 eHe(this->blockSize_);
00538 ginnersep(*this->eta_,*this->AX_,eg);
00539 ginnersep(*this->eta_,*Heta,eHe);
00540 if (this->leftMost_) {
00541 for (int j=0; j<this->blockSize_; ++j) {
00542 eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
00543 }
00544 }
00545 else {
00546 for (int j=0; j<this->blockSize_; ++j) {
00547 eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
00548 }
00549 }
00550 this->om_->stream(Debug)
00551 << " Debugging checks: IRTR inner iteration " << innerIters_ << endl
00552 << " >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
00553 for (int j=0; j<this->blockSize_; ++j) {
00554 this->om_->stream(Debug)
00555 << " >> m_X(eta_" << j << ") : " << eg[j] << endl;
00556 }
00557 }
00558
00559
00560
00561 if (!trncstep.empty()) {
00562 break;
00563 }
00564
00565
00566
00567
00568
00569 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
00570 {
00571
00572 TimeMonitor lcltimer( *this->timerOrtho_ );
00573 this->orthman_->projectGen(
00574 *this->R_,
00575 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00576 tuple<PSDM>(null),
00577 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00578 }
00579
00580
00581
00582 MagnitudeType r_norm = MAT::squareroot(ginner(*this->R_,*this->R_));
00583
00584
00585
00586
00587
00588
00589
00590 if (this->om_->isVerbosity(Debug)) {
00591 this->om_->stream(Debug)
00592 << " >> |r" << innerIters_ << "| : " << r_norm << endl;
00593 }
00594 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
00595 if (tconv <= kconv) {
00596 innerStop_ = THETA_CONVERGENCE;
00597 }
00598 else {
00599 innerStop_ = KAPPA_CONVERGENCE;
00600 }
00601 break;
00602 }
00603
00604
00605
00606 zold_rold = z_r;
00607 if (this->hasPrec_)
00608 {
00609 TimeMonitor prectimer( *this->timerPrec_ );
00610 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
00611 this->counterPrec_ += this->blockSize_;
00612
00613 {
00614 TimeMonitor orthtimer( *this->timerOrtho_ );
00615 this->orthman_->projectGen(
00616 *this->Z_,
00617 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false,
00618 tuple<PSDM>(null),
00619 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00620 }
00621 ginnersep(*this->R_,*this->Z_,z_r);
00622 }
00623 else {
00624 ginnersep(*this->R_,z_r);
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640 for (unsigned int j=0; j<beta.size(); ++j) {
00641 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
00642 }
00643 {
00644 std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
00645 MVT::MvScale(*this->delta_,beta_comp);
00646 }
00647 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
00648
00649 }
00650
00653 if (innerIters_ > d) innerIters_ = d;
00654
00655 this->om_->stream(Debug)
00656 << " >> stop reason is " << stopReasons_[innerStop_] << endl
00657 << endl;
00658
00659 }
00660
00661
00663
00664 template <class ScalarType, class MV, class OP>
00665 void IRTR<ScalarType,MV,OP>::iterate() {
00666
00667 using Teuchos::RCP;
00668 using Teuchos::null;
00669 using Teuchos::tuple;
00670 using Teuchos::TimeMonitor;
00671 using std::endl;
00672 typedef Teuchos::RCP<const MV> PCMV;
00673 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
00674
00675
00676
00677
00678 if (this->initialized_ == false) {
00679 this->initialize();
00680 }
00681
00682 Teuchos::SerialDenseMatrix<int,ScalarType> AA(this->blockSize_,this->blockSize_),
00683 BB(this->blockSize_,this->blockSize_),
00684 S(this->blockSize_,this->blockSize_);
00685
00686
00687 innerIters_ = -1;
00688 innerStop_ = UNINITIALIZED;
00689
00690
00691 while (this->tester_->checkStatus(this) != Passed) {
00692
00693
00694 if (this->om_->isVerbosity(Debug)) {
00695 this->currentStatus( this->om_->stream(Debug) );
00696 }
00697 else if (this->om_->isVerbosity(IterationDetails)) {
00698 this->currentStatus( this->om_->stream(IterationDetails) );
00699 }
00700
00701
00702 this->iter_++;
00703
00704
00705 solveTRSubproblem();
00706
00707
00708 if (this->om_->isVerbosity( Debug ) ) {
00709 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00710
00711 chk.checkBR = true;
00712 chk.checkEta = true;
00713 chk.checkAEta = true;
00714 chk.checkBEta = true;
00715 this->om_->print( Debug, this->accuracyCheck(chk, "in iterate() after solveTRSubproblem()") );
00716 this->om_->stream(Debug)
00717 << " >> norm(Eta) : " << MAT::squareroot(ginner(*this->eta_)) << endl
00718 << endl;
00719 }
00720
00721
00722
00723 {
00724 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
00725 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*this->delta_);
00726 MVT::MvAddMv(ONE,*this->AX_,ONE,*this->Aeta_,*this->Adelta_);
00727 if (this->hasBOp_) {
00728 MVT::MvAddMv(ONE,*this->BX_,ONE,*this->Beta_,*this->Bdelta_);
00729 }
00730 }
00731
00732
00733 if (this->om_->isVerbosity( Debug ) ) {
00734
00735
00736
00737 Teuchos::SerialDenseMatrix<int,ScalarType> XE(this->blockSize_,this->blockSize_),
00738 E(this->blockSize_,this->blockSize_);
00739 MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,XE);
00740 MVT::MvTransMv(ONE,*this->eta_,*this->Beta_,E);
00741 this->om_->stream(Debug)
00742 << " >> Error in AX+AEta == A(X+Eta) : " << Utils::errorEquality(*this->delta_,*this->Adelta_,this->AOp_) << endl
00743 << " >> Error in BX+BEta == B(X+Eta) : " << Utils::errorEquality(*this->delta_,*this->Bdelta_,this->BOp_) << endl
00744 << " >> norm( (X+Eta)^T B (X+Eta) ) : " << XE.normFrobenius() << endl
00745 << " >> norm( Eta^T B Eta ) : " << E.normFrobenius() << endl
00746 << endl;
00747 }
00748
00749
00750
00751
00752
00753 MagnitudeType oldfx = this->fx_;
00754 std::vector<MagnitudeType> oldtheta(this->theta_);
00755 int rank, ret;
00756 rank = this->blockSize_;
00757 {
00758 TimeMonitor lcltimer( *this->timerLocalProj_ );
00759 MVT::MvTransMv(ONE,*this->delta_,*this->Adelta_,AA);
00760 MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,BB);
00761 }
00762 this->om_->stream(Debug) << "AA: " << std::endl << AA << std::endl;;
00763 this->om_->stream(Debug) << "BB: " << std::endl << BB << std::endl;;
00764 {
00765 TimeMonitor lcltimer( *this->timerDS_ );
00766 ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1);
00767 }
00768 this->om_->stream(Debug) << "S: " << std::endl << S << std::endl;;
00769 TEST_FOR_EXCEPTION(ret != 0,std::logic_error,"Anasazi::IRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret);
00770 TEST_FOR_EXCEPTION(rank != this->blockSize_,RTRRitzFailure,"Anasazi::IRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank);
00771
00772
00773
00774
00775 {
00776 TimeMonitor lcltimer( *this->timerSort_ );
00777 std::vector<int> order(this->blockSize_);
00778
00779 this->sm_->sort(this->theta_, Teuchos::rcpFromRef(order), this->blockSize_);
00780
00781 Utils::permuteVectors(order,S);
00782 }
00783
00784
00785 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
00786
00787
00788
00789 if (this->om_->isVerbosity( Debug ) ) {
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800 MagnitudeType rhonum, rhoden, mxeta;
00801 std::vector<MagnitudeType> eBe(this->blockSize_);
00802 ginnersep(*this->eta_,*this->Beta_,eBe);
00803
00804
00805 rhonum = oldfx - this->fx_;
00806
00807
00808 rhoden = -2.0*ginner(*this->AX_ ,*this->eta_)
00809 -ginner(*this->Aeta_,*this->eta_);
00810 for (int i=0; i<this->blockSize_; ++i) {
00811 rhoden += eBe[i]*oldtheta[i];
00812 }
00813 mxeta = oldfx - rhoden;
00814 this->rho_ = rhonum / rhoden;
00815 this->om_->stream(Debug)
00816 << " >> old f(x) is : " << oldfx << endl
00817 << " >> new f(x) is : " << this->fx_ << endl
00818 << " >> m_x(eta) is : " << mxeta << endl
00819 << " >> rhonum is : " << rhonum << endl
00820 << " >> rhoden is : " << rhoden << endl
00821 << " >> rho is : " << this->rho_ << endl;
00822
00823 for (int j=0; j<this->blockSize_; ++j) {
00824 this->om_->stream(Debug)
00825 << " >> rho[" << j << "] : " << 1.0/(1.0+eBe[j]) << endl;
00826 }
00827 }
00828
00829
00830
00831
00832
00833
00834 {
00835
00836 this->X_ = Teuchos::null;
00837 this->BX_ = Teuchos::null;
00838
00839 std::vector<int> ind(this->blockSize_);
00840 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
00841 Teuchos::RCP<MV> X, BX;
00842 X = MVT::CloneView(*this->V_,ind);
00843 if (this->hasBOp_) {
00844 BX = MVT::CloneView(*this->BV_,ind);
00845 }
00846
00847 {
00848 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
00849 MVT::MvTimesMatAddMv(ONE,*this->delta_,S,ZERO,*X);
00850 MVT::MvTimesMatAddMv(ONE,*this->Adelta_,S,ZERO,*this->AX_);
00851 if (this->hasBOp_) {
00852 MVT::MvTimesMatAddMv(ONE,*this->Bdelta_,S,ZERO,*BX);
00853 }
00854 }
00855
00856 X = Teuchos::null;
00857 BX = Teuchos::null;
00858 this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
00859 this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
00860 }
00861
00862
00863
00864 {
00865 TimeMonitor lcltimer( *this->timerCompRes_ );
00866 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
00867 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00868 MVT::MvScale( *this->R_, theta_comp );
00869 MVT::MvAddMv( ONE, *this->AX_, -ONE, *this->R_, *this->R_ );
00870 }
00871
00872
00873 this->Rnorms_current_ = false;
00874 this->R2norms_current_ = false;
00875
00876
00877
00878
00879 if (this->om_->isVerbosity( Debug ) ) {
00880
00881 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00882 chk.checkX = true;
00883 chk.checkAX = true;
00884 chk.checkBX = true;
00885 chk.checkR = true;
00886 this->om_->print( Debug, this->accuracyCheck(chk, "after local update") );
00887 }
00888 else if (this->om_->isVerbosity( OrthoDetails )) {
00889 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00890 chk.checkX = true;
00891 chk.checkR = true;
00892 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after local update") );
00893 }
00894
00895 }
00896
00897 }
00898
00899
00901
00902 template <class ScalarType, class MV, class OP>
00903 void
00904 IRTR<ScalarType,MV,OP>::currentStatus(std::ostream &os)
00905 {
00906 using std::endl;
00907 os.setf(std::ios::scientific, std::ios::floatfield);
00908 os.precision(6);
00909 os <<endl;
00910 os <<"================================================================================" << endl;
00911 os << endl;
00912 os <<" IRTR Solver Status" << endl;
00913 os << endl;
00914 os <<"The solver is "<<(this->initialized_ ? "initialized." : "not initialized.") << endl;
00915 os <<"The number of iterations performed is " << this->iter_ << endl;
00916 os <<"The current block size is " << this->blockSize_ << endl;
00917 os <<"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
00918 os <<"The number of operations A*x is " << this->counterAOp_ << endl;
00919 os <<"The number of operations B*x is " << this->counterBOp_ << endl;
00920 os <<"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
00921 os <<"The number of operations Prec*x is " << this->counterPrec_ << endl;
00922 os <<"Parameter rho_prime is " << rho_prime_ << endl;
00923 os <<"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
00924 os <<"Number of inner iterations was " << innerIters_ << endl;
00925 os <<"f(x) is " << this->fx_ << endl;
00926
00927 os.setf(std::ios_base::right, std::ios_base::adjustfield);
00928
00929 if (this->initialized_) {
00930 os << endl;
00931 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
00932 os << std::setw(20) << "Eigenvalue"
00933 << std::setw(20) << "Residual(B)"
00934 << std::setw(20) << "Residual(2)"
00935 << endl;
00936 os <<"--------------------------------------------------------------------------------"<<endl;
00937 for (int i=0; i<this->blockSize_; i++) {
00938 os << std::setw(20) << this->theta_[i];
00939 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
00940 else os << std::setw(20) << "not current";
00941 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
00942 else os << std::setw(20) << "not current";
00943 os << endl;
00944 }
00945 }
00946 os <<"================================================================================" << endl;
00947 os << endl;
00948 }
00949
00950
00951 }
00952
00953 #endif // ANASAZI_IRTR_HPP