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