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