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