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

Generated on Wed Oct 21 14:28:10 2009 for Anasazi by  doxygen 1.5.9