Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTraceMinRitzOp.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
14#ifndef TRACEMIN_RITZ_OP_HPP
15#define TRACEMIN_RITZ_OP_HPP
16
17#include "AnasaziConfigDefs.hpp"
18#include "AnasaziMatOrthoManager.hpp"
19#include "AnasaziMinres.hpp"
20#include "AnasaziTraceMinBase.hpp"
21
22#ifdef HAVE_ANASAZI_BELOS
23 #include "BelosMultiVecTraits.hpp"
24 #include "BelosLinearProblem.hpp"
25 #include "BelosPseudoBlockGmresSolMgr.hpp"
26 #include "BelosOperator.hpp"
27 #ifdef HAVE_ANASAZI_TPETRA
28 #include "BelosTpetraAdapter.hpp"
29 #endif
30 #ifdef HAVE_ANASAZI_EPETRA
31 #include "BelosEpetraAdapter.hpp"
32 #endif
33#endif
34
35#include "Teuchos_RCP.hpp"
36#include "Teuchos_SerialDenseSolver.hpp"
37#include "Teuchos_ScalarTraitsDecl.hpp"
38
39
40using Teuchos::RCP;
41
42namespace Anasazi {
43namespace Experimental {
44
45
46
49// Abstract base class for all operators
52template <class Scalar, class MV, class OP>
53class TraceMinOp
54{
55public:
56 virtual ~TraceMinOp() { };
57 virtual void Apply(const MV& X, MV& Y) const =0;
58 virtual void removeIndices(const std::vector<int>& indicesToRemove) =0;
59};
60
61
62
65// Defines a projector
66// Applies P_i to each individual vector x_i
69template <class Scalar, class MV, class OP>
70class TraceMinProjOp
71{
72 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
73 const Scalar ONE;
74
75public:
76 // Constructors
77 TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
78 TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs);
79
80 // Destructor
81 ~TraceMinProjOp();
82
83 // Applies the projector to a multivector
84 void Apply(const MV& X, MV& Y) const;
85
86 // Called by MINRES when certain vectors converge
87 void removeIndices(const std::vector<int>& indicesToRemove);
88
89private:
90 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
91 Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
92 Teuchos::RCP<const OP> B_;
93
94#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
95 RCP<Teuchos::Time> ProjTime_;
96#endif
97};
98
99
102// This class defines an operator A + \sigma B
103// This is used to apply shifts within TraceMin
106template <class Scalar, class MV, class OP>
107class TraceMinRitzOp : public TraceMinOp<Scalar,MV,OP>
108{
109 template <class Scalar_, class MV_, class OP_> friend class TraceMinProjRitzOp;
110 template <class Scalar_, class MV_, class OP_> friend class TraceMinProjRitzOpWithPrec;
111 template <class Scalar_, class MV_, class OP_> friend class TraceMinProjectedPrecOp;
112
113 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
114 typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
115 const Scalar ONE;
116 const Scalar ZERO;
117
118public:
119 // constructors for standard/generalized EVP
120 TraceMinRitzOp(const Teuchos::RCP<const OP>& A, const Teuchos::RCP<const OP>& B = Teuchos::null, const Teuchos::RCP<const OP>& Prec = Teuchos::null);
121
122 // Destructor
123 ~TraceMinRitzOp() { };
124
125 // sets the Ritz shift
126 void setRitzShifts(std::vector<Scalar> shifts) {ritzShifts_ = shifts;};
127
128 Scalar getRitzShift(const int subscript) { return ritzShifts_[subscript]; };
129
130 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > getPrec() { return Prec_; };
131
132 // sets the tolerances for the inner solves
133 void setInnerTol(const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
134
135 void getInnerTol(std::vector<Scalar>& tolerances) const { tolerances = tolerances_; };
136
137 void setMaxIts(const int maxits) { maxits_ = maxits; };
138
139 int getMaxIts() const { return maxits_; };
140
141 // applies A+\sigma B to a vector
142 void Apply(const MV& X, MV& Y) const;
143
144 // returns (A+\sigma B)\X
145 void ApplyInverse(const MV& X, MV& Y);
146
147 void removeIndices(const std::vector<int>& indicesToRemove);
148
149private:
150 Teuchos::RCP<const OP> A_, B_;
151 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Prec_;
152
153 int maxits_;
154 std::vector<Scalar> ritzShifts_;
155 std::vector<Scalar> tolerances_;
156
157 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > > solver_;
158
159#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
160 RCP<Teuchos::Time> PetraMultTime_, AopMultTime_;
161#endif
162};
163
164
165
168// Defines an operator P (A + \sigma B) P
169// Used for TraceMin with the projected iterative solver
172template <class Scalar, class MV, class OP>
173class TraceMinProjRitzOp : public TraceMinOp<Scalar,MV,OP>
174{
175 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
176 typedef Anasazi::OperatorTraits<Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > OPT;
177
178public:
179 // constructors for standard/generalized EVP
180 TraceMinProjRitzOp(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, const Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
181 TraceMinProjRitzOp(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, const Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs);
182
183 // applies P (A+\sigma B) P to a vector
184 void Apply(const MV& X, MV& Y) const;
185
186 // returns P(A+\sigma B)P\X
187 // this is not const due to the clumsiness with amSolving
188 void ApplyInverse(const MV& X, MV& Y);
189
190 void removeIndices(const std::vector<int>& indicesToRemove);
191
192private:
193 Teuchos::RCP< TraceMinRitzOp<Scalar,MV,OP> > Op_;
194 Teuchos::RCP< TraceMinProjOp<Scalar,MV,OP> > projector_;
195
196 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> > > solver_;
197};
198
199
200
203// Defines a preconditioner to be used with our projection method
204// Because we're using projected CG/minres/gmres, this preconditioner has to do projection as well
207// TODO: Should this be public?
208template <class Scalar, class MV, class OP>
209class TraceMinProjectedPrecOp : public TraceMinOp<Scalar,MV,OP>
210{
211 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
212 typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
213 const Scalar ONE;
214
215public:
216 // constructors for standard/generalized EVP
217 TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
218 TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs);
219
220 ~TraceMinProjectedPrecOp();
221
222 void Apply(const MV& X, MV& Y) const;
223
224 void removeIndices(const std::vector<int>& indicesToRemove);
225
226private:
227 Teuchos::RCP<const OP> Op_;
228 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
229
230 Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
231 Teuchos::RCP<const OP> B_;
232};
233
234
235
238// Defines a preconditioner to be used with our projection method
239// Because we're using projected CG/minres/gmres, this preconditioner has to do projection as well
242#ifdef HAVE_ANASAZI_BELOS
243template <class Scalar, class MV, class OP>
244class TraceMinProjRitzOpWithPrec : public TraceMinOp<Scalar,MV,OP>
245{
246 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
247 typedef Anasazi::OperatorTraits<Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > OPT;
248 const Scalar ONE;
249
250public:
251 // constructors for standard/generalized EVP
252 TraceMinProjRitzOpWithPrec(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
253 TraceMinProjRitzOpWithPrec(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs);
254
255 void Apply(const MV& X, MV& Y) const;
256
257 // returns P(A+\sigma B)P\X
258 // this is not const due to the clumsiness with amSolving
259 void ApplyInverse(const MV& X, MV& Y);
260
261 void removeIndices(const std::vector<int>& indicesToRemove);
262
263private:
264 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Op_;
265 Teuchos::RCP<TraceMinProjectedPrecOp<Scalar,MV,OP> > Prec_;
266
267 Teuchos::RCP< Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> > > solver_;
268 Teuchos::RCP< Belos::LinearProblem<Scalar,MV,TraceMinOp<Scalar,MV,OP> > > problem_;
269};
270#endif
271
272}} // end of namespace
273
274
275
278// Operator traits classes
279// Required to use user-defined operators with a Krylov solver in Belos
282namespace Anasazi
283{
284template <class Scalar, class MV, class OP>
285class OperatorTraits <Scalar, MV, Experimental::TraceMinOp<Scalar,MV,OP> >
286 {
287 public:
288 static void
289 Apply (const Experimental::TraceMinOp<Scalar,MV,OP>& Op,
290 const MV& x,
291 MV& y) {Op.Apply(x,y);};
292
294 static bool
295 HasApplyTranspose (const Experimental::TraceMinOp<Scalar,MV,OP>& Op) {return true;};
296 };
297}
298
299
300#ifdef HAVE_ANASAZI_BELOS
301namespace Belos
302{
303template <class Scalar, class MV, class OP>
304class OperatorTraits <Scalar, MV, Anasazi::Experimental::TraceMinOp<Scalar,MV,OP> >
305 {
306 public:
307 static void
308 Apply (const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op,
309 const MV& x,
310 MV& y, Belos::ETrans trans = NOTRANS) {Op.Apply(x,y);};
311
313 static bool
314 HasApplyTranspose (const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op) {return true;};
315 };
316}
317#endif
318
319
320
321namespace Anasazi
322{
323template <class Scalar, class MV, class OP>
324class OperatorTraits <Scalar, MV, Experimental::TraceMinRitzOp<Scalar,MV,OP> >
325 {
326 public:
327 static void
328 Apply (const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op,
329 const MV& x,
330 MV& y) {Op.Apply(x,y);};
331
333 static bool
334 HasApplyTranspose (const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op) {return true;};
335 };
336}
337
338
339
340namespace Anasazi
341{
342template <class Scalar, class MV, class OP>
343class OperatorTraits <Scalar, MV, Experimental::TraceMinProjRitzOp<Scalar,MV,OP> >
344 {
345 public:
346 static void
347 Apply (const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op,
348 const MV& x,
349 MV& y) {Op.Apply(x,y);};
350
352 static bool
353 HasApplyTranspose (const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op) {return true;};
354 };
355}
356
357
358
359namespace Anasazi
360{
361template <class Scalar, class MV, class OP>
362class OperatorTraits <Scalar, MV, Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP> >
363 {
364 public:
365 static void
366 Apply (const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op,
367 const MV& x,
368 MV& y) {Op.Apply(x,y);};
369
371 static bool
372 HasApplyTranspose (const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op) {return true;};
373 };
374}
375
376
377
378namespace Anasazi {
379namespace Experimental {
382// TraceMinProjOp implementations
385template <class Scalar, class MV, class OP>
386TraceMinProjOp<Scalar,MV,OP>::TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
387ONE(Teuchos::ScalarTraits<Scalar>::one())
388{
389#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
390 ProjTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinProjOp::Apply()");
391#endif
392
393 B_ = B;
394 orthman_ = orthman;
395 if(orthman_ != Teuchos::null && B_ != Teuchos::null)
396 orthman_->setOp(Teuchos::null);
397
398 // Make it so X'BBX = I
399 // If there is no B, this step is unnecessary because X'X = I already
400 if(B_ != Teuchos::null)
401 {
402 int nvec = MVT::GetNumberVecs(*X);
403
404 if(orthman_ != Teuchos::null)
405 {
406 // Want: X'X = I NOT X'MX = I
407 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*X);
408 orthman_->normalizeMat(*helperMV);
409 projVecs_.push_back(helperMV);
410 }
411 else
412 {
413 std::vector<Scalar> normvec(nvec);
414 MVT::MvNorm(*X,normvec);
415 for(int i=0; i<nvec; i++)
416 normvec[i] = ONE/normvec[i];
417 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*X);
418 MVT::MvScale(*helperMV,normvec);
419 projVecs_.push_back(helperMV);
420 }
421 }
422 else
423 projVecs_.push_back(MVT::CloneCopy(*X));
424}
425
426
427template <class Scalar, class MV, class OP>
428TraceMinProjOp<Scalar,MV,OP>::TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs) :
429ONE(Teuchos::ScalarTraits<Scalar>::one())
430{
431#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
432 ProjTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinProjOp::Apply()");
433#endif
434
435 B_ = B;
436 orthman_ = orthman;
437 if(B_ != Teuchos::null)
438 orthman_->setOp(Teuchos::null);
439
440 projVecs_ = auxVecs;
441
442 // Make it so X'BBX = I
443 // If there is no B, this step is unnecessary because X'X = I already
444 if(B_ != Teuchos::null)
445 {
446 // Want: X'X = I NOT X'MX = I
447 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*X);
448 orthman_->normalizeMat(*helperMV);
449 projVecs_.push_back(helperMV);
450
451 }
452 else
453 projVecs_.push_back(MVT::CloneCopy(*X));
454}
455
456
457// Destructor - make sure to reset the operator in the ortho manager
458template <class Scalar, class MV, class OP>
459TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
460{
461 if(orthman_ != Teuchos::null)
462 orthman_->setOp(B_);
463}
464
465
466// Compute Px = x - proj proj'x
467template <class Scalar, class MV, class OP>
468void TraceMinProjOp<Scalar,MV,OP>::Apply(const MV& X, MV& Y) const
469{
470#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
471 Teuchos::TimeMonitor lcltimer( *ProjTime_ );
472#endif
473
474 if(orthman_ != Teuchos::null)
475 {
476 MVT::Assign(X,Y);
477 orthman_->projectMat(Y,projVecs_);
478 }
479 else
480 {
481 int nvec = MVT::GetNumberVecs(X);
482 std::vector<Scalar> dotProducts(nvec);
483 MVT::MvDot(*projVecs_[0],X,dotProducts);
484 Teuchos::RCP<MV> helper = MVT::CloneCopy(*projVecs_[0]);
485 MVT::MvScale(*helper,dotProducts);
486 MVT::MvAddMv(ONE,X,-ONE,*helper,Y);
487 }
488}
489
490
491
492template <class Scalar, class MV, class OP>
493void TraceMinProjOp<Scalar,MV,OP>::removeIndices(const std::vector<int>& indicesToRemove)
494{
495 if (orthman_ == Teuchos::null) {
496 const int nprojvecs = projVecs_.size();
497 const int nvecs = MVT::GetNumberVecs(*projVecs_[nprojvecs-1]);
498 const int numRemoving = indicesToRemove.size();
499 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
500
501 for (int i=0; i<nvecs; i++) {
502 helper[i] = i;
503 }
504
505 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
506
507 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*projVecs_[nprojvecs-1],indicesToLeave);
508 projVecs_[nprojvecs-1] = helperMV;
509 }
510}
511
512
515// TraceMinRitzOp implementations
518template <class Scalar, class MV, class OP>
519TraceMinRitzOp<Scalar,MV,OP>::TraceMinRitzOp(const Teuchos::RCP<const OP>& A, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<const OP>& Prec) :
520ONE(Teuchos::ScalarTraits<Scalar>::one()),
521ZERO(Teuchos::ScalarTraits<Scalar>::zero())
522{
523 A_ = A;
524 B_ = B;
525 // TODO: maxits should not be hard coded
526 maxits_ = 200;
527
528#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
529 PetraMultTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinRitzOp: *Petra::Apply()");
530 AopMultTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinRitzOp::Apply()");
531#endif
532
533 // create the operator for my minres solver
534 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(this);
535 linSolOp.release();
536
537 // TODO: This should support left and right prec
538 if(Prec != Teuchos::null)
539 Prec_ = Teuchos::rcp( new TraceMinRitzOp<Scalar,MV,OP>(Prec) );
540
541 // create my minres solver
542 solver_ = Teuchos::rcp( new PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> >(linSolOp,Prec_) );
543}
544
545
546
547template <class Scalar, class MV, class OP>
548void TraceMinRitzOp<Scalar,MV,OP>::Apply(const MV& X, MV& Y) const
549{
550 int nvecs = MVT::GetNumberVecs(X);
551
552#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
553 Teuchos::TimeMonitor outertimer( *AopMultTime_ );
554#endif
555
556 // Y := A*X
557 {
558#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
559 Teuchos::TimeMonitor lcltimer( *PetraMultTime_ );
560#endif
561
562 OPT::Apply(*A_,X,Y);
563 }
564
565 // If we are a preconditioner, we're not using shifts
566 if(ritzShifts_.size() > 0)
567 {
568 // Get the indices of nonzero Ritz shifts
569 std::vector<int> nonzeroRitzIndices;
570 nonzeroRitzIndices.reserve(nvecs);
571 for(int i=0; i<nvecs; i++)
572 {
573 if(ritzShifts_[i] != ZERO)
574 nonzeroRitzIndices.push_back(i);
575 }
576
577 // Handle Ritz shifts
578 int numRitzShifts = nonzeroRitzIndices.size();
579 if(numRitzShifts > 0)
580 {
581 // Get pointers to the appropriate parts of X and Y
582 Teuchos::RCP<const MV> ritzX = MVT::CloneView(X,nonzeroRitzIndices);
583 Teuchos::RCP<MV> ritzY = MVT::CloneViewNonConst(Y,nonzeroRitzIndices);
584
585 // Get the nonzero Ritz shifts
586 std::vector<Scalar> nonzeroRitzShifts(numRitzShifts);
587 for(int i=0; i<numRitzShifts; i++)
588 nonzeroRitzShifts[i] = ritzShifts_[nonzeroRitzIndices[i]];
589
590 // Compute Y := AX + ritzShift BX
591 if(B_ != Teuchos::null)
592 {
593 Teuchos::RCP<MV> BX = MVT::Clone(*ritzX,numRitzShifts);
594 OPT::Apply(*B_,*ritzX,*BX);
595 MVT::MvScale(*BX,nonzeroRitzShifts);
596 MVT::MvAddMv(ONE,*ritzY,-ONE,*BX,*ritzY);
597 }
598 // Compute Y := AX + ritzShift X
599 else
600 {
601 Teuchos::RCP<MV> scaledX = MVT::CloneCopy(*ritzX);
602 MVT::MvScale(*scaledX,nonzeroRitzShifts);
603 MVT::MvAddMv(ONE,*ritzY,-ONE,*scaledX,*ritzY);
604 }
605 }
606 }
607}
608
609
610
611template <class Scalar, class MV, class OP>
612void TraceMinRitzOp<Scalar,MV,OP>::ApplyInverse(const MV& X, MV& Y)
613{
614 int nvecs = MVT::GetNumberVecs(X);
615 std::vector<int> indices(nvecs);
616 for(int i=0; i<nvecs; i++)
617 indices[i] = i;
618
619 Teuchos::RCP<const MV> rcpX = MVT::CloneView(X,indices);
620 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
621
622 // Solve the linear system A*Y = X
623 solver_->setTol(tolerances_);
624 solver_->setMaxIter(maxits_);
625
626 // Set solution and RHS
627 solver_->setSol(rcpY);
628 solver_->setRHS(rcpX);
629
630 // Solve the linear system
631 solver_->solve();
632}
633
634
635
636template <class Scalar, class MV, class OP>
637void TraceMinRitzOp<Scalar,MV,OP>::removeIndices(const std::vector<int>& indicesToRemove)
638{
639 int nvecs = tolerances_.size();
640 int numRemoving = indicesToRemove.size();
641 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
642 std::vector<Scalar> helperS(nvecs-numRemoving);
643
644 for(int i=0; i<nvecs; i++)
645 helper[i] = i;
646
647 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
648
649 for(int i=0; i<nvecs-numRemoving; i++)
650 helperS[i] = ritzShifts_[indicesToLeave[i]];
651 ritzShifts_ = helperS;
652
653 for(int i=0; i<nvecs-numRemoving; i++)
654 helperS[i] = tolerances_[indicesToLeave[i]];
655 tolerances_ = helperS;
656}
657
658
659
662// TraceMinProjRitzOp implementations
665template <class Scalar, class MV, class OP>
666TraceMinProjRitzOp<Scalar,MV,OP>::TraceMinProjRitzOp(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, const Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman)
667{
668 Op_ = Op;
669
670 // Create the projector object
671 projector_ = Teuchos::rcp( new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman) );
672
673 // create the operator for my minres solver
674 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(this);
675 linSolOp.release();
676
677 // create my minres solver
678 solver_ = Teuchos::rcp( new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
679}
680
681
682template <class Scalar, class MV, class OP>
683TraceMinProjRitzOp<Scalar,MV,OP>::TraceMinProjRitzOp(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, const Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs)
684{
685 Op_ = Op;
686
687 // Create the projector object
688 projector_ = Teuchos::rcp( new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman, auxVecs) );
689
690 // create the operator for my minres solver
691 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(this);
692 linSolOp.release();
693
694 // create my minres solver
695 solver_ = Teuchos::rcp( new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
696}
697
698
699
700// Y:= P (A+\sigma B) P X
701template <class Scalar, class MV, class OP>
702void TraceMinProjRitzOp<Scalar,MV,OP>::Apply(const MV& X, MV& Y) const
703{
704 int nvecs = MVT::GetNumberVecs(X);
705
706// // compute PX
707// Teuchos::RCP<MV> PX = MVT::Clone(X,nvecs);
708// projector_->Apply(X,*PX);
709
710 // compute (A+\sigma B) P X
711 Teuchos::RCP<MV> APX = MVT::Clone(X,nvecs);
712 OPT::Apply(*Op_,X,*APX);
713
714 // compute Y := P (A+\sigma B) P X
715 projector_->Apply(*APX,Y);
716}
717
718
719
720template <class Scalar, class MV, class OP>
721void TraceMinProjRitzOp<Scalar,MV,OP>::ApplyInverse(const MV& X, MV& Y)
722{
723 int nvecs = MVT::GetNumberVecs(X);
724 std::vector<int> indices(nvecs);
725 for(int i=0; i<nvecs; i++)
726 indices[i] = i;
727
728 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
729 Teuchos::RCP<MV> PX = MVT::Clone(X,nvecs);
730 projector_->Apply(X,*PX);
731
732 // Solve the linear system A*Y = X
733 solver_->setTol(Op_->tolerances_);
734 solver_->setMaxIter(Op_->maxits_);
735
736 // Set solution and RHS
737 solver_->setSol(rcpY);
738 solver_->setRHS(PX);
739
740 // Solve the linear system
741 solver_->solve();
742}
743
744
745
746template <class Scalar, class MV, class OP>
747void TraceMinProjRitzOp<Scalar,MV,OP>::removeIndices(const std::vector<int>& indicesToRemove)
748{
749 Op_->removeIndices(indicesToRemove);
750
751 projector_->removeIndices(indicesToRemove);
752}
753
754
755
756
757#ifdef HAVE_ANASAZI_BELOS
760// TraceMinProjRitzOpWithPrec implementations
763template <class Scalar, class MV, class OP>
764TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::TraceMinProjRitzOpWithPrec(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
765ONE(Teuchos::ScalarTraits<Scalar>::one())
766{
767 Op_ = Op;
768
769 // create the operator for the Belos solver
770 Teuchos::RCP<TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(this);
771 linSolOp.release();
772
773 // Create the linear problem
774 problem_ = Teuchos::rcp(new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
775
776 // Set the operator
777 problem_->setOperator(linSolOp);
778
779 // Set the preconditioner
780 // TODO: This does not support right preconditioning
781 Prec_ = Teuchos::rcp( new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman) );
782// problem_->setLeftPrec(Prec_);
783
784 // create the pseudoblock gmres solver
785 // minres has trouble with the projected preconditioner
786 solver_ = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
787}
788
789
790template <class Scalar, class MV, class OP>
791TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::TraceMinProjRitzOpWithPrec(const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs) :
792ONE(Teuchos::ScalarTraits<Scalar>::one())
793{
794 Op_ = Op;
795
796 // create the operator for the Belos solver
797 Teuchos::RCP<const TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(this);
798 linSolOp.release();
799
800 // Create the linear problem
801 problem_ = Teuchos::rcp(new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
802
803 // Set the operator
804 problem_->setOperator(linSolOp);
805
806 // Set the preconditioner
807 // TODO: This does not support right preconditioning
808 Prec_ = Teuchos::rcp( new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman, auxVecs) );
809// problem_->setLeftPrec(Prec_);
810
811 // create the pseudoblock gmres solver
812 // minres has trouble with the projected preconditioner
813 solver_ = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
814}
815
816
817template <class Scalar, class MV, class OP>
818void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::Apply(const MV& X, MV& Y) const
819{
820 int nvecs = MVT::GetNumberVecs(X);
821 RCP<MV> Ydot = MVT::Clone(Y,nvecs);
822 OPT::Apply(*Op_,X,*Ydot);
823 Prec_->Apply(*Ydot,Y);
824}
825
826
827template <class Scalar, class MV, class OP>
828void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::ApplyInverse(const MV& X, MV& Y)
829{
830 int nvecs = MVT::GetNumberVecs(X);
831 std::vector<int> indices(nvecs);
832 for(int i=0; i<nvecs; i++)
833 indices[i] = i;
834
835 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
836 Teuchos::RCP<MV> rcpX = MVT::Clone(X,nvecs);
837
838 Prec_->Apply(X,*rcpX);
839
840 // Create the linear problem
841 problem_->setProblem(rcpY,rcpX);
842
843 // Set the problem for the solver
844 solver_->setProblem(problem_);
845
846 // Set up the parameters for gmres
847 // TODO: Accept maximum number of iterations
848 // TODO: Make sure single shift really means single shift
849 // TODO: Look into fixing my problem with the deflation quorum
850 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(new Teuchos::ParameterList());
851 pl->set("Convergence Tolerance", Op_->tolerances_[0]);
852 pl->set("Block Size", nvecs);
853// pl->set("Verbosity", Belos::IterationDetails + Belos::StatusTestDetails + Belos::Debug);
854// pl->set("Output Frequency", 1);
855 pl->set("Maximum Iterations", Op_->getMaxIts());
856 pl->set("Num Blocks", Op_->getMaxIts());
857 solver_->setParameters(pl);
858
859 // Solve the linear system
860 solver_->solve();
861}
862
863
864template <class Scalar, class MV, class OP>
865void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::removeIndices(const std::vector<int>& indicesToRemove)
866{
867 Op_->removeIndices(indicesToRemove);
868
869 Prec_->removeIndices(indicesToRemove);
870}
871#endif
872
873
876// TraceMinProjectedPrecOp implementations
879template <class Scalar, class MV, class OP>
880TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
881ONE (Teuchos::ScalarTraits<Scalar>::one())
882{
883 Op_ = Op;
884 orthman_ = orthman;
885
886 int nvecs = MVT::GetNumberVecs(*projVecs);
887 Teuchos::RCP<MV> helperMV = MVT::Clone(*projVecs,nvecs);
888
889 // Compute Prec \ projVecs
890 OPT::Apply(*Op_,*projVecs,*helperMV);
891
892 if(orthman_ != Teuchos::null)
893 {
894 // Set the operator for the inner products
895 B_ = orthman_->getOp();
896 orthman_->setOp(Op_);
897
898 Teuchos::RCP<MV> locProjVecs = MVT::CloneCopy(*projVecs);
899
900 // Normalize the vectors such that Y' Prec \ Y = I
901 const int rank = orthman_->normalizeMat (*locProjVecs, Teuchos::null, helperMV);
902
903 // FIXME (mfh 08 Aug 2014) Write a named exception for this case.
904 TEUCHOS_TEST_FOR_EXCEPTION(rank != nvecs, std::runtime_error, "Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
905
906 orthman_->setOp(Teuchos::null);
907 }
908 else
909 {
910 std::vector<Scalar> dotprods(nvecs);
911 MVT::MvDot(*projVecs,*helperMV,dotprods);
912
913 for(int i=0; i<nvecs; i++)
914 dotprods[i] = ONE/sqrt(dotprods[i]);
915
916 MVT::MvScale(*helperMV, dotprods);
917 }
918
919 projVecs_.push_back(helperMV);
920}
921
922
923template <class Scalar, class MV, class OP>
924TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs) :
925ONE(Teuchos::ScalarTraits<Scalar>::one())
926{
927 Op_ = Op;
928 orthman_ = orthman;
929
930 int nvecs;
931 Teuchos::RCP<MV> locProjVecs;
932
933 // Add the aux vecs to the projector
934 if(auxVecs.size() > 0)
935 {
936 // Get the total number of vectors
937 nvecs = MVT::GetNumberVecs(*projVecs);
938 for(int i=0; i<auxVecs.size(); i++)
939 nvecs += MVT::GetNumberVecs(*auxVecs[i]);
940
941 // Allocate space for all of them
942 locProjVecs = MVT::Clone(*projVecs, nvecs);
943
944 // Copy the vectors over
945 int startIndex = 0;
946 std::vector<int> locind(nvecs);
947
948 locind.resize(MVT::GetNumberVecs(*projVecs));
949 for (size_t i = 0; i<locind.size(); i++) {
950 locind[i] = startIndex + i;
951 }
952 startIndex += locind.size();
953 MVT::SetBlock(*projVecs,locind,*locProjVecs);
954
955 for (size_t i=0; i < static_cast<size_t> (auxVecs.size ()); ++i)
956 {
957 locind.resize(MVT::GetNumberVecs(*auxVecs[i]));
958 for(size_t j=0; j<locind.size(); j++) locind[j] = startIndex + j;
959 startIndex += locind.size();
960 MVT::SetBlock(*auxVecs[i],locind,*locProjVecs);
961 }
962 }
963 else
964 {
965 // Copy the vectors over
966 nvecs = MVT::GetNumberVecs(*projVecs);
967 locProjVecs = MVT::CloneCopy(*projVecs);
968 }
969
970 Teuchos::RCP<MV> helperMV = MVT::Clone(*projVecs,nvecs);
971
972 // Compute Prec \ projVecs
973 OPT::Apply(*Op_,*locProjVecs,*helperMV);
974
975 // Set the operator for the inner products
976 B_ = orthman_->getOp();
977 orthman_->setOp(Op_);
978
979 // Normalize the vectors such that Y' Prec \ Y = I
980 const int rank = orthman_->normalizeMat(*locProjVecs,Teuchos::null,helperMV);
981
982 projVecs_.push_back(helperMV);
983
984// helperMV->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
985
986 // FIXME (mfh 08 Aug 2014) Write a named exception for this case.
987 TEUCHOS_TEST_FOR_EXCEPTION(
988 rank != nvecs, std::runtime_error,
989 "Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
990
991 orthman_->setOp(Teuchos::null);
992}
993
994
995template <class Scalar, class MV, class OP>
996TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
997{
998 if(orthman_ != Teuchos::null)
999 orthman_->setOp(B_);
1000}
1001
1002
1003
1004template <class Scalar, class MV, class OP>
1005void TraceMinProjectedPrecOp<Scalar,MV,OP>::Apply(const MV& X, MV& Y) const
1006{
1007 int nvecsX = MVT::GetNumberVecs(X);
1008
1009 if(orthman_ != Teuchos::null)
1010 {
1011 // Y = M\X - proj proj' X
1012 int nvecsP = MVT::GetNumberVecs(*projVecs_[0]);
1013 OPT::Apply(*Op_,X,Y);
1014
1015 Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> > projX = Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,Scalar>(nvecsP,nvecsX));
1016
1017 MVT::MvTransMv(ONE, *projVecs_[0], X, *projX);
1018
1019 MVT::MvTimesMatAddMv(-ONE, *projVecs_[0], *projX, ONE, Y);
1020 }
1021 else
1022 {
1023 Teuchos::RCP<MV> MX = MVT::Clone(X,nvecsX);
1024 OPT::Apply(*Op_,X,*MX);
1025
1026 std::vector<Scalar> dotprods(nvecsX);
1027 MVT::MvDot(*projVecs_[0], X, dotprods);
1028
1029 Teuchos::RCP<MV> helper = MVT::CloneCopy(*projVecs_[0]);
1030 MVT::MvScale(*helper, dotprods);
1031 MVT::MvAddMv(ONE, *MX, -ONE, *helper, Y);
1032 }
1033}
1034
1035
1036template <class Scalar, class MV, class OP>
1037void TraceMinProjectedPrecOp<Scalar,MV,OP>::removeIndices(const std::vector<int>& indicesToRemove)
1038{
1039 if(orthman_ == Teuchos::null)
1040 {
1041 int nvecs = MVT::GetNumberVecs(*projVecs_[0]);
1042 int numRemoving = indicesToRemove.size();
1043 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
1044
1045 for(int i=0; i<nvecs; i++)
1046 helper[i] = i;
1047
1048 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
1049
1050 Teuchos::RCP<const MV> helperMV = MVT::CloneCopy(*projVecs_[0],indicesToLeave);
1051 projVecs_[0] = helperMV;
1052 }
1053}
1054
1055}} // end of namespace
1056
1057#endif
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
AnasaziMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
AnasaziTraceMinBase.hpp
Abstract base class for trace minimization eigensolvers.
Anasazi::MatOrthoManager
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Definition AnasaziMatOrthoManager.hpp:44
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition AnasaziMultiVecTraits.hpp:95
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition AnasaziOperatorTraits.hpp:53
Anasazi::OperatorTraits::Apply
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Definition AnasaziOperatorTraits.hpp:60
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Generated by doxygen 1.9.8