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