|
Anasazi
Version of the Day
|
|
Anasazi
Version of the Day
|
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 // TODO: add randomization 00041 // TODO: add expensive debug checking on Teuchos_Debug 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 // Convenience typedefs 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 // these correspond to above 00111 std::vector<std::string> stopReasons_; 00112 // 00113 // Consts 00114 // 00115 const MagnitudeType ZERO; 00116 const MagnitudeType ONE; 00117 // 00118 // Internal methods 00119 // 00121 void solveTRSubproblem(); 00122 // 00123 // rho_prime 00124 MagnitudeType rho_prime_; 00125 // 00126 // norm of initial gradient: this is used for scaling 00127 MagnitudeType normgradf0_; 00128 // 00129 // tr stopping reason 00130 trRetType innerStop_; 00131 // 00132 // number of inner iterations 00133 int innerIters_, totalInnerIters_; 00134 }; 00135 00136 00137 00138 00140 // Constructor 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 totalInnerIters_(0) 00154 { 00155 // set up array of stop reasons 00156 stopReasons_.push_back("n/a"); 00157 stopReasons_.push_back("maximum iterations"); 00158 stopReasons_.push_back("negative curvature"); 00159 stopReasons_.push_back("exceeded TR"); 00160 stopReasons_.push_back("kappa convergence"); 00161 stopReasons_.push_back("theta convergence"); 00162 00163 rho_prime_ = params.get("Rho Prime",0.5); 00164 TEUCHOS_TEST_FOR_EXCEPTION(rho_prime_ <= 0 || rho_prime_ >= 1,std::invalid_argument, 00165 "Anasazi::SIRTR::constructor: rho_prime must be in (0,1)."); 00166 } 00167 00168 00170 // TR subproblem solver 00171 // 00172 // FINISH: 00173 // define pre- and post-conditions 00174 // 00175 // POST: 00176 // delta_,Adelta_,Hdelta_ undefined 00177 // 00178 template <class ScalarType, class MV, class OP> 00179 void SIRTR<ScalarType,MV,OP>::solveTRSubproblem() { 00180 00181 // return one of: 00182 // MAXIMUM_ITERATIONS 00183 // NEGATIVE_CURVATURE 00184 // EXCEEDED_TR 00185 // KAPPA_CONVERGENCE 00186 // THETA_CONVERGENCE 00187 00188 using Teuchos::RCP; 00189 using Teuchos::tuple; 00190 using Teuchos::null; 00191 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00192 using Teuchos::TimeMonitor; 00193 #endif 00194 using std::endl; 00195 typedef Teuchos::RCP<const MV> PCMV; 00196 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM; 00197 00198 innerStop_ = MAXIMUM_ITERATIONS; 00199 00200 const int n = MVT::GetVecLength(*this->eta_); 00201 const int p = this->blockSize_; 00202 const int d = n*p - (p*p+p)/2; 00203 00204 // We have the following: 00205 // 00206 // X'*B*X = I 00207 // X'*A*X = theta_ 00208 // 00209 // We desire to remain in the trust-region: 00210 // { eta : rho_Y(eta) \geq rho_prime } 00211 // where 00212 // rho_Y(eta) = 1/(1+eta'*B*eta) 00213 // Therefore, the trust-region is 00214 // { eta : eta'*B*eta \leq 1/rho_prime - 1 } 00215 // 00216 const double D2 = ONE/rho_prime_ - ONE; 00217 00218 std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p); 00219 std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p); 00220 MagnitudeType r0_norm; 00221 00222 MVT::MvInit(*this->eta_ ,0.0); 00223 00224 // 00225 // R_ contains direct residuals: 00226 // R_ = A X_ - B X_ diag(theta_) 00227 // 00228 // r0 = grad f(X) = 2 P_BX A X = 2 P_BX (A X - B X diag(theta_) = 2 proj(R_) 00229 // We will do this in place. 00230 // For seeking the rightmost, we want to maximize f 00231 // This is the same as minimizing -f 00232 // Substitute all f with -f here. In particular, 00233 // grad -f(X) = -grad f(X) 00234 // Hess -f(X) = -Hess f(X) 00235 // 00236 { 00237 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00238 TimeMonitor lcltimer( *this->timerOrtho_ ); 00239 #endif 00240 this->orthman_->projectGen( 00241 *this->R_, // operating on R 00242 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I 00243 tuple<PSDM>(null), // don't care about coeffs 00244 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V 00245 if (this->leftMost_) { 00246 MVT::MvScale(*this->R_,2.0); 00247 } 00248 else { 00249 MVT::MvScale(*this->R_,-2.0); 00250 } 00251 } 00252 00253 r0_norm = MAT::squareroot( RTRBase<ScalarType,MV,OP>::ginner(*this->R_) ); 00254 // 00255 // kappa (linear) convergence 00256 // theta (superlinear) convergence 00257 // 00258 MagnitudeType kconv = r0_norm * this->conv_kappa_; 00259 // FINISH: consider inserting some scaling here 00260 // MagnitudeType tconv = r0_norm * MAT::pow(r0_norm/normgradf0_,this->conv_theta_); 00261 MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE); 00262 if (this->om_->isVerbosity(Debug)) { 00263 this->om_->stream(Debug) 00264 << " >> |r0| : " << r0_norm << endl 00265 << " >> kappa conv : " << kconv << endl 00266 << " >> theta conv : " << tconv << endl; 00267 } 00268 00269 // 00270 // For Olsen preconditioning, the preconditioner is 00271 // Z = P_{Prec^-1 BX, BX} Prec^-1 R 00272 // for efficiency, we compute Prec^-1 BX once here for use later 00273 // Otherwise, we don't need PBX 00274 if (this->hasPrec_ && this->olsenPrec_) 00275 { 00276 std::vector<int> ind(this->blockSize_); 00277 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i; 00278 Teuchos::RCP<MV> PBX = MVT::CloneViewNonConst(*this->PBV_,ind); 00279 { 00280 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00281 TimeMonitor prectimer( *this->timerPrec_ ); 00282 #endif 00283 OPT::Apply(*this->Prec_,*this->BX_,*PBX); 00284 this->counterPrec_ += this->blockSize_; 00285 } 00286 PBX = Teuchos::null; 00287 } 00288 00289 // Z = P_{Prec^-1 BV, BV} Prec^-1 R 00290 // Prec^-1 BV in PBV 00291 // or 00292 // Z = P_{BV,BV} Prec^-1 R 00293 if (this->hasPrec_) 00294 { 00295 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00296 TimeMonitor prectimer( *this->timerPrec_ ); 00297 #endif 00298 OPT::Apply(*this->Prec_,*this->R_,*this->Z_); 00299 this->counterPrec_ += this->blockSize_; 00300 // the orthogonalization time counts under Ortho and under Preconditioning 00301 if (this->olsenPrec_) { 00302 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00303 TimeMonitor orthtimer( *this->timerOrtho_ ); 00304 #endif 00305 this->orthman_->projectGen( 00306 *this->Z_, // operating on Z 00307 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false, // P_{PBV,V}, B inner product, and <PBV,V>_B != I 00308 tuple<PSDM>(null), // don't care about coeffs 00309 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*PBV or B*Z, but do have B*V 00310 } 00311 else { 00312 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00313 TimeMonitor orthtimer( *this->timerOrtho_ ); 00314 #endif 00315 this->orthman_->projectGen( 00316 *this->Z_, // operating on Z 00317 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I 00318 tuple<PSDM>(null), // don't care about coeffs 00319 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V 00320 } 00321 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r); 00322 } 00323 else { 00324 // Z = R 00325 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_); 00326 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r); 00327 } 00328 00329 if (this->om_->isVerbosity( Debug )) { 00330 // Check that gradient is B-orthogonal to X 00331 typename RTRBase<ScalarType,MV,OP>::CheckList chk; 00332 chk.checkBR = true; 00333 if (this->hasPrec_) chk.checkZ = true; 00334 this->om_->print( Debug, this->accuracyCheck(chk, "after computing gradient") ); 00335 } 00336 else if (this->om_->isVerbosity( OrthoDetails )) { 00337 // Check that gradient is B-orthogonal to X 00338 typename RTRBase<ScalarType,MV,OP>::CheckList chk; 00339 chk.checkBR = true; 00340 if (this->hasPrec_) chk.checkZ = true; 00341 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after computing gradient") ); 00342 } 00343 00344 // delta = -z 00345 MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_); 00346 00347 if (this->om_->isVerbosity(Debug)) { 00348 // compute the model at eta 00349 // we need Heta, which requires A*eta and B*eta 00350 // we also need A*X 00351 // use Z for storage of these 00352 std::vector<MagnitudeType> eAx(this->blockSize_), 00353 d_eAe(this->blockSize_), 00354 d_eBe(this->blockSize_), 00355 d_mxe(this->blockSize_); 00356 // compute AX and <eta,AX> 00357 { 00358 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00359 TimeMonitor lcltimer( *this->timerAOp_ ); 00360 #endif 00361 OPT::Apply(*this->AOp_,*this->X_,*this->Z_); 00362 this->counterAOp_ += this->blockSize_; 00363 } 00364 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx); 00365 // compute A*eta and <eta,A*eta> 00366 { 00367 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00368 TimeMonitor lcltimer( *this->timerAOp_ ); 00369 #endif 00370 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_); 00371 this->counterAOp_ += this->blockSize_; 00372 } 00373 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe); 00374 // compute B*eta and <eta,B*eta> 00375 if (this->hasBOp_) { 00376 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00377 TimeMonitor lcltimer( *this->timerBOp_ ); 00378 #endif 00379 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_); 00380 this->counterBOp_ += this->blockSize_; 00381 } 00382 else { 00383 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_); 00384 } 00385 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe); 00386 // compute model: 00387 // m_x(eta) = theta + 2*eta'*A*x + eta'*A*eta - eta'*B*eta*theta 00388 if (this->leftMost_) { 00389 for (int j=0; j<this->blockSize_; ++j) { 00390 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j]; 00391 } 00392 } 00393 else { 00394 for (int j=0; j<this->blockSize_; ++j) { 00395 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j]; 00396 } 00397 } 00398 this->om_->stream(Debug) 00399 << " Debugging checks: SIRTR inner iteration " << innerIters_ << endl 00400 << " >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl; 00401 for (int j=0; j<this->blockSize_; ++j) { 00402 this->om_->stream(Debug) 00403 << " >> m_X(eta_" << j << ") : " << d_mxe[j] << endl; 00404 } 00405 } 00406 00409 // the inner/tCG loop 00410 for (innerIters_=1; innerIters_<=d; ++innerIters_) { 00411 00412 // 00413 // [Hdelta,Adelta,Bdelta] = Hess*delta = 2 Proj(A*delta - B*delta*X'*A*X) 00414 // X'*A*X = diag(theta), so that 00415 // (B*delta)*diag(theta) can be done on the cheap 00416 // 00417 { 00418 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00419 TimeMonitor lcltimer( *this->timerAOp_ ); 00420 #endif 00421 OPT::Apply(*this->AOp_,*this->delta_,*this->Z_); 00422 this->counterAOp_ += this->blockSize_; 00423 } 00424 if (this->hasBOp_) { 00425 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00426 TimeMonitor lcltimer( *this->timerBOp_ ); 00427 #endif 00428 OPT::Apply(*this->BOp_,*this->delta_,*this->Hdelta_); 00429 this->counterBOp_ += this->blockSize_; 00430 } 00431 else { 00432 MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Hdelta_); 00433 } 00434 // while we have B*delta, compute <eta,B*delta> and <delta,B*delta> 00435 // these will be needed below 00436 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_ ,*this->Hdelta_,eBd); 00437 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,dBd); 00438 // put 2*A*d - 2*B*d*theta --> Hd 00439 { 00440 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end()); 00441 MVT::MvScale(*this->Hdelta_,theta_comp); 00442 } 00443 if (this->leftMost_) { 00444 MVT::MvAddMv( 2.0,*this->Z_,-2.0,*this->Hdelta_,*this->Hdelta_); 00445 } 00446 else { 00447 MVT::MvAddMv(-2.0,*this->Z_, 2.0,*this->Hdelta_,*this->Hdelta_); 00448 } 00449 // apply projector 00450 { 00451 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00452 TimeMonitor lcltimer( *this->timerOrtho_ ); 00453 #endif 00454 this->orthman_->projectGen( 00455 *this->Hdelta_, // operating on Hdelta 00456 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I 00457 tuple<PSDM>(null), // don't care about coeffs 00458 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V 00459 } 00460 RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,d_Hd); 00461 00462 00463 // compute update step 00464 for (unsigned int j=0; j<alpha.size(); ++j) 00465 { 00466 alpha[j] = z_r[j]/d_Hd[j]; 00467 } 00468 00469 // compute new B-norms 00470 for (unsigned int j=0; j<alpha.size(); ++j) 00471 { 00472 new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j]; 00473 } 00474 00475 if (this->om_->isVerbosity(Debug)) { 00476 for (unsigned int j=0; j<alpha.size(); j++) { 00477 this->om_->stream(Debug) 00478 << " >> z_r[" << j << "] : " << z_r[j] 00479 << " d_Hd[" << j << "] : " << d_Hd[j] << endl 00480 << " >> eBe[" << j << "] : " << eBe[j] 00481 << " neweBe[" << j << "] : " << new_eBe[j] << endl 00482 << " >> eBd[" << j << "] : " << eBd[j] 00483 << " dBd[" << j << "] : " << dBd[j] << endl; 00484 } 00485 } 00486 00487 // check truncation criteria: negative curvature or exceeded trust-region 00488 std::vector<int> trncstep; 00489 trncstep.reserve(p); 00490 // trncstep will contain truncated step, due to 00491 // negative curvature or exceeding implicit trust-region 00492 bool atleastonenegcur = false; 00493 for (unsigned int j=0; j<d_Hd.size(); ++j) { 00494 if (d_Hd[j] <= 0) { 00495 trncstep.push_back(j); 00496 atleastonenegcur = true; 00497 } 00498 else if (new_eBe[j] > D2) { 00499 trncstep.push_back(j); 00500 } 00501 } 00502 00503 if (!trncstep.empty()) 00504 { 00505 // compute step to edge of trust-region, for trncstep vectors 00506 if (this->om_->isVerbosity(Debug)) { 00507 for (unsigned int j=0; j<trncstep.size(); ++j) { 00508 this->om_->stream(Debug) 00509 << " >> alpha[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl; 00510 } 00511 } 00512 for (unsigned int j=0; j<trncstep.size(); ++j) { 00513 int jj = trncstep[j]; 00514 alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj]; 00515 } 00516 if (this->om_->isVerbosity(Debug)) { 00517 for (unsigned int j=0; j<trncstep.size(); ++j) { 00518 this->om_->stream(Debug) 00519 << " >> tau[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl; 00520 } 00521 } 00522 if (atleastonenegcur) { 00523 innerStop_ = NEGATIVE_CURVATURE; 00524 } 00525 else { 00526 innerStop_ = EXCEEDED_TR; 00527 } 00528 } 00529 00530 // compute new eta = eta + alpha*delta 00531 // we need delta*diag(alpha) 00532 // do this in situ in delta_ and friends (we will note this for below) 00533 // then set eta_ = eta_ + delta_ 00534 { 00535 std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end()); 00536 MVT::MvScale(*this->delta_,alpha_comp); 00537 MVT::MvScale(*this->Hdelta_,alpha_comp); 00538 } 00539 MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_); 00540 00541 // store new eBe 00542 eBe = new_eBe; 00543 00544 // 00545 // print some debugging info 00546 if (this->om_->isVerbosity(Debug)) { 00547 // compute the model at eta 00548 // we need Heta, which requires A*eta and B*eta 00549 // we also need A*X 00550 // use Z for storage of these 00551 std::vector<MagnitudeType> eAx(this->blockSize_), 00552 d_eAe(this->blockSize_), 00553 d_eBe(this->blockSize_), 00554 d_mxe(this->blockSize_); 00555 // compute AX and <eta,AX> 00556 { 00557 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00558 TimeMonitor lcltimer( *this->timerAOp_ ); 00559 #endif 00560 OPT::Apply(*this->AOp_,*this->X_,*this->Z_); 00561 this->counterAOp_ += this->blockSize_; 00562 } 00563 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,eAx); 00564 // compute A*eta and <eta,A*eta> 00565 { 00566 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00567 TimeMonitor lcltimer( *this->timerAOp_ ); 00568 #endif 00569 OPT::Apply(*this->AOp_,*this->eta_,*this->Z_); 00570 this->counterAOp_ += this->blockSize_; 00571 } 00572 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eAe); 00573 // compute B*eta and <eta,B*eta> 00574 if (this->hasBOp_) { 00575 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00576 TimeMonitor lcltimer( *this->timerBOp_ ); 00577 #endif 00578 OPT::Apply(*this->BOp_,*this->eta_,*this->Z_); 00579 this->counterBOp_ += this->blockSize_; 00580 } 00581 else { 00582 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*this->Z_); 00583 } 00584 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Z_,d_eBe); 00585 // compute model: 00586 // m_x(eta) = theta + 2*eta'*A*x + eta'*A*eta - eta'*B*eta*theta 00587 if (this->leftMost_) { 00588 for (int j=0; j<this->blockSize_; ++j) { 00589 d_mxe[j] = this->theta_[j] + 2*eAx[j] + d_eAe[j] - d_eBe[j]*this->theta_[j]; 00590 } 00591 } 00592 else { 00593 for (int j=0; j<this->blockSize_; ++j) { 00594 d_mxe[j] = -this->theta_[j] - 2*eAx[j] - d_eAe[j] + d_eBe[j]*this->theta_[j]; 00595 } 00596 } 00597 this->om_->stream(Debug) 00598 << " Debugging checks: SIRTR inner iteration " << innerIters_ << endl 00599 << " >> m_X(eta) : " << std::accumulate(d_mxe.begin(),d_mxe.end(),0.0) << endl; 00600 for (int j=0; j<this->blockSize_; ++j) { 00601 this->om_->stream(Debug) 00602 << " >> m_X(eta_" << j << ") : " << d_mxe[j] << endl; 00603 } 00604 } 00605 00606 // 00607 // if we found negative curvature or exceeded trust-region, then quit 00608 if (!trncstep.empty()) { 00609 break; 00610 } 00611 00612 // update gradient of m 00613 // R = R + Hdelta*diag(alpha) 00614 // however, Hdelta_ already stores Hdelta*diag(alpha) 00615 // so just add them 00616 MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_); 00617 { 00618 // re-tangentialize r 00619 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00620 TimeMonitor lcltimer( *this->timerOrtho_ ); 00621 #endif 00622 this->orthman_->projectGen( 00623 *this->R_, // operating on R 00624 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I 00625 tuple<PSDM>(null), // don't care about coeffs 00626 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V 00627 } 00628 00629 // 00630 // check convergence 00631 MagnitudeType r_norm = MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->R_,*this->R_)); 00632 00633 // 00634 // check local convergece 00635 // 00636 // kappa (linear) convergence 00637 // theta (superlinear) convergence 00638 // 00639 if (this->om_->isVerbosity(Debug)) { 00640 this->om_->stream(Debug) 00641 << " >> |r" << innerIters_ << "| : " << r_norm << endl; 00642 } 00643 if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) { 00644 if (tconv <= kconv) { 00645 innerStop_ = THETA_CONVERGENCE; 00646 } 00647 else { 00648 innerStop_ = KAPPA_CONVERGENCE; 00649 } 00650 break; 00651 } 00652 00653 // Z = P_{Prec^-1 BV, BV} Prec^-1 R 00654 // Prec^-1 BV in PBV 00655 // or 00656 // Z = P_{BV,BV} Prec^-1 R 00657 zold_rold = z_r; 00658 if (this->hasPrec_) 00659 { 00660 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00661 TimeMonitor prectimer( *this->timerPrec_ ); 00662 #endif 00663 OPT::Apply(*this->Prec_,*this->R_,*this->Z_); 00664 this->counterPrec_ += this->blockSize_; 00665 // the orthogonalization time counts under Ortho and under Preconditioning 00666 if (this->olsenPrec_) { 00667 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00668 TimeMonitor orthtimer( *this->timerOrtho_ ); 00669 #endif 00670 this->orthman_->projectGen( 00671 *this->Z_, // operating on Z 00672 tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false, // P_{PBV,V}, B inner product, and <PBV,V>_B != I 00673 tuple<PSDM>(null), // don't care about coeffs 00674 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*PBV or B*Z, but do have B*V 00675 } 00676 else { 00677 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00678 TimeMonitor orthtimer( *this->timerOrtho_ ); 00679 #endif 00680 this->orthman_->projectGen( 00681 *this->Z_, // operating on Z 00682 tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I 00683 tuple<PSDM>(null), // don't care about coeffs 00684 null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V 00685 } 00686 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r); 00687 } 00688 else { 00689 // Z = R 00690 MVT::MvAddMv(ONE,*this->R_,ZERO,*this->R_,*this->Z_); 00691 RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r); 00692 } 00693 00694 // compute new search direction 00695 // below, we need to perform 00696 // delta = -Z + delta*diag(beta) 00697 // however, delta_ currently stores delta*diag(alpha) 00698 // therefore, set 00699 // beta_ to beta/alpha 00700 // so that 00701 // delta_ = delta_*diag(beta_) 00702 // will in fact result in 00703 // delta_ = delta_*diag(beta_) 00704 // = delta*diag(alpha)*diag(beta/alpha) 00705 // = delta*diag(beta) 00706 // i hope this is numerically sound... 00707 for (unsigned int j=0; j<beta.size(); ++j) { 00708 beta[j] = z_r[j]/(zold_rold[j]*alpha[j]); 00709 } 00710 { 00711 std::vector<ScalarType> beta_comp(beta.begin(),beta.end()); 00712 MVT::MvScale(*this->delta_,beta_comp); 00713 } 00714 MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_); 00715 00716 } 00717 // end of the inner iteration loop 00720 if (innerIters_ > d) innerIters_ = d; 00721 00722 this->om_->stream(Debug) 00723 << " >> stop reason is " << stopReasons_[innerStop_] << endl 00724 << endl; 00725 00726 } // end of solveTRSubproblem 00727 00728 00729 #define SIRTR_GET_TEMP_MV(mv,workspace) \ 00730 { \ 00731 TEUCHOS_TEST_FOR_EXCEPTION(workspace.size() == 0,std::logic_error,"SIRTR: Request for workspace could not be honored."); \ 00732 mv = workspace.back(); \ 00733 workspace.pop_back(); \ 00734 } 00735 00736 #define SIRTR_RELEASE_TEMP_MV(mv,workspace) \ 00737 { \ 00738 workspace.push_back(mv); \ 00739 mv = Teuchos::null; \ 00740 } 00741 00742 00744 // Eigensolver iterate() method 00745 template <class ScalarType, class MV, class OP> 00746 void SIRTR<ScalarType,MV,OP>::iterate() { 00747 00748 using Teuchos::RCP; 00749 using Teuchos::null; 00750 using Teuchos::tuple; 00751 using Teuchos::TimeMonitor; 00752 using std::endl; 00753 // typedef Teuchos::RCP<const MV> PCMV; // unused 00754 // typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM; // unused 00755 00756 // 00757 // Allocate/initialize data structures 00758 // 00759 if (this->initialized_ == false) { 00760 this->initialize(); 00761 } 00762 00763 Teuchos::SerialDenseMatrix<int,ScalarType> AA(this->blockSize_,this->blockSize_), 00764 BB(this->blockSize_,this->blockSize_), 00765 S(this->blockSize_,this->blockSize_); 00766 00767 // we will often exploit temporarily unused storage for workspace 00768 // in order to keep it straight and make for clearer code, 00769 // we will put pointers to available multivectors into the following vector 00770 // when we need them, we get them out, using a meaningfully-named pointer 00771 // when we're done, we put them back 00772 std::vector< RCP<MV> > workspace; 00773 // we only have 7 multivectors, so that is more than the maximum number that 00774 // we could use for temp storage 00775 workspace.reserve(7); 00776 00777 // set iteration details to invalid, as they don't have any meaning right now 00778 innerIters_ = -1; 00779 innerStop_ = UNINITIALIZED; 00780 00781 // allocate temporary space 00782 while (this->tester_->checkStatus(this) != Passed) { 00783 00784 // Print information on current status 00785 if (this->om_->isVerbosity(Debug)) { 00786 this->currentStatus( this->om_->stream(Debug) ); 00787 } 00788 else if (this->om_->isVerbosity(IterationDetails)) { 00789 this->currentStatus( this->om_->stream(IterationDetails) ); 00790 } 00791 00792 // increment iteration counter 00793 this->iter_++; 00794 00795 // solve the trust-region subproblem 00796 solveTRSubproblem(); 00797 totalInnerIters_ += innerIters_; 00798 00799 // perform debugging on eta et al. 00800 if (this->om_->isVerbosity( Debug ) ) { 00801 typename RTRBase<ScalarType,MV,OP>::CheckList chk; 00802 // this is the residual of the model, should still be in the tangent plane 00803 chk.checkBR = true; 00804 chk.checkEta = true; 00805 this->om_->print( Debug, this->accuracyCheck(chk, "in iterate() after solveTRSubproblem()") ); 00806 } 00807 00808 00809 // 00810 // multivectors X, BX (if hasB) and eta contain meaningful information that we need below 00811 // the others will be sacrificed to temporary storage 00812 // we are not allowed to reference these anymore, RELEASE_TEMP_MV will clear the pointers 00813 // the RCP in workspace will keep the MV alive, we will get the MVs back 00814 // as we need them using GET_TEMP_MV 00815 // 00816 // this strategy doesn't cost us much, and it keeps us honest 00817 // 00818 TEUCHOS_TEST_FOR_EXCEPTION(workspace.size() != 0,std::logic_error,"SIRTR::iterate(): workspace list should be empty."); 00819 SIRTR_RELEASE_TEMP_MV(this->delta_ ,workspace); // workspace size is 1 00820 SIRTR_RELEASE_TEMP_MV(this->Hdelta_,workspace); // workspace size is 2 00821 SIRTR_RELEASE_TEMP_MV(this->R_ ,workspace); // workspace size is 3 00822 SIRTR_RELEASE_TEMP_MV(this->Z_ ,workspace); // workspace size is 4 00823 00824 00825 // compute the retraction of eta: R_X(eta) = X+eta 00826 // we must accept, but we will work out of temporary so that we can multiply back into X below 00827 RCP<MV> XpEta; 00828 SIRTR_GET_TEMP_MV(XpEta,workspace); // workspace size is 3 00829 { 00830 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00831 TimeMonitor lcltimer( *this->timerLocalUpdate_ ); 00832 #endif 00833 MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*XpEta); 00834 } 00835 00836 // 00837 // perform rayleigh-ritz for XpEta = X+eta 00838 // save an old copy of f(X) for rho analysis below 00839 // 00840 MagnitudeType oldfx = this->fx_; 00841 int rank, ret; 00842 rank = this->blockSize_; 00843 // compute AA = (X+eta)'*A*(X+eta) 00844 // get temporarily storage for A*(X+eta) 00845 RCP<MV> AXpEta; 00846 SIRTR_GET_TEMP_MV(AXpEta,workspace); // workspace size is 2 00847 { 00848 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00849 TimeMonitor lcltimer( *this->timerAOp_ ); 00850 #endif 00851 OPT::Apply(*this->AOp_,*XpEta,*AXpEta); 00852 this->counterAOp_ += this->blockSize_; 00853 } 00854 // project A onto X+eta 00855 { 00856 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00857 TimeMonitor lcltimer( *this->timerLocalProj_ ); 00858 #endif 00859 MVT::MvTransMv(ONE,*XpEta,*AXpEta,AA); 00860 } 00861 // compute BB = (X+eta)'*B*(X+eta) 00862 // get temporary storage for B*(X+eta) 00863 RCP<MV> BXpEta; 00864 if (this->hasBOp_) { 00865 SIRTR_GET_TEMP_MV(BXpEta,workspace); // workspace size is 1 00866 { 00867 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00868 TimeMonitor lcltimer( *this->timerBOp_ ); 00869 #endif 00870 OPT::Apply(*this->BOp_,*XpEta,*BXpEta); 00871 this->counterBOp_ += this->blockSize_; 00872 } 00873 // project B onto X+eta 00874 { 00875 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00876 TimeMonitor lcltimer( *this->timerLocalProj_ ); 00877 #endif 00878 MVT::MvTransMv(ONE,*XpEta,*BXpEta,BB); 00879 } 00880 } 00881 else { 00882 // project I onto X+eta 00883 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00884 TimeMonitor lcltimer( *this->timerLocalProj_ ); 00885 #endif 00886 MVT::MvTransMv(ONE,*XpEta,*XpEta,BB); 00887 } 00888 this->om_->stream(Debug) << "AA: " << std::endl << AA << std::endl;; 00889 this->om_->stream(Debug) << "BB: " << std::endl << BB << std::endl;; 00890 // do the direct solve 00891 // save old theta first 00892 std::vector<MagnitudeType> oldtheta(this->theta_); 00893 { 00894 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00895 TimeMonitor lcltimer( *this->timerDS_ ); 00896 #endif 00897 ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1); 00898 } 00899 this->om_->stream(Debug) << "S: " << std::endl << S << std::endl;; 00900 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); 00901 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,RTRRitzFailure,"Anasazi::SIRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank); 00902 00903 // 00904 // order the projected ritz values and vectors 00905 // this ensures that the ritz vectors produced below are ordered 00906 { 00907 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00908 TimeMonitor lcltimer( *this->timerSort_ ); 00909 #endif 00910 std::vector<int> order(this->blockSize_); 00911 // sort the first blockSize_ values in theta_ 00912 this->sm_->sort(this->theta_, Teuchos::rcpFromRef(order), this->blockSize_); // don't catch exception 00913 // apply the same ordering to the primitive ritz vectors 00914 Utils::permuteVectors(order,S); 00915 } 00916 // 00917 // update f(x) 00918 this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO); 00919 00920 // 00921 // if debugging, do rho analysis before overwriting X,AX,BX 00922 RCP<MV> AX; 00923 SIRTR_GET_TEMP_MV(AX,workspace); // workspace size is 0 00924 if (this->om_->isVerbosity( Debug ) ) { 00925 // 00926 // compute rho 00927 // f(X) - f(X+eta) f(X) - f(X+eta) 00928 // rho = ----------------- = ------------------------- 00929 // m(0) - m(eta) -<2AX,eta> - .5*<Heta,eta> 00930 MagnitudeType rhonum, rhoden, mxeta; 00931 // 00932 // compute rhonum 00933 rhonum = oldfx - this->fx_; 00934 00935 // 00936 // compute rhoden = -<eta,gradfx> - 0.5 <eta,H*eta> 00937 // = -2.0*<eta,AX> - <eta,A*eta> + <eta,B*eta*theta> 00938 // in three steps (3) (1) (2) 00939 // 00940 // first, compute seconder-order decrease in model (steps 1 and 2) 00941 // get temp storage for second order decrease of model 00942 // 00943 // do the first-order decrease last, because we need AX below 00944 { 00945 // compute A*eta and then <eta,A*eta> 00946 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00947 TimeMonitor lcltimer( *this->timerAOp_ ); 00948 #endif 00949 OPT::Apply(*this->AOp_,*this->eta_,*AX); 00950 this->counterAOp_ += this->blockSize_; 00951 } 00952 // compute A part of second order decrease into rhoden 00953 rhoden = -RTRBase<ScalarType,MV,OP>::ginner(*this->eta_,*AX); 00954 if (this->hasBOp_) { 00955 // compute B*eta into AX 00956 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00957 TimeMonitor lcltimer( *this->timerBOp_ ); 00958 #endif 00959 OPT::Apply(*this->BOp_,*this->eta_,*AX); 00960 this->counterBOp_ += this->blockSize_; 00961 } 00962 else { 00963 // put B*eta==eta into AX 00964 MVT::MvAddMv(ONE,*this->eta_,ZERO,*this->eta_,*AX); 00965 } 00966 // we need this below for computing individual rho, get it now 00967 std::vector<MagnitudeType> eBe(this->blockSize_); 00968 RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*AX,eBe); 00969 // scale B*eta by theta 00970 { 00971 std::vector<ScalarType> oldtheta_complex(oldtheta.begin(),oldtheta.end()); 00972 MVT::MvScale( *AX, oldtheta_complex); 00973 } 00974 // accumulate B part of second order decrease into rhoden 00975 rhoden += RTRBase<ScalarType,MV,OP>::ginner(*this->eta_,*AX); 00976 { 00977 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00978 TimeMonitor lcltimer( *this->timerAOp_ ); 00979 #endif 00980 OPT::Apply(*this->AOp_,*this->X_,*AX); 00981 this->counterAOp_ += this->blockSize_; 00982 } 00983 // accumulate first-order decrease of model into rhoden 00984 rhoden += -2.0*RTRBase<ScalarType,MV,OP>::ginner(*AX,*this->eta_); 00985 00986 mxeta = oldfx - rhoden; 00987 this->rho_ = rhonum / rhoden; 00988 this->om_->stream(Debug) 00989 << " >> old f(x) is : " << oldfx << endl 00990 << " >> new f(x) is : " << this->fx_ << endl 00991 << " >> m_x(eta) is : " << mxeta << endl 00992 << " >> rhonum is : " << rhonum << endl 00993 << " >> rhoden is : " << rhoden << endl 00994 << " >> rho is : " << this->rho_ << endl; 00995 // compute individual rho 00996 for (int j=0; j<this->blockSize_; ++j) { 00997 this->om_->stream(Debug) 00998 << " >> rho[" << j << "] : " << 1.0/(1.0+eBe[j]) << endl; 00999 } 01000 } 01001 01002 // compute Ritz vectors back into X,BX,AX 01003 { 01004 // release const views to X, BX 01005 this->X_ = Teuchos::null; 01006 this->BX_ = Teuchos::null; 01007 // get non-const views 01008 std::vector<int> ind(this->blockSize_); 01009 for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i; 01010 Teuchos::RCP<MV> X, BX; 01011 X = MVT::CloneViewNonConst(*this->V_,ind); 01012 if (this->hasBOp_) { 01013 BX = MVT::CloneViewNonConst(*this->BV_,ind); 01014 } 01015 // compute ritz vectors, A,B products into X,AX,BX 01016 { 01017 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 01018 TimeMonitor lcltimer( *this->timerLocalUpdate_ ); 01019 #endif 01020 MVT::MvTimesMatAddMv(ONE,* XpEta,S,ZERO,*X); 01021 MVT::MvTimesMatAddMv(ONE,*AXpEta,S,ZERO,*AX); 01022 if (this->hasBOp_) { 01023 MVT::MvTimesMatAddMv(ONE,*BXpEta,S,ZERO,*BX); 01024 } 01025 } 01026 // clear non-const views, restore const views 01027 X = Teuchos::null; 01028 BX = Teuchos::null; 01029 this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind); 01030 this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind); 01031 } 01032 // 01033 // return XpEta and BXpEta to temp storage 01034 SIRTR_RELEASE_TEMP_MV(XpEta,workspace); // workspace size is 1 01035 SIRTR_RELEASE_TEMP_MV(AXpEta,workspace); // workspace size is 2 01036 if (this->hasBOp_) { 01037 SIRTR_RELEASE_TEMP_MV(BXpEta,workspace); // workspace size is 3 01038 } 01039 01040 // 01041 // solveTRSubproblem destroyed R, we must recompute it 01042 // compute R = AX - BX*theta 01043 // 01044 // get R back from temp storage 01045 SIRTR_GET_TEMP_MV(this->R_,workspace); // workspace size is 2 01046 { 01047 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 01048 TimeMonitor lcltimer( *this->timerCompRes_ ); 01049 #endif 01050 MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ ); 01051 { 01052 std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end()); 01053 MVT::MvScale( *this->R_, theta_comp ); 01054 } 01055 MVT::MvAddMv( ONE, *AX, -ONE, *this->R_, *this->R_ ); 01056 } 01057 // 01058 // R has been updated; mark the norms as out-of-date 01059 this->Rnorms_current_ = false; 01060 this->R2norms_current_ = false; 01061 01062 // 01063 // we are done with AX, release it 01064 SIRTR_RELEASE_TEMP_MV(AX,workspace); // workspace size is 3 01065 // 01066 // get data back for delta, Hdelta and Z 01067 SIRTR_GET_TEMP_MV(this->delta_,workspace); // workspace size is 2 01068 SIRTR_GET_TEMP_MV(this->Hdelta_,workspace); // workspace size is 1 01069 SIRTR_GET_TEMP_MV(this->Z_,workspace); // workspace size is 0 01070 01071 // 01072 // When required, monitor some orthogonalities 01073 if (this->om_->isVerbosity( Debug ) ) { 01074 // Check almost everything here 01075 typename RTRBase<ScalarType,MV,OP>::CheckList chk; 01076 chk.checkX = true; 01077 chk.checkBX = true; 01078 chk.checkR = true; 01079 this->om_->print( Debug, this->accuracyCheck(chk, "after local update") ); 01080 } 01081 else if (this->om_->isVerbosity( OrthoDetails )) { 01082 typename RTRBase<ScalarType,MV,OP>::CheckList chk; 01083 chk.checkX = true; 01084 chk.checkR = true; 01085 this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after local update") ); 01086 } 01087 01088 } // end while (statusTest == false) 01089 01090 } // end of iterate() 01091 01092 01094 // Print the current status of the solver 01095 template <class ScalarType, class MV, class OP> 01096 void 01097 SIRTR<ScalarType,MV,OP>::currentStatus(std::ostream &os) 01098 { 01099 using std::endl; 01100 os.setf(std::ios::scientific, std::ios::floatfield); 01101 os.precision(6); 01102 os <<endl; 01103 os <<"================================================================================" << endl; 01104 os << endl; 01105 os <<" SIRTR Solver Status" << endl; 01106 os << endl; 01107 os <<"The solver is "<<(this->initialized_ ? "initialized." : "not initialized.") << endl; 01108 os <<"The number of iterations performed is " << this->iter_ << endl; 01109 os <<"The current block size is " << this->blockSize_ << endl; 01110 os <<"The number of auxiliary vectors is " << this->numAuxVecs_ << endl; 01111 os <<"The number of operations A*x is " << this->counterAOp_ << endl; 01112 os <<"The number of operations B*x is " << this->counterBOp_ << endl; 01113 os <<"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl; 01114 os <<"The number of operations Prec*x is " << this->counterPrec_ << endl; 01115 os <<"Parameter rho_prime is " << rho_prime_ << endl; 01116 os <<"Inner stopping condition was " << stopReasons_[innerStop_] << endl; 01117 os <<"Number of inner iterations was " << innerIters_ << endl; 01118 os <<"Total number of inner iterations is " << totalInnerIters_ << endl; 01119 os <<"f(x) is " << this->fx_ << endl; 01120 01121 os.setf(std::ios_base::right, std::ios_base::adjustfield); 01122 01123 if (this->initialized_) { 01124 os << endl; 01125 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl; 01126 os << std::setw(20) << "Eigenvalue" 01127 << std::setw(20) << "Residual(B)" 01128 << std::setw(20) << "Residual(2)" 01129 << endl; 01130 os <<"--------------------------------------------------------------------------------"<<endl; 01131 for (int i=0; i<this->blockSize_; i++) { 01132 os << std::setw(20) << this->theta_[i]; 01133 if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i]; 01134 else os << std::setw(20) << "not current"; 01135 if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i]; 01136 else os << std::setw(20) << "not current"; 01137 os << endl; 01138 } 01139 } 01140 os <<"================================================================================" << endl; 01141 os << endl; 01142 } 01143 01144 01145 } // end Anasazi namespace 01146 01147 #endif // ANASAZI_SIRTR_HPP
1.7.6.1
