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