Anasazi Version of the Day
AnasaziIRTR.hpp
Go to the documentation of this file.
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                                 &params 
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                                 &params
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;
00732     typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends