Anasazi  Version of the Day
AnasaziIRTR.hpp
Go to the documentation of this file.
1 
3 
4 #ifndef ANASAZI_IRTR_HPP
5 #define ANASAZI_IRTR_HPP
6 
7 #include "AnasaziTypes.hpp"
8 #include "AnasaziRTRBase.hpp"
9 
10 #include "AnasaziEigensolver.hpp"
13 #include "Teuchos_ScalarTraits.hpp"
14 
15 #include "Teuchos_LAPACK.hpp"
16 #include "Teuchos_BLAS.hpp"
19 #include "Teuchos_TimeMonitor.hpp"
20 
40 // TODO: add randomization
41 // TODO: add expensive debug checking on Teuchos_Debug
42 
43 namespace Anasazi {
44 
45  template <class ScalarType, class MV, class OP>
46  class IRTR : public RTRBase<ScalarType,MV,OP> {
47  public:
48 
50 
51 
65  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
69  );
70 
72  virtual ~IRTR() {};
73 
75 
77 
78 
80  void iterate();
81 
83 
85 
86 
88  void currentStatus(std::ostream &os);
89 
91 
92  private:
93  //
94  // Convenience typedefs
95  //
96  typedef SolverUtils<ScalarType,MV,OP> Utils;
100  typedef typename SCT::magnitudeType MagnitudeType;
102  enum trRetType {
103  UNINITIALIZED = 0,
104  MAXIMUM_ITERATIONS,
105  NEGATIVE_CURVATURE,
106  EXCEEDED_TR,
107  KAPPA_CONVERGENCE,
108  THETA_CONVERGENCE
109  };
110  // these correspond to above
111  std::vector<std::string> stopReasons_;
112  //
113  // Consts
114  //
115  const MagnitudeType ZERO;
116  const MagnitudeType ONE;
117  //
118  // Internal methods
119  //
121  void solveTRSubproblem();
122  //
123  // rho_prime
124  MagnitudeType rho_prime_;
125  //
126  // norm of initial gradient: this is used for scaling
127  MagnitudeType normgradf0_;
128  //
129  // tr stopping reason
130  trRetType innerStop_;
131  //
132  // number of inner iterations
133  int innerIters_, totalInnerIters_;
134  //
135  // use 2D subspace acceleration of X+Eta to generate new iterate?
136  bool useSA_;
137  };
138 
139 
140 
141 
143  // Constructor
144  template <class ScalarType, class MV, class OP>
148  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
151  Teuchos::ParameterList &params
152  ) :
153  RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,"IRTR",false),
154  ZERO(MAT::zero()),
155  ONE(MAT::one()),
156  totalInnerIters_(0)
157  {
158  // set up array of stop reasons
159  stopReasons_.push_back("n/a");
160  stopReasons_.push_back("maximum iterations");
161  stopReasons_.push_back("negative curvature");
162  stopReasons_.push_back("exceeded TR");
163  stopReasons_.push_back("kappa convergence");
164  stopReasons_.push_back("theta convergence");
165 
166  rho_prime_ = params.get("Rho Prime",0.5);
167  TEUCHOS_TEST_FOR_EXCEPTION(rho_prime_ <= 0 || rho_prime_ >= 1,std::invalid_argument,
168  "Anasazi::IRTR::constructor: rho_prime must be in (0,1).");
169 
170  useSA_ = params.get<bool>("Use SA",false);
171  }
172 
173 
175  // TR subproblem solver
176  //
177  // FINISH:
178  // define pre- and post-conditions
179  //
180  // POST:
181  // delta_,Adelta_,Bdelta_,Hdelta_ undefined
182  //
183  template <class ScalarType, class MV, class OP>
185 
186  // return one of:
187  // MAXIMUM_ITERATIONS
188  // NEGATIVE_CURVATURE
189  // EXCEEDED_TR
190  // KAPPA_CONVERGENCE
191  // THETA_CONVERGENCE
192 
193  using Teuchos::RCP;
194  using Teuchos::tuple;
195  using Teuchos::null;
196 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
197  using Teuchos::TimeMonitor;
198 #endif
199  using std::endl;
200  typedef Teuchos::RCP<const MV> PCMV;
202 
203  innerStop_ = MAXIMUM_ITERATIONS;
204 
205  const int n = MVT::GetGlobalLength(*this->eta_);
206  const int p = this->blockSize_;
207  const int d = n*p - (p*p+p)/2;
208 
209  // We have the following:
210  //
211  // X'*B*X = I
212  // X'*A*X = theta_
213  //
214  // We desire to remain in the trust-region:
215  // { eta : rho_Y(eta) \geq rho_prime }
216  // where
217  // rho_Y(eta) = 1/(1+eta'*B*eta)
218  // Therefore, the trust-region is
219  // { eta : eta'*B*eta \leq 1/rho_prime - 1 }
220  //
221  const double D2 = ONE/rho_prime_ - ONE;
222 
223  std::vector<MagnitudeType> d_Hd(p), alpha(p), beta(p), z_r(p), zold_rold(p);
224  std::vector<MagnitudeType> eBe(p), eBd(p), dBd(p), new_eBe(p);
225  MagnitudeType r0_norm;
226 
227  MVT::MvInit(*this->eta_ ,0.0);
228  MVT::MvInit(*this->Aeta_,0.0);
229  if (this->hasBOp_) {
230  MVT::MvInit(*this->Beta_,0.0);
231  }
232 
233  //
234  // R_ contains direct residuals:
235  // R_ = A X_ - B X_ diag(theta_)
236  //
237  // r0 = grad f(X) = 2 P_BX A X = 2 P_BX (A X - B X diag(theta_) = 2 proj(R_)
238  // We will do this in place.
239  // For seeking the rightmost, we want to maximize f
240  // This is the same as minimizing -f
241  // Substitute all f with -f here. In particular,
242  // grad -f(X) = -grad f(X)
243  // Hess -f(X) = -Hess f(X)
244  //
245  {
246 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
247  TimeMonitor lcltimer( *this->timerOrtho_ );
248 #endif
249  this->orthman_->projectGen(
250  *this->R_, // operating on R
251  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
252  tuple<PSDM>(null), // don't care about coeffs
253  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
254  if (this->leftMost_) {
255  MVT::MvScale(*this->R_,2.0);
256  }
257  else {
258  MVT::MvScale(*this->R_,-2.0);
259  }
260  }
261 
262  r0_norm = MAT::squareroot( RTRBase<ScalarType,MV,OP>::ginner(*this->R_) );
263  //
264  // kappa (linear) convergence
265  // theta (superlinear) convergence
266  //
267  MagnitudeType kconv = r0_norm * this->conv_kappa_;
268  // FINISH: consider inserting some scaling here
269  // MagnitudeType tconv = r0_norm * MAT::pow(r0_norm/normgradf0_,this->conv_theta_);
270  MagnitudeType tconv = MAT::pow(r0_norm,this->conv_theta_+ONE);
271  if (this->om_->isVerbosity(Debug)) {
272  this->om_->stream(Debug)
273  << " >> |r0| : " << r0_norm << endl
274  << " >> kappa conv : " << kconv << endl
275  << " >> theta conv : " << tconv << endl;
276  }
277 
278  //
279  // For Olsen preconditioning, the preconditioner is
280  // Z = P_{Prec^-1 BX, BX} Prec^-1 R
281  // for efficiency, we compute Prec^-1 BX once here for use later
282  // Otherwise, we don't need PBX
283  if (this->hasPrec_ && this->olsenPrec_)
284  {
285  std::vector<int> ind(this->blockSize_);
286  for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
287  Teuchos::RCP<MV> PBX = MVT::CloneViewNonConst(*this->PBV_,ind);
288  {
289 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
290  TimeMonitor prectimer( *this->timerPrec_ );
291 #endif
292  OPT::Apply(*this->Prec_,*this->BX_,*PBX);
293  this->counterPrec_ += this->blockSize_;
294  }
295  PBX = Teuchos::null;
296  }
297 
298  // Z = P_{Prec^-1 BV, BV} Prec^-1 R
299  // Prec^-1 BV in PBV
300  // or
301  // Z = P_{BV,BV} Prec^-1 R
302  if (this->hasPrec_)
303  {
304 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
305  TimeMonitor prectimer( *this->timerPrec_ );
306 #endif
307  OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
308  this->counterPrec_ += this->blockSize_;
309  // the orthogonalization time counts under Ortho and under Preconditioning
310  if (this->olsenPrec_) {
311 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
312  TimeMonitor orthtimer( *this->timerOrtho_ );
313 #endif
314  this->orthman_->projectGen(
315  *this->Z_, // operating on Z
316  tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false, // P_{PBV,V}, B inner product, and <PBV,V>_B != I
317  tuple<PSDM>(null), // don't care about coeffs
318  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*PBV or B*Z, but do have B*V
319  }
320  else {
321 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
322  TimeMonitor orthtimer( *this->timerOrtho_ );
323 #endif
324  this->orthman_->projectGen(
325  *this->Z_, // operating on Z
326  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
327  tuple<PSDM>(null), // don't care about coeffs
328  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
329  }
330  RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
331  }
332  else {
333  RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
334  }
335 
336  if (this->om_->isVerbosity( Debug )) {
337  // Check that gradient is B-orthogonal to X
338  typename RTRBase<ScalarType,MV,OP>::CheckList chk;
339  chk.checkBR = true;
340  if (this->hasPrec_) chk.checkZ = true;
341  this->om_->print( Debug, this->accuracyCheck(chk, "after computing gradient") );
342  }
343  else if (this->om_->isVerbosity( OrthoDetails )) {
344  // Check that gradient is B-orthogonal to X
345  typename RTRBase<ScalarType,MV,OP>::CheckList chk;
346  chk.checkBR = true;
347  if (this->hasPrec_) chk.checkZ = true;
348  this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after computing gradient") );
349  }
350 
351  // delta = -z
352  MVT::MvAddMv(-ONE,*this->Z_,ZERO,*this->Z_,*this->delta_);
353 
354  if (this->om_->isVerbosity(Debug)) {
355  RCP<MV> Heta = MVT::Clone(*this->eta_,this->blockSize_);
356  // put 2*A*e - 2*B*e*theta --> He
357  MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
358  std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
359  MVT::MvScale(*Heta,theta_comp);
360  if (this->leftMost_) {
361  MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta); // Heta = 2*Aeta + (-2)*Beta*theta
362  }
363  else {
364  MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta); // Heta = (-2)*Aeta + 2*Beta*theta
365  }
366  // apply projector
367  {
368 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
369  TimeMonitor lcltimer( *this->timerOrtho_ );
370 #endif
371  this->orthman_->projectGen(
372  *Heta, // operating on Heta
373  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
374  tuple<PSDM>(null), // don't care about coeffs
375  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
376  }
377 
378  std::vector<MagnitudeType> eg(this->blockSize_),
379  eHe(this->blockSize_);
380  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->AX_,eg);
381  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*Heta,eHe);
382  if (this->leftMost_) {
383  for (int j=0; j<this->blockSize_; ++j) {
384  eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
385  }
386  }
387  else {
388  for (int j=0; j<this->blockSize_; ++j) {
389  eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
390  }
391  }
392  this->om_->stream(Debug)
393  << " Debugging checks: IRTR inner iteration " << innerIters_ << endl
394  << " >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
395  for (int j=0; j<this->blockSize_; ++j) {
396  this->om_->stream(Debug)
397  << " >> m_X(eta_" << j << ") : " << eg[j] << endl;
398  }
399  }
400 
403  // the inner/tCG loop
404  for (innerIters_=1; innerIters_<=d; ++innerIters_) {
405 
406  //
407  // [Hdelta,Adelta,Bdelta] = Hess*delta = 2 Proj(A*delta - B*delta*X'*A*X)
408  // X'*A*X = diag(theta), so that
409  // (B*delta)*diag(theta) can be done on the cheap
410  //
411  {
412 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
413  TimeMonitor lcltimer( *this->timerAOp_ );
414 #endif
415  OPT::Apply(*this->AOp_,*this->delta_,*this->Adelta_);
416  this->counterAOp_ += this->blockSize_;
417  }
418  if (this->hasBOp_) {
419 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
420  TimeMonitor lcltimer( *this->timerBOp_ );
421 #endif
422  OPT::Apply(*this->BOp_,*this->delta_,*this->Bdelta_);
423  this->counterBOp_ += this->blockSize_;
424  }
425  else {
426  MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->Bdelta_);
427  }
428  // compute <eta,B*delta> and <delta,B*delta>
429  // these will be needed below
430  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_ ,*this->Bdelta_,eBd);
431  RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Bdelta_,dBd);
432  // put 2*A*d - 2*B*d*theta --> Hd
433  MVT::MvAddMv(ONE,*this->Bdelta_,ZERO,*this->Bdelta_,*this->Hdelta_);
434  {
435  std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
436  MVT::MvScale(*this->Hdelta_,theta_comp);
437  }
438  if (this->leftMost_) {
439  MVT::MvAddMv( 2.0,*this->Adelta_,-2.0,*this->Hdelta_,*this->Hdelta_);
440  }
441  else {
442  MVT::MvAddMv(-2.0,*this->Adelta_, 2.0,*this->Hdelta_,*this->Hdelta_);
443  }
444  // apply projector
445  {
446 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
447  TimeMonitor lcltimer( *this->timerOrtho_ );
448 #endif
449  this->orthman_->projectGen(
450  *this->Hdelta_, // operating on Hdelta
451  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
452  tuple<PSDM>(null), // don't care about coeffs
453  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
454  }
455  RTRBase<ScalarType,MV,OP>::ginnersep(*this->delta_,*this->Hdelta_,d_Hd);
456 
457 
458  // compute update step
459  for (unsigned int j=0; j<alpha.size(); ++j)
460  {
461  alpha[j] = z_r[j]/d_Hd[j];
462  }
463 
464  // compute new B-norms
465  for (unsigned int j=0; j<alpha.size(); ++j)
466  {
467  new_eBe[j] = eBe[j] + 2*alpha[j]*eBd[j] + alpha[j]*alpha[j]*dBd[j];
468  }
469 
470  if (this->om_->isVerbosity(Debug)) {
471  for (unsigned int j=0; j<alpha.size(); j++) {
472  this->om_->stream(Debug)
473  << " >> z_r[" << j << "] : " << z_r[j]
474  << " d_Hd[" << j << "] : " << d_Hd[j] << endl
475  << " >> eBe[" << j << "] : " << eBe[j]
476  << " neweBe[" << j << "] : " << new_eBe[j] << endl
477  << " >> eBd[" << j << "] : " << eBd[j]
478  << " dBd[" << j << "] : " << dBd[j] << endl;
479  }
480  }
481 
482  // check truncation criteria: negative curvature or exceeded trust-region
483  std::vector<int> trncstep;
484  trncstep.reserve(p);
485  // trncstep will contain truncated step, due to
486  // negative curvature or exceeding implicit trust-region
487  bool atleastonenegcur = false;
488  for (unsigned int j=0; j<d_Hd.size(); ++j) {
489  if (d_Hd[j] <= 0) {
490  trncstep.push_back(j);
491  atleastonenegcur = true;
492  }
493  else if (new_eBe[j] > D2) {
494  trncstep.push_back(j);
495  }
496  }
497 
498  if (!trncstep.empty())
499  {
500  // compute step to edge of trust-region, for trncstep vectors
501  if (this->om_->isVerbosity(Debug)) {
502  for (unsigned int j=0; j<trncstep.size(); ++j) {
503  this->om_->stream(Debug)
504  << " >> alpha[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
505  }
506  }
507  for (unsigned int j=0; j<trncstep.size(); ++j) {
508  int jj = trncstep[j];
509  alpha[jj] = ( -eBd[jj] + MAT::squareroot(eBd[jj]*eBd[jj] + dBd[jj]*(D2-eBe[jj]) ) ) / dBd[jj];
510  }
511  if (this->om_->isVerbosity(Debug)) {
512  for (unsigned int j=0; j<trncstep.size(); ++j) {
513  this->om_->stream(Debug)
514  << " >> tau[" << trncstep[j] << "] : " << alpha[trncstep[j]] << endl;
515  }
516  }
517  if (atleastonenegcur) {
518  innerStop_ = NEGATIVE_CURVATURE;
519  }
520  else {
521  innerStop_ = EXCEEDED_TR;
522  }
523  }
524 
525  // compute new eta = eta + alpha*delta
526  // we need delta*diag(alpha)
527  // do this in situ in delta_ and friends (we will note this for below)
528  // then set eta_ = eta_ + delta_
529  {
530  std::vector<ScalarType> alpha_comp(alpha.begin(),alpha.end());
531  MVT::MvScale(*this->delta_,alpha_comp);
532  MVT::MvScale(*this->Adelta_,alpha_comp);
533  if (this->hasBOp_) {
534  MVT::MvScale(*this->Bdelta_,alpha_comp);
535  }
536  MVT::MvScale(*this->Hdelta_,alpha_comp);
537  }
538  MVT::MvAddMv(ONE,*this->delta_ ,ONE,*this->eta_ ,*this->eta_);
539  MVT::MvAddMv(ONE,*this->Adelta_,ONE,*this->Aeta_,*this->Aeta_);
540  if (this->hasBOp_) {
541  MVT::MvAddMv(ONE,*this->Bdelta_,ONE,*this->Beta_,*this->Beta_);
542  }
543 
544  // store new eBe
545  eBe = new_eBe;
546 
547  //
548  // print some debugging info
549  if (this->om_->isVerbosity(Debug)) {
550  RCP<MV> Heta = MVT::Clone(*this->eta_,this->blockSize_);
551  // put 2*A*e - 2*B*e*theta --> He
552  MVT::MvAddMv(ONE,*this->Beta_,ZERO,*this->Beta_,*Heta);
553  {
554  std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
555  MVT::MvScale(*Heta,theta_comp);
556  }
557  if (this->leftMost_) {
558  MVT::MvAddMv( 2.0,*this->Aeta_,-2.0,*Heta,*Heta);
559  }
560  else {
561  MVT::MvAddMv(-2.0,*this->Aeta_, 2.0,*Heta,*Heta);
562  }
563  // apply projector
564  {
565 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
566  TimeMonitor lcltimer( *this->timerOrtho_ );
567 #endif
568  this->orthman_->projectGen(
569  *Heta, // operating on Heta
570  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
571  tuple<PSDM>(null), // don't care about coeffs
572  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
573  }
574 
575  std::vector<MagnitudeType> eg(this->blockSize_),
576  eHe(this->blockSize_);
577  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->AX_,eg);
578  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*Heta,eHe);
579  if (this->leftMost_) {
580  for (int j=0; j<this->blockSize_; ++j) {
581  eg[j] = this->theta_[j] + 2*eg[j] + .5*eHe[j];
582  }
583  }
584  else {
585  for (int j=0; j<this->blockSize_; ++j) {
586  eg[j] = -this->theta_[j] - 2*eg[j] + .5*eHe[j];
587  }
588  }
589  this->om_->stream(Debug)
590  << " Debugging checks: IRTR inner iteration " << innerIters_ << endl
591  << " >> m_X(eta) : " << std::accumulate(eg.begin(),eg.end(),0.0) << endl;
592  for (int j=0; j<this->blockSize_; ++j) {
593  this->om_->stream(Debug)
594  << " >> m_X(eta_" << j << ") : " << eg[j] << endl;
595  }
596  }
597 
598  //
599  // if we found negative curvature or exceeded trust-region, then quit
600  if (!trncstep.empty()) {
601  break;
602  }
603 
604  // update gradient of m
605  // R = R + Hdelta*diag(alpha)
606  // however, Hdelta_ already stores Hdelta*diag(alpha)
607  // so just add them
608  MVT::MvAddMv(ONE,*this->Hdelta_,ONE,*this->R_,*this->R_);
609  {
610  // re-tangentialize r
611 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
612  TimeMonitor lcltimer( *this->timerOrtho_ );
613 #endif
614  this->orthman_->projectGen(
615  *this->R_, // operating on R
616  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
617  tuple<PSDM>(null), // don't care about coeffs
618  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
619  }
620 
621  //
622  // check convergence
623  MagnitudeType r_norm = MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->R_,*this->R_));
624 
625  //
626  // check local convergece
627  //
628  // kappa (linear) convergence
629  // theta (superlinear) convergence
630  //
631  if (this->om_->isVerbosity(Debug)) {
632  this->om_->stream(Debug)
633  << " >> |r" << innerIters_ << "| : " << r_norm << endl;
634  }
635  if ( r_norm <= ANASAZI_MIN(tconv,kconv) ) {
636  if (tconv <= kconv) {
637  innerStop_ = THETA_CONVERGENCE;
638  }
639  else {
640  innerStop_ = KAPPA_CONVERGENCE;
641  }
642  break;
643  }
644 
645  // Z = P_{Prec^-1 BV, BV} Prec^-1 R
646  // Prec^-1 BV in PBV
647  // or
648  // Z = P_{BV,BV} Prec^-1 R
649  zold_rold = z_r;
650  if (this->hasPrec_)
651  {
652 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
653  TimeMonitor prectimer( *this->timerPrec_ );
654 #endif
655  OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
656  this->counterPrec_ += this->blockSize_;
657  // the orthogonalization time counts under Ortho and under Preconditioning
658  if (this->olsenPrec_) {
659 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
660  TimeMonitor orthtimer( *this->timerOrtho_ );
661 #endif
662  this->orthman_->projectGen(
663  *this->Z_, // operating on Z
664  tuple<PCMV>(this->PBV_),tuple<PCMV>(this->V_),false, // P_{PBV,V}, B inner product, and <PBV,V>_B != I
665  tuple<PSDM>(null), // don't care about coeffs
666  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*PBV or B*Z, but do have B*V
667  }
668  else {
669 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
670  TimeMonitor orthtimer( *this->timerOrtho_ );
671 #endif
672  this->orthman_->projectGen(
673  *this->Z_, // operating on Z
674  tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false, // P_{BV,V}, and <BV,V>_B != I
675  tuple<PSDM>(null), // don't care about coeffs
676  null,tuple<PCMV>(null), tuple<PCMV>(this->BV_)); // don't have B*BV, but do have B*V
677  }
678  RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,*this->Z_,z_r);
679  }
680  else {
681  RTRBase<ScalarType,MV,OP>::ginnersep(*this->R_,z_r);
682  }
683 
684  // compute new search direction
685  // below, we need to perform
686  // delta = -Z + delta*diag(beta)
687  // however, delta_ currently stores delta*diag(alpha)
688  // therefore, set
689  // beta_ to beta/alpha
690  // so that
691  // delta_ = delta_*diag(beta_)
692  // will in fact result in
693  // delta_ = delta_*diag(beta_)
694  // = delta*diag(alpha)*diag(beta/alpha)
695  // = delta*diag(beta)
696  // i hope this is numerically sound...
697  for (unsigned int j=0; j<beta.size(); ++j) {
698  beta[j] = z_r[j]/(zold_rold[j]*alpha[j]);
699  }
700  {
701  std::vector<ScalarType> beta_comp(beta.begin(),beta.end());
702  MVT::MvScale(*this->delta_,beta_comp);
703  }
704  MVT::MvAddMv(-ONE,*this->Z_,ONE,*this->delta_,*this->delta_);
705 
706  }
707  // end of the inner iteration loop
710  if (innerIters_ > d) innerIters_ = d;
711 
712  this->om_->stream(Debug)
713  << " >> stop reason is " << stopReasons_[innerStop_] << endl
714  << endl;
715 
716  } // end of solveTRSubproblem
717 
718 
720  // Eigensolver iterate() method
721  template <class ScalarType, class MV, class OP>
723 
724  using Teuchos::RCP;
725  using Teuchos::null;
726  using Teuchos::tuple;
727 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
728  using Teuchos::TimeMonitor;
729 #endif
730  using std::endl;
731  //typedef Teuchos::RCP<const MV> PCMV; // unused
732  //typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM; // unused
733 
734  //
735  // Allocate/initialize data structures
736  //
737  if (this->initialized_ == false) {
738  this->initialize();
739  }
740 
742  if (useSA_ == true) {
743  AA.reshape(2*this->blockSize_,2*this->blockSize_);
744  BB.reshape(2*this->blockSize_,2*this->blockSize_);
745  S.reshape(2*this->blockSize_,2*this->blockSize_);
746  }
747  else {
748  AA.reshape(this->blockSize_,this->blockSize_);
749  BB.reshape(this->blockSize_,this->blockSize_);
750  S.reshape(this->blockSize_,this->blockSize_);
751  }
752 
753  // set iteration details to invalid, as they don't have any meaning right now
754  innerIters_ = -1;
755  innerStop_ = UNINITIALIZED;
756 
757  // allocate temporary space
758  while (this->tester_->checkStatus(this) != Passed) {
759 
760  // Print information on current status
761  if (this->om_->isVerbosity(Debug)) {
762  this->currentStatus( this->om_->stream(Debug) );
763  }
764  else if (this->om_->isVerbosity(IterationDetails)) {
765  this->currentStatus( this->om_->stream(IterationDetails) );
766  }
767 
768  // increment iteration counter
769  this->iter_++;
770 
771  // solve the trust-region subproblem
772  solveTRSubproblem();
773  totalInnerIters_ += innerIters_;
774 
775  // perform debugging on eta et al.
776  if (this->om_->isVerbosity( Debug ) ) {
778  // this is the residual of the model, should still be in the tangent plane
779  chk.checkBR = true;
780  chk.checkEta = true;
781  chk.checkAEta = true;
782  chk.checkBEta = true;
783  this->om_->print( Debug, this->accuracyCheck(chk, "in iterate() after solveTRSubproblem()") );
784  this->om_->stream(Debug)
785  << " >> norm(Eta) : " << MAT::squareroot(RTRBase<ScalarType,MV,OP>::ginner(*this->eta_)) << endl
786  << endl;
787  }
788 
789  if (useSA_ == false) {
790  // compute the retraction of eta: R_X(eta) = X+eta
791  // we must accept, but we will work out of delta so that we can multiply back into X below
792  {
793 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
794  TimeMonitor lcltimer( *this->timerLocalUpdate_ );
795 #endif
796  MVT::MvAddMv(ONE,*this->X_,ONE,*this->eta_,*this->delta_);
797  MVT::MvAddMv(ONE,*this->AX_,ONE,*this->Aeta_,*this->Adelta_);
798  if (this->hasBOp_) {
799  MVT::MvAddMv(ONE,*this->BX_,ONE,*this->Beta_,*this->Bdelta_);
800  }
801  }
802  // perform some debugging on X+eta
803  if (this->om_->isVerbosity( Debug ) ) {
804  // X^T B X = I
805  // X^T B Eta = 0
806  // (X+Eta)^T B (X+Eta) = I + Eta^T B Eta
807  Teuchos::SerialDenseMatrix<int,ScalarType> XE(this->blockSize_,this->blockSize_),
808  E(this->blockSize_,this->blockSize_);
809  MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,XE);
810  MVT::MvTransMv(ONE,*this->eta_,*this->Beta_,E);
811  this->om_->stream(Debug)
812  << " >> Error in AX+AEta == A(X+Eta) : " << Utils::errorEquality(*this->delta_,*this->Adelta_,this->AOp_) << endl
813  << " >> Error in BX+BEta == B(X+Eta) : " << Utils::errorEquality(*this->delta_,*this->Bdelta_,this->BOp_) << endl
814  << " >> norm( (X+Eta)^T B (X+Eta) ) : " << XE.normFrobenius() << endl
815  << " >> norm( Eta^T B Eta ) : " << E.normFrobenius() << endl
816  << endl;
817  }
818  }
819 
820  //
821  // perform rayleigh-ritz for X+eta or [X,eta] according to useSA_
822  // save an old copy of f(X) for rho analysis below
823  //
824  MagnitudeType oldfx = this->fx_;
825  std::vector<MagnitudeType> oldtheta(this->theta_), newtheta(AA.numRows());
826  int rank, ret;
827  rank = AA.numRows();
828  if (useSA_ == true) {
829 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
830  TimeMonitor lcltimer( *this->timerLocalProj_ );
831 #endif
832  Teuchos::SerialDenseMatrix<int,ScalarType> AA11(Teuchos::View,AA,this->blockSize_,this->blockSize_,0,0),
833  AA12(Teuchos::View,AA,this->blockSize_,this->blockSize_,0,this->blockSize_),
834  AA22(Teuchos::View,AA,this->blockSize_,this->blockSize_,this->blockSize_,this->blockSize_);
835  Teuchos::SerialDenseMatrix<int,ScalarType> BB11(Teuchos::View,BB,this->blockSize_,this->blockSize_,0,0),
836  BB12(Teuchos::View,BB,this->blockSize_,this->blockSize_,0,this->blockSize_),
837  BB22(Teuchos::View,BB,this->blockSize_,this->blockSize_,this->blockSize_,this->blockSize_);
838  MVT::MvTransMv(ONE,*this->X_ ,*this->AX_ ,AA11);
839  MVT::MvTransMv(ONE,*this->X_ ,*this->Aeta_,AA12);
840  MVT::MvTransMv(ONE,*this->eta_,*this->Aeta_,AA22);
841  MVT::MvTransMv(ONE,*this->X_ ,*this->BX_ ,BB11);
842  MVT::MvTransMv(ONE,*this->X_ ,*this->Beta_,BB12);
843  MVT::MvTransMv(ONE,*this->eta_,*this->Beta_,BB22);
844  }
845  else {
846 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
847  TimeMonitor lcltimer( *this->timerLocalProj_ );
848 #endif
849  MVT::MvTransMv(ONE,*this->delta_,*this->Adelta_,AA);
850  MVT::MvTransMv(ONE,*this->delta_,*this->Bdelta_,BB);
851  }
852  this->om_->stream(Debug) << "AA: " << std::endl << AA << std::endl;;
853  this->om_->stream(Debug) << "BB: " << std::endl << BB << std::endl;;
854  {
855 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
856  TimeMonitor lcltimer( *this->timerDS_ );
857 #endif
858  ret = Utils::directSolver(AA.numRows(),AA,Teuchos::rcpFromRef(BB),S,newtheta,rank,1);
859  }
860  this->om_->stream(Debug) << "S: " << std::endl << S << std::endl;;
861  TEUCHOS_TEST_FOR_EXCEPTION(ret != 0,std::logic_error,"Anasazi::IRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret);
862  TEUCHOS_TEST_FOR_EXCEPTION(rank != AA.numRows(),RTRRitzFailure,"Anasazi::IRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank);
863 
864  //
865  // order the projected ritz values and vectors
866  // this ensures that the ritz vectors produced below are ordered
867  {
868 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
869  TimeMonitor lcltimer( *this->timerSort_ );
870 #endif
871  std::vector<int> order(newtheta.size());
872  // sort the values in newtheta
873  this->sm_->sort(newtheta, Teuchos::rcpFromRef(order), -1); // don't catch exception
874  // apply the same ordering to the primitive ritz vectors
875  Utils::permuteVectors(order,S);
876  }
877  //
878  // save the first blockSize values into this->theta_
879  std::copy(newtheta.begin(), newtheta.begin()+this->blockSize_, this->theta_.begin());
880  //
881  // update f(x)
882  this->fx_ = std::accumulate(this->theta_.begin(),this->theta_.end(),ZERO);
883 
884  //
885  // if debugging, do rho analysis before overwriting X,AX,BX
886  if (this->om_->isVerbosity( Debug ) ) {
887  //
888  // compute rho
889  // f(X) - f(newX) f(X) - f(newX)
890  // rho = ---------------- = ---------------------------
891  // m(0) - m(eta) -<2AX,eta> - .5*<Heta,eta>
892  //
893  // f(X) - f(newX)
894  // = ---------------------------------------
895  // -<2AX,eta> - <eta,Aeta> + <eta,Beta XAX>
896  //
897  MagnitudeType rhonum, rhoden, mxeta;
898  std::vector<MagnitudeType> eBe(this->blockSize_);
899  RTRBase<ScalarType,MV,OP>::ginnersep(*this->eta_,*this->Beta_,eBe);
900  //
901  // compute rhonum
902  rhonum = oldfx - this->fx_;
903  //
904  // compute rhoden
905  rhoden = -2.0*RTRBase<ScalarType,MV,OP>::ginner(*this->AX_ ,*this->eta_)
906  -RTRBase<ScalarType,MV,OP>::ginner(*this->Aeta_,*this->eta_);
907  for (int i=0; i<this->blockSize_; ++i) {
908  rhoden += eBe[i]*oldtheta[i];
909  }
910  mxeta = oldfx - rhoden;
911  this->rho_ = rhonum / rhoden;
912  this->om_->stream(Debug)
913  << " >> old f(x) is : " << oldfx << endl
914  << " >> new f(x) is : " << this->fx_ << endl
915  << " >> m_x(eta) is : " << mxeta << endl
916  << " >> rhonum is : " << rhonum << endl
917  << " >> rhoden is : " << rhoden << endl
918  << " >> rho is : " << this->rho_ << endl;
919  // compute individual rho
920  for (int j=0; j<this->blockSize_; ++j) {
921  this->om_->stream(Debug)
922  << " >> rho[" << j << "] : " << 1.0/(1.0+eBe[j]) << endl;
923  }
924  }
925 
926  // form new X as Ritz vectors, using the primitive Ritz vectors in S and
927  // either [X eta] or X+eta
928  // we will clear the const views of X,BX into V,BV and
929  // work from non-const temporary views
930  {
931  // release const views to X, BX
932  // get non-const views
933  this->X_ = Teuchos::null;
934  this->BX_ = Teuchos::null;
935  std::vector<int> ind(this->blockSize_);
936  for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
937  Teuchos::RCP<MV> X, BX;
938  X = MVT::CloneViewNonConst(*this->V_,ind);
939  if (this->hasBOp_) {
940  BX = MVT::CloneViewNonConst(*this->BV_,ind);
941  }
942  if (useSA_ == false) {
943  // multiply delta=(X+eta),Adelta=...,Bdelta=...
944  // by primitive Ritz vectors back into X,AX,BX
945  // compute ritz vectors, A,B products into X,AX,BX
946 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
947  TimeMonitor lcltimer( *this->timerLocalUpdate_ );
948 #endif
949  MVT::MvTimesMatAddMv(ONE,*this->delta_,S,ZERO,*X);
950  MVT::MvTimesMatAddMv(ONE,*this->Adelta_,S,ZERO,*this->AX_);
951  if (this->hasBOp_) {
952  MVT::MvTimesMatAddMv(ONE,*this->Bdelta_,S,ZERO,*BX);
953  }
954  }
955  else {
956  // compute ritz vectors, A,B products into X,AX,BX
957  // currently, X in X and eta in eta
958  // compute each result into delta, then copy to appropriate place
959  // decompose S into [Sx;Se]
960  Teuchos::SerialDenseMatrix<int,ScalarType> Sx(Teuchos::View,S,this->blockSize_,this->blockSize_,0,0),
961  Se(Teuchos::View,S,this->blockSize_,this->blockSize_,this->blockSize_,0);
962 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
963  TimeMonitor lcltimer( *this->timerLocalUpdate_ );
964 #endif
965  // X = [X eta] S = X*Sx + eta*Se
966  MVT::MvTimesMatAddMv(ONE,* X ,Sx,ZERO,*this->delta_);
967  MVT::MvTimesMatAddMv(ONE,*this->eta_,Se,ONE ,*this->delta_);
968  MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*X);
969  // AX = [AX Aeta] S = AX*Sx + Aeta*Se
970  MVT::MvTimesMatAddMv(ONE,*this->AX_ ,Sx,ZERO,*this->delta_);
971  MVT::MvTimesMatAddMv(ONE,*this->Aeta_,Se,ONE ,*this->delta_);
972  MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*this->AX_);
973  if (this->hasBOp_) {
974  // BX = [BX Beta] S = BX*Sx + Beta*Se
975  MVT::MvTimesMatAddMv(ONE,* BX ,Sx,ZERO,*this->delta_);
976  MVT::MvTimesMatAddMv(ONE,*this->Beta_,Se,ONE ,*this->delta_);
977  MVT::MvAddMv(ONE,*this->delta_,ZERO,*this->delta_,*BX);
978  }
979  }
980  // clear non-const views, restore const views
981  X = Teuchos::null;
982  BX = Teuchos::null;
983  this->X_ = MVT::CloneView(static_cast<const MV&>(*this->V_ ),ind);
984  this->BX_ = MVT::CloneView(static_cast<const MV&>(*this->BV_),ind);
985  }
986 
987  //
988  // update residual, R = AX - BX*theta
989  {
990 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
991  TimeMonitor lcltimer( *this->timerCompRes_ );
992 #endif
993  MVT::MvAddMv( ONE, *this->BX_, ZERO, *this->BX_, *this->R_ );
994  std::vector<ScalarType> theta_comp(this->theta_.begin(),this->theta_.end());
995  MVT::MvScale( *this->R_, theta_comp );
996  MVT::MvAddMv( ONE, *this->AX_, -ONE, *this->R_, *this->R_ );
997  }
998  //
999  // R has been updated; mark the norms as out-of-date
1000  this->Rnorms_current_ = false;
1001  this->R2norms_current_ = false;
1002 
1003 
1004  //
1005  // When required, monitor some orthogonalities
1006  if (this->om_->isVerbosity( Debug ) ) {
1007  // Check almost everything here
1009  chk.checkX = true;
1010  chk.checkAX = true;
1011  chk.checkBX = true;
1012  chk.checkR = true;
1013  this->om_->print( Debug, this->accuracyCheck(chk, "after local update") );
1014  }
1015  else if (this->om_->isVerbosity( OrthoDetails )) {
1017  chk.checkX = true;
1018  chk.checkR = true;
1019  this->om_->print( OrthoDetails, this->accuracyCheck(chk, "after local update") );
1020  }
1021 
1022  } // end while (statusTest == false)
1023 
1024  } // end of iterate()
1025 
1026 
1028  // Print the current status of the solver
1029  template <class ScalarType, class MV, class OP>
1030  void
1032  {
1033  using std::endl;
1034  os.setf(std::ios::scientific, std::ios::floatfield);
1035  os.precision(6);
1036  os <<endl;
1037  os <<"================================================================================" << endl;
1038  os << endl;
1039  os <<" IRTR Solver Status" << endl;
1040  os << endl;
1041  os <<"The solver is "<<(this->initialized_ ? "initialized." : "not initialized.") << endl;
1042  os <<"The number of iterations performed is " << this->iter_ << endl;
1043  os <<"The current block size is " << this->blockSize_ << endl;
1044  os <<"The number of auxiliary vectors is " << this->numAuxVecs_ << endl;
1045  os <<"The number of operations A*x is " << this->counterAOp_ << endl;
1046  os <<"The number of operations B*x is " << this->counterBOp_ << endl;
1047  os <<"The number of operations B*x by the orthomanager is " << this->orthman_->getOpCounter() << endl;
1048  os <<"The number of operations Prec*x is " << this->counterPrec_ << endl;
1049  os <<"Parameter rho_prime is " << rho_prime_ << endl;
1050  os <<"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
1051  os <<"Number of inner iterations was " << innerIters_ << endl;
1052  os <<"Total number of inner iterations is " << totalInnerIters_ << endl;
1053  os <<"f(x) is " << this->fx_ << endl;
1054 
1055  os.setf(std::ios_base::right, std::ios_base::adjustfield);
1056 
1057  if (this->initialized_) {
1058  os << endl;
1059  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1060  os << std::setw(20) << "Eigenvalue"
1061  << std::setw(20) << "Residual(B)"
1062  << std::setw(20) << "Residual(2)"
1063  << endl;
1064  os <<"--------------------------------------------------------------------------------"<<endl;
1065  for (int i=0; i<this->blockSize_; i++) {
1066  os << std::setw(20) << this->theta_[i];
1067  if (this->Rnorms_current_) os << std::setw(20) << this->Rnorms_[i];
1068  else os << std::setw(20) << "not current";
1069  if (this->R2norms_current_) os << std::setw(20) << this->R2norms_[i];
1070  else os << std::setw(20) << "not current";
1071  os << endl;
1072  }
1073  }
1074  os <<"================================================================================" << endl;
1075  os << endl;
1076  }
1077 
1078 
1079 } // end Anasazi namespace
1080 
1081 #endif // ANASAZI_IRTR_HPP
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
This class defines the interface required by an eigensolver and status test class to compute solution...
Base class for Implicit Riemannian Trust-Region solvers.
Declaration of basic traits for the multivector type.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers. The class provides the interfaces shared by the IRTR solvers (e.g., getState() and initialize()) as well as the shared implementations (e.g., inner products).
Virtual base class which defines basic traits for the operator type.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi&#39;s templated, static class providing utilities for the solvers.
OrdinalType numRows() const
Output managers remove the need for the eigensolver to know any information about the required output...
IRTR(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
IRTR constructor with eigenproblem, solver utilities, and parameter list of solver options...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
void iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
int reshape(OrdinalType numRows, OrdinalType numCols)
Types and exceptions used within Anasazi solvers and interfaces.
virtual ~IRTR()
IRTR destructor
Definition: AnasaziIRTR.hpp:72
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi&#39;s solvers.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
void currentStatus(std::ostream &os)
Impemements Eigensolver. This method requests that the solver print out its current status to screen...