00001
00002
00003
00004 #ifndef ANASAZI_SIRTR_HPP
00005 #define ANASAZI_SIRTR_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 SIRTR : public RTRBase<ScalarType,MV,OP> {
00047 public:
00048
00050
00051
00063 SIRTR( 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 ~SIRTR() {};
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 SIRTR<ScalarType,MV,OP>::SIRTR(
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,"SIRTR",true),
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::SIRTR::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 SIRTR<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
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 {
00234 TimeMonitor lcltimer( *this->timerOrtho_ );
00235 this->orthman_->projectGen(
00236 *this->R_,
00237 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00238 tuple<PSDM>(null),
00239 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00240 if (this->leftMost_) {
00241 MVT::MvScale(*this->R_,2.0);
00242 }
00243 else {
00244 MVT::MvScale(*this->R_,-2.0);
00245 }
00246 }
00247
00248 r0_norm = MAT::squareroot( ginner(*this->R_) );
00249
00250
00251
00252
00253 MagnitudeType kconv = r0_norm * this->conv_kappa_;
00254
00255
00256 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
00257 if (this->om_->isVerbosity(Debug)) {
00258 this->om_->stream(Debug)
00259 << " >> |r0| : " << r0_norm << endl
00260 << " >> kappa conv : " << kconv << endl
00261 << " >> theta conv : " << tconv << endl;
00262 }
00263
00264
00265
00266
00267
00268 if (this->hasPrec_)
00269 {
00270 std::vector<int> ind(this->blockSize_);
00271 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
00272 Teuchos::RCP<MV> PBX = MVT::CloneView(*this->PBV_,ind);
00273 {
00274 TimeMonitor prectimer( *this->timerPrec_ );
00275 OPT::Apply(*this->Prec_,*this->BX_,*PBX);
00276 this->counterPrec_ += this->blockSize_;
00277 }
00278 PBX = Teuchos::null;
00279 }
00280
00281
00282
00283 if (this->hasPrec_)
00284 {
00285 TimeMonitor prectimer( *this->timerPrec_ );
00286 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
00287 this->counterPrec_ += this->blockSize_;
00288
00289 {
00290 TimeMonitor orthtimer( *this->timerOrtho_ );
00291 this->orthman_->projectGen(
00292 *this->Z_,
00293 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false,
00294 tuple<PSDM>(null),
00295 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00296 }
00297 ginnersep(*this->R_,*this->Z_,z_r);
00298 }
00299 else {
00300
00301 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
00302 ginnersep(*this->R_,z_r);
00303 }
00304
00305 if (this->om_->isVerbosity( Debug )) {
00306
00307 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00308 chk.checkBR = true;
00309 if (this->hasPrec_) chk.checkZ = true;
00310 this->om_->print( Debug, this->accuracyCheck(chk, "after computing gradient") );
00311 }
00312 else if (this->om_->isVerbosity( OrthoDetails )) {
00313
00314 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00315 chk.checkBR = true;
00316 if (this->hasPrec_) chk.checkZ = true;
00317 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after computing gradient") );
00318 }
00319
00320
00321 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
00322
00323 if (this->om_->isVerbosity(Debug)) {
00324
00325
00326
00327
00328 std::vector<MagnitudeType> eAx(this->blockSize_),
00329 d_eAe(this->blockSize_),
00330 d_eBe(this->blockSize_),
00331 d_mxe(this->blockSize_);
00332
00333 {
00334 TimeMonitor lcltimer( *this->timerAOp_ );
00335 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
00336 this->counterAOp_ += this->blockSize_;
00337 }
00338 ginnersep(*this->eta_,*this->Z_,eAx);
00339
00340 {
00341 TimeMonitor lcltimer( *this->timerAOp_ );
00342 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
00343 this->counterAOp_ += this->blockSize_;
00344 }
00345 ginnersep(*this->eta_,*this->Z_,d_eAe);
00346
00347 if (this->hasBOp_) {
00348 TimeMonitor lcltimer( *this->timerBOp_ );
00349 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
00350 this->counterBOp_ += this->blockSize_;
00351 }
00352 else {
00353 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
00354 }
00355 ginnersep(*this->eta_,*this->Z_,d_eBe);
00356
00357
00358 if (this->leftMost_) {
00359 for (int j=0; j<this->blockSize_; ++j) {
00360 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
00361 }
00362 }
00363 else {
00364 for (int j=0; j<this->blockSize_; ++j) {
00365 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
00366 }
00367 }
00368 this->om_->stream(Debug)
00369 << " Debugging checks: SIRTR inner iteration " << innerIters_ << endl
00370 << " >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
00371 for (int j=0; j<this->blockSize_; ++j) {
00372 this->om_->stream(Debug)
00373 << " >> m_X(eta_" << j << ") : " << d_mxe[j] << endl;
00374 }
00375 }
00376
00379
00380 for (innerIters_=1; innerIters_<=d; ++innerIters_) {
00381
00382
00383
00384
00385
00386
00387 {
00388 TimeMonitor lcltimer( *this->timerAOp_ );
00389 OPT::Apply(*this->AOp_,*this->delta_,*this->Z_);
00390 this->counterAOp_ += this->blockSize_;
00391 }
00392 if (this->hasBOp_) {
00393 TimeMonitor lcltimer( *this->timerBOp_ );
00394 OPT::Apply(*this->BOp_,*this->delta_,*this->Hdelta_);
00395 this->counterBOp_ += this->blockSize_;
00396 }
00397 else {
00398 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Hdelta_);
00399 }
00400
00401
00402 ginnersep(*this->eta_ ,*this->Hdelta_,eBd);
00403 ginnersep(*this->delta_,*this->Hdelta_,dBd);
00404
00405 {
00406 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00407 MVT::MvScale(*this->Hdelta_,theta_comp);
00408 }
00409 if (this->leftMost_) {
00410 MVT::MvAddMv( 2.0,*this->Z_,-2.0,*this->Hdelta_,*this->Hdelta_);
00411 }
00412 else {
00413 MVT::MvAddMv(-2.0,*this->Z_, 2.0,*this->Hdelta_,*this->Hdelta_);
00414 }
00415
00416 {
00417 TimeMonitor lcltimer( *this->timerOrtho_ );
00418 this->orthman_->projectGen(
00419 *this->Hdelta_,
00420 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00421 tuple<PSDM>(null),
00422 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00423 }
00424 ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
00425
00426
00427
00428 for (unsigned int j=0; j<alpha.size(); ++j)
00429 {
00430 alpha[j] = z_r[j]/d_Hd[j];
00431 }
00432
00433
00434 for (unsigned int j=0; j<alpha.size(); ++j)
00435 {
00436 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
00437 }
00438
00439 if (this->om_->isVerbosity(Debug)) {
00440 for (unsigned int j=0; j<alpha.size(); j++) {
00441 this->om_->stream(Debug)
00442 << " >> z_r[" << j << "] : " << z_r[j]
00443 << " d_Hd[" << j << "] : " << d_Hd[j] << endl
00444 << " >> eBe[" << j << "] : " << eBe[j]
00445 << " neweBe[" << j << "] : " << new_eBe[j] << endl
00446 << " >> eBd[" << j << "] : " << eBd[j]
00447 << " dBd[" << j << "] : " << dBd[j] << endl;
00448 }
00449 }
00450
00451
00452 std::vector<int> trncstep;
00453 trncstep.reserve(p);
00454
00455
00456 bool atleastonenegcur = false;
00457 for (unsigned int j=0; j<d_Hd.size(); ++j) {
00458 if (d_Hd[j] <= 0) {
00459 trncstep.push_back(j);
00460 atleastonenegcur = true;
00461 }
00462 else if (new_eBe[j] > D2) {
00463 trncstep.push_back(j);
00464 }
00465 }
00466
00467 if (!trncstep.empty())
00468 {
00469
00470 if (this->om_->isVerbosity(Debug)) {
00471 for (unsigned int j=0; j<trncstep.size(); ++j) {
00472 this->om_->stream(Debug)
00473 << " >> alpha[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
00474 }
00475 }
00476 for (unsigned int j=0; j<trncstep.size(); ++j) {
00477 int jj = trncstep[j];
00478 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
00479 }
00480 if (this->om_->isVerbosity(Debug)) {
00481 for (unsigned int j=0; j<trncstep.size(); ++j) {
00482 this->om_->stream(Debug)
00483 << " >> tau[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
00484 }
00485 }
00486 if (atleastonenegcur) {
00487 innerStop_ = NEGATIVE_CURVATURE;
00488 }
00489 else {
00490 innerStop_ = EXCEEDED_TR;
00491 }
00492 }
00493
00494
00495
00496
00497
00498 {
00499 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
00500 MVT::MvScale(*this->delta_,alpha_comp);
00501 MVT::MvScale(*this->Hdelta_,alpha_comp);
00502 }
00503 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
00504
00505
00506 eBe = new_eBe;
00507
00508
00509
00510 if (this->om_->isVerbosity(Debug)) {
00511
00512
00513
00514
00515 std::vector<MagnitudeType> eAx(this->blockSize_),
00516 d_eAe(this->blockSize_),
00517 d_eBe(this->blockSize_),
00518 d_mxe(this->blockSize_);
00519
00520 {
00521 TimeMonitor lcltimer( *this->timerAOp_ );
00522 OPT::Apply(*this->AOp_,*this->X_,*this->Z_);
00523 this->counterAOp_ += this->blockSize_;
00524 }
00525 ginnersep(*this->eta_,*this->Z_,eAx);
00526
00527 {
00528 TimeMonitor lcltimer( *this->timerAOp_ );
00529 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_);
00530 this->counterAOp_ += this->blockSize_;
00531 }
00532 ginnersep(*this->eta_,*this->Z_,d_eAe);
00533
00534 if (this->hasBOp_) {
00535 TimeMonitor lcltimer( *this->timerBOp_ );
00536 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_);
00537 this->counterBOp_ += this->blockSize_;
00538 }
00539 else {
00540 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_);
00541 }
00542 ginnersep(*this->eta_,*this->Z_,d_eBe);
00543
00544
00545 if (this->leftMost_) {
00546 for (int j=0; j<this->blockSize_; ++j) {
00547 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j];
00548 }
00549 }
00550 else {
00551 for (int j=0; j<this->blockSize_; ++j) {
00552 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j];
00553 }
00554 }
00555 this->om_->stream(Debug)
00556 << " Debugging checks: SIRTR inner iteration " << innerIters_ << endl
00557 << " >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl;
00558 for (int j=0; j<this->blockSize_; ++j) {
00559 this->om_->stream(Debug)
00560 << " >> m_X(eta_" << j << ") : " << d_mxe[j] << endl;
00561 }
00562 }
00563
00564
00565
00566 if (!trncstep.empty()) {
00567 break;
00568 }
00569
00570
00571
00572
00573
00574 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
00575 {
00576
00577 TimeMonitor lcltimer( *this->timerOrtho_ );
00578 this->orthman_->projectGen(
00579 *this->R_,
00580 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,
00581 tuple<PSDM>(null),
00582 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00583 }
00584
00585
00586
00587 MagnitudeType r_norm = MAT::squareroot(ginner(*this->R_,*this->R_));
00588
00589
00590
00591
00592
00593
00594
00595 if (this->om_->isVerbosity(Debug)) {
00596 this->om_->stream(Debug)
00597 << " >> |r" << innerIters_ << "| : " << r_norm << endl;
00598 }
00599 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
00600 if (tconv <= kconv) {
00601 innerStop_ = THETA_CONVERGENCE;
00602 }
00603 else {
00604 innerStop_ = KAPPA_CONVERGENCE;
00605 }
00606 break;
00607 }
00608
00609
00610
00611 zold_rold = z_r;
00612 if (this->hasPrec_)
00613 {
00614 TimeMonitor prectimer( *this->timerPrec_ );
00615 OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
00616 this->counterPrec_ += this->blockSize_;
00617
00618 {
00619 TimeMonitor orthtimer( *this->timerOrtho_ );
00620 this->orthman_->projectGen(
00621 *this->Z_,
00622 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false,
00623 tuple<PSDM>(null),
00624 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));
00625 }
00626 ginnersep(*this->R_,*this->Z_,z_r);
00627 }
00628 else {
00629
00630 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_);
00631 ginnersep(*this->R_,z_r);
00632 }
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 for (unsigned int j=0; j<beta.size(); ++j) {
00648 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
00649 }
00650 {
00651 std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
00652 MVT::MvScale(*this->delta_,beta_comp);
00653 }
00654 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
00655
00656 }
00657
00660 if (innerIters_ > d) innerIters_ = d;
00661
00662 this->om_->stream(Debug)
00663 << " >> stop reason is " << stopReasons_[innerStop_] << endl
00664 << endl;
00665
00666 }
00667
00668
00669 #define SIRTR_GET_TEMP_MV(mv,workspace) \
00670 { \
00671 TEST_FOR_EXCEPTION(workspace.size() == 0,std::logic_error,"SIRTR: Request for workspace could not be honored."); \
00672 mv = workspace.back(); \
00673 workspace.pop_back(); \
00674 }
00675
00676 #define SIRTR_RELEASE_TEMP_MV(mv,workspace) \
00677 { \
00678 workspace.push_back(mv); \
00679 mv = Teuchos::null; \
00680 }
00681
00682
00684
00685 template <class ScalarType, class MV, class OP>
00686 void SIRTR<ScalarType,MV,OP>::iterate() {
00687
00688 using Teuchos::RCP;
00689 using Teuchos::null;
00690 using Teuchos::tuple;
00691 using Teuchos::TimeMonitor;
00692 using std::endl;
00693 typedef Teuchos::RCP<const MV> PCMV;
00694 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
00695
00696
00697
00698
00699 if (this->initialized_ == false) {
00700 this->initialize();
00701 }
00702
00703 Teuchos::SerialDenseMatrix<int,ScalarType> AA(this->blockSize_,this->blockSize_),
00704 BB(this->blockSize_,this->blockSize_),
00705 S(this->blockSize_,this->blockSize_);
00706
00707
00708
00709
00710
00711
00712 std::vector< RCP<MV> > workspace;
00713
00714
00715 workspace.reserve(7);
00716
00717
00718 innerIters_ = -1;
00719 innerStop_ = UNINITIALIZED;
00720
00721
00722 while (this->tester_->checkStatus(this) != Passed) {
00723
00724
00725 if (this->om_->isVerbosity(Debug)) {
00726 this->currentStatus( this->om_->stream(Debug) );
00727 }
00728 else if (this->om_->isVerbosity(IterationDetails)) {
00729 this->currentStatus( this->om_->stream(IterationDetails) );
00730 }
00731
00732
00733 this->iter_++;
00734
00735
00736 solveTRSubproblem();
00737
00738
00739 if (this->om_->isVerbosity( Debug ) ) {
00740 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00741
00742 chk.checkBR = true;
00743 chk.checkEta = true;
00744 this->om_->print( Debug, this->accuracyCheck(chk, "in iterate() after solveTRSubproblem()") );
00745 }
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 TEST_FOR_EXCEPTION(workspace.size() != 0,std::logic_error,"SIRTR::iterate(): workspace list should be empty.");
00758 SIRTR_RELEASE_TEMP_MV(this->delta_ ,workspace);
00759 SIRTR_RELEASE_TEMP_MV(this->Hdelta_,workspace);
00760 SIRTR_RELEASE_TEMP_MV(this->R_ ,workspace);
00761 SIRTR_RELEASE_TEMP_MV(this->Z_ ,workspace);
00762
00763
00764
00765
00766 RCP<MV> XpEta;
00767 SIRTR_GET_TEMP_MV(XpEta,workspace);
00768 {
00769 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
00770 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*XpEta);
00771 }
00772
00773
00774
00775
00776
00777 MagnitudeType oldfx = this->fx_;
00778 int rank, ret;
00779 rank = this->blockSize_;
00780
00781
00782 RCP<MV> AXpEta;
00783 SIRTR_GET_TEMP_MV(AXpEta,workspace);
00784 {
00785 TimeMonitor lcltimer( *this->timerAOp_ );
00786 OPT::Apply(*this->AOp_,*XpEta,*AXpEta);
00787 this->counterAOp_ += this->blockSize_;
00788 }
00789
00790 {
00791 TimeMonitor lcltimer( *this->timerLocalProj_ );
00792 MVT::MvTransMv(ONE,*XpEta,*AXpEta,AA);
00793 }
00794
00795
00796 RCP<MV> BXpEta;
00797 if (this->hasBOp_) {
00798 SIRTR_GET_TEMP_MV(BXpEta,workspace);
00799 {
00800 TimeMonitor lcltimer( *this->timerBOp_ );
00801 OPT::Apply(*this->BOp_,*XpEta,*BXpEta);
00802 this->counterBOp_ += this->blockSize_;
00803 }
00804
00805 {
00806 TimeMonitor lcltimer( *this->timerLocalProj_ );
00807 MVT::MvTransMv(ONE,*XpEta,*BXpEta,BB);
00808 }
00809 }
00810 else {
00811
00812 TimeMonitor lcltimer( *this->timerLocalProj_ );
00813 MVT::MvTransMv(ONE,*XpEta,*XpEta,BB);
00814 }
00815 this->om_->stream(Debug) << "AA: " << std::endl << AA << std::endl;;
00816 this->om_->stream(Debug) << "BB: " << std::endl << BB << std::endl;;
00817
00818
00819 std::vector<MagnitudeType> oldtheta(this->theta_);
00820 {
00821 TimeMonitor lcltimer( *this->timerDS_ );
00822 ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1);
00823 }
00824 this->om_->stream(Debug) << "S: " << std::endl << S << std::endl;;
00825 TEST_FOR_EXCEPTION(ret != 0,std::logic_error,"Anasazi::SIRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret);
00826 TEST_FOR_EXCEPTION(rank != this->blockSize_,RTRRitzFailure,"Anasazi::SIRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank);
00827
00828
00829
00830
00831 {
00832 TimeMonitor lcltimer( *this->timerSort_ );
00833 std::vector<int> order(this->blockSize_);
00834
00835 this->sm_->sort(this->theta_, Teuchos::rcpFromRef(order), this->blockSize_);
00836
00837 Utils::permuteVectors(order,S);
00838 }
00839
00840
00841 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
00842
00843
00844
00845 RCP<MV> AX;
00846 SIRTR_GET_TEMP_MV(AX,workspace);
00847 if (this->om_->isVerbosity( Debug ) ) {
00848
00849
00850
00851
00852
00853 MagnitudeType rhonum, rhoden, mxeta;
00854
00855
00856 rhonum = oldfx - this->fx_;
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867 {
00868
00869 TimeMonitor lcltimer( *this->timerAOp_ );
00870 OPT::Apply(*this->AOp_,*this->eta_,*AX);
00871 this->counterAOp_ += this->blockSize_;
00872 }
00873
00874 rhoden = -ginner(*this->eta_,*AX);
00875 if (this->hasBOp_) {
00876
00877 TimeMonitor lcltimer( *this->timerBOp_ );
00878 OPT::Apply(*this->BOp_,*this->eta_,*AX);
00879 this->counterBOp_ += this->blockSize_;
00880 }
00881 else {
00882
00883 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*AX);
00884 }
00885
00886 std::vector<MagnitudeType> eBe(this->blockSize_);
00887 ginnersep(*this->eta_,*AX,eBe);
00888
00889 {
00890 std::vector<ScalarType> oldtheta_complex(oldtheta.begin(),oldtheta.end());
00891 MVT::MvScale( *AX, oldtheta_complex);
00892 }
00893
00894 rhoden += ginner(*this->eta_,*AX);
00895 {
00896 TimeMonitor lcltimer( *this->timerAOp_ );
00897 OPT::Apply(*this->AOp_,*this->X_,*AX);
00898 this->counterAOp_ += this->blockSize_;
00899 }
00900
00901 rhoden += -2.0*ginner(*AX,*this->eta_);
00902
00903 mxeta = oldfx - rhoden;
00904 this->rho_ = rhonum / rhoden;
00905 this->om_->stream(Debug)
00906 << " >> old f(x) is : " << oldfx << endl
00907 << " >> new f(x) is : " << this->fx_ << endl
00908 << " >> m_x(eta) is : " << mxeta << endl
00909 << " >> rhonum is : " << rhonum << endl
00910 << " >> rhoden is : " << rhoden << endl
00911 << " >> rho is : " << this->rho_ << endl;
00912
00913 for (int j=0; j<this->blockSize_; ++j) {
00914 this->om_->stream(Debug)
00915 << " >> rho[" << j << "] : " << 1.0/(1.0+eBe[j]) << endl;
00916 }
00917 }
00918
00919
00920 {
00921
00922 this->X_ = Teuchos::null;
00923 this->BX_ = Teuchos::null;
00924
00925 std::vector<int> ind(this->blockSize_);
00926 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
00927 Teuchos::RCP<MV> X, BX;
00928 X = MVT::CloneView(*this->V_,ind);
00929 if (this->hasBOp_) {
00930 BX = MVT::CloneView(*this->BV_,ind);
00931 }
00932
00933 {
00934 TimeMonitor lcltimer( *this->timerLocalUpdate_ );
00935 MVT::MvTimesMatAddMv(ONE,* XpEta,S,ZERO,*X);
00936 MVT::MvTimesMatAddMv(ONE,*AXpEta,S,ZERO,*AX);
00937 if (this->hasBOp_) {
00938 MVT::MvTimesMatAddMv(ONE,*BXpEta,S,ZERO,*BX);
00939 }
00940 }
00941
00942 X = Teuchos::null;
00943 BX = Teuchos::null;
00944 this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
00945 this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
00946 }
00947
00948
00949 SIRTR_RELEASE_TEMP_MV(XpEta,workspace);
00950 SIRTR_RELEASE_TEMP_MV(AXpEta,workspace);
00951 if (this->hasBOp_) {
00952 SIRTR_RELEASE_TEMP_MV(BXpEta,workspace);
00953 }
00954
00955
00956
00957
00958
00959
00960 SIRTR_GET_TEMP_MV(this->R_,workspace);
00961 {
00962 TimeMonitor lcltimer( *this->timerCompRes_ );
00963 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
00964 {
00965 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
00966 MVT::MvScale( *this->R_, theta_comp );
00967 }
00968 MVT::MvAddMv( ONE, *AX, -ONE, *this->R_, *this->R_ );
00969 }
00970
00971
00972 this->Rnorms_current_ = false;
00973 this->R2norms_current_ = false;
00974
00975
00976
00977 SIRTR_RELEASE_TEMP_MV(AX,workspace);
00978
00979
00980 SIRTR_GET_TEMP_MV(this->delta_,workspace);
00981 SIRTR_GET_TEMP_MV(this->Hdelta_,workspace);
00982 SIRTR_GET_TEMP_MV(this->Z_,workspace);
00983
00984
00985
00986 if (this->om_->isVerbosity( Debug ) ) {
00987
00988 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00989 chk.checkX = true;
00990 chk.checkBX = true;
00991 chk.checkR = true;
00992 this->om_->print( Debug, this->accuracyCheck(chk, "after local update") );
00993 }
00994 else if (this->om_->isVerbosity( OrthoDetails )) {
00995 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
00996 chk.checkX = true;
00997 chk.checkR = true;
00998 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after local update") );
00999 }
01000
01001 }
01002
01003 }
01004
01005
01007
01008 template <class ScalarType, class MV, class OP>
01009 void
01010 SIRTR<ScalarType,MV,OP>::currentStatus(std::ostream &os)
01011 {
01012 using std::endl;
01013 os.setf(std::ios::scientific, std::ios::floatfield);
01014 os.precision(6);
01015 os <<endl;
01016 os <<"================================================================================" << endl;
01017 os << endl;
01018 os <<" SIRTR Solver Status" << endl;
01019 os << endl;
01020 os <<"The solver is "<<(this->initialized_ ? "initialized." : "not initialized.") << endl;
01021 os <<"The number of iterations performed is " << this->iter_ << endl;
01022 os <<"The current block size is " << this->blockSize_ << endl;
01023 os <<"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
01024 os <<"The number of operations A*x is " << this->counterAOp_ << endl;
01025 os <<"The number of operations B*x is " << this->counterBOp_ << endl;
01026 os <<"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
01027 os <<"The number of operations Prec*x is " << this->counterPrec_ << endl;
01028 os <<"Parameter rho_prime is " << rho_prime_ << endl;
01029 os <<"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
01030 os <<"Number of inner iterations was " << innerIters_ << endl;
01031 os <<"f(x) is " << this->fx_ << endl;
01032
01033 os.setf(std::ios_base::right, std::ios_base::adjustfield);
01034
01035 if (this->initialized_) {
01036 os << endl;
01037 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
01038 os << std::setw(20) << "Eigenvalue"
01039 << std::setw(20) << "Residual(B)"
01040 << std::setw(20) << "Residual(2)"
01041 << endl;
01042 os <<"--------------------------------------------------------------------------------"<<endl;
01043 for (int i=0; i<this->blockSize_; i++) {
01044 os << std::setw(20) << this->theta_[i];
01045 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
01046 else os << std::setw(20) << "not current";
01047 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
01048 else os << std::setw(20) << "not current";
01049 os << endl;
01050 }
01051 }
01052 os <<"================================================================================" << endl;
01053 os << endl;
01054 }
01055
01056
01057 }
01058
01059 #endif // ANASAZI_SIRTR_HPP