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

Generated on Thu Dec 17 11:02:21 2009 for Anasazi by  doxygen 1.5.9