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