Anasazi Version of the Day
Loading...
Searching...
No Matches
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
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 SIRTR : public RTRBase<ScalarType,MV,OP> {
87 public:
88
90
91
103 SIRTR( 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 ~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;
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
176
177
178
180 // Constructor
181 template <class ScalarType, class MV, class OP>
183 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
184 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
185 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
186 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
187 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
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;
236 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PSDM;
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
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.
This class defines the interface required by an eigensolver and status test class to compute solution...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Output managers remove the need for the eigensolver to know any information about the required output...
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...
virtual ~SIRTR()
SIRTR destructor
void currentStatus(std::ostream &os)
Impemements Eigensolver. This method requests that the solver print out its current status to screen.
void iterate()
Impemements Eigensolver. The outer IRTR iteration. See RTRBase::iterate().
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.
Anasazi's templated, static class providing utilities for the solvers.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
@ IterationDetails