Anasazi Version of the Day
Loading...
Searching...
No Matches
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
53#include "Teuchos_ScalarTraits.hpp"
54
55#include "Teuchos_LAPACK.hpp"
56#include "Teuchos_BLAS.hpp"
57#include "Teuchos_SerialDenseMatrix.hpp"
58#include "Teuchos_ParameterList.hpp"
59#include "Teuchos_TimeMonitor.hpp"
60
80// TODO: add randomization
81// TODO: add expensive debug checking on Teuchos_Debug
82
83namespace Anasazi {
84
85 template <class ScalarType, class MV, class OP>
86 class IRTR : public RTRBase<ScalarType,MV,OP> {
87 public:
88
90
91
103 IRTR( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
104 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
105 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
106 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
107 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
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;
138 typedef OperatorTraits<ScalarType,MV,OP> OPT;
139 typedef Teuchos::ScalarTraits<ScalarType> SCT;
140 typedef typename SCT::magnitudeType MagnitudeType;
141 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
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>
186 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
187 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
188 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
189 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
190 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
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;
241 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
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
781 Teuchos::SerialDenseMatrix<int,ScalarType> AA, BB, S;
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 ) ) {
817 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
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 << printMat(AA) << std::endl;;
893 this->om_->stream(Debug) << "BB: " << std::endl << printMat(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 << printMat(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
1048 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
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 )) {
1056 typename RTRBase<ScalarType,MV,OP>::CheckList chk;
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
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Base class for Implicit Riemannian Trust-Region solvers.
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.
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.
virtual ~IRTR()
IRTR destructor
void iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers....
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.