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