Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
40 using Teuchos::RCP;
41 
42 namespace Anasazi {
43 namespace Experimental {
44 
45 
46 
49 // Abstract base class for all operators
52 template <class Scalar, class MV, class OP>
53 class TraceMinOp
54 {
55 public:
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
69 template <class Scalar, class MV, class OP>
70 class TraceMinProjOp
71 {
72  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
73  const Scalar ONE;
74 
75 public:
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 
89 private:
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
106 template <class Scalar, class MV, class OP>
107 class 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 
118 public:
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 
149 private:
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
172 template <class Scalar, class MV, class OP>
173 class 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 
178 public:
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 
192 private:
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?
208 template <class Scalar, class MV, class OP>
209 class 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 
215 public:
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 
226 private:
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
243 template <class Scalar, class MV, class OP>
244 class 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 
250 public:
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 
263 private:
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
282 namespace Anasazi
283 {
284 template <class Scalar, class MV, class OP>
285 class 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
301 namespace Belos
302 {
303 template <class Scalar, class MV, class OP>
304 class 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 
321 namespace Anasazi
322 {
323 template <class Scalar, class MV, class OP>
324 class 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 
340 namespace Anasazi
341 {
342 template <class Scalar, class MV, class OP>
343 class 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 
359 namespace Anasazi
360 {
361 template <class Scalar, class MV, class OP>
362 class 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 
378 namespace Anasazi {
379 namespace Experimental {
382 // TraceMinProjOp implementations
385 template <class Scalar, class MV, class OP>
386 TraceMinProjOp<Scalar,MV,OP>::TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
387 ONE(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 
427 template <class Scalar, class MV, class OP>
428 TraceMinProjOp<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) :
429 ONE(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
458 template <class Scalar, class MV, class OP>
459 TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
460 {
461  if(orthman_ != Teuchos::null)
462  orthman_->setOp(B_);
463 }
464 
465 
466 // Compute Px = x - proj proj'x
467 template <class Scalar, class MV, class OP>
468 void 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 
492 template <class Scalar, class MV, class OP>
493 void 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
518 template <class Scalar, class MV, class OP>
519 TraceMinRitzOp<Scalar,MV,OP>::TraceMinRitzOp(const Teuchos::RCP<const OP>& A, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<const OP>& Prec) :
520 ONE(Teuchos::ScalarTraits<Scalar>::one()),
521 ZERO(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 
547 template <class Scalar, class MV, class OP>
548 void 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 
611 template <class Scalar, class MV, class OP>
612 void 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 
636 template <class Scalar, class MV, class OP>
637 void 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
665 template <class Scalar, class MV, class OP>
666 TraceMinProjRitzOp<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 
682 template <class Scalar, class MV, class OP>
683 TraceMinProjRitzOp<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
701 template <class Scalar, class MV, class OP>
702 void 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 
720 template <class Scalar, class MV, class OP>
721 void 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 
746 template <class Scalar, class MV, class OP>
747 void 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
758 // TraceMinProjRitzOpWithPrec implementations
763 template <class Scalar, class MV, class OP>
764 TraceMinProjRitzOpWithPrec<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) :
765 ONE(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 
790 template <class Scalar, class MV, class OP>
791 TraceMinProjRitzOpWithPrec<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) :
792 ONE(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 
817 template <class Scalar, class MV, class OP>
818 void 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 
827 template <class Scalar, class MV, class OP>
828 void 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 
864 template <class Scalar, class MV, class OP>
865 void 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
879 template <class Scalar, class MV, class OP>
880 TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
881 ONE (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 
923 template <class Scalar, class MV, class OP>
924 TraceMinProjectedPrecOp<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) :
925 ONE(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 
995 template <class Scalar, class MV, class OP>
996 TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
997 {
998  if(orthman_ != Teuchos::null)
999  orthman_->setOp(B_);
1000 }
1001 
1002 
1003 
1004 template <class Scalar, class MV, class OP>
1005 void 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 
1036 template <class Scalar, class MV, class OP>
1037 void 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
AnasaziTraceMinBase.hpp
Abstract base class for trace minimization eigensolvers.
AnasaziMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition: AnasaziOperatorTraits.hpp:52
Teuchos_ScalarTraitsDecl.hpp
Anasazi::MultiVecTraits::MvDot
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
Definition: AnasaziMultiVecTraits.hpp:222
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
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
Teuchos_SerialDenseSolver.hpp
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp:95
Anasazi::MultiVecTraits::CloneCopy
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Definition: AnasaziMultiVecTraits.hpp:111
Anasazi::MultiVecTraits::MvScale
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Definition: AnasaziMultiVecTraits.hpp:204
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Anasazi::MultiVecTraits::MvNorm
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Definition: AnasaziMultiVecTraits.hpp:232
Teuchos::RCP::release
Ptr< T > release()
Teuchos::ParameterList::setParameters
ParameterList & setParameters(const ParameterList &source)
Anasazi::MultiVecTraits::SetBlock
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
Definition: AnasaziMultiVecTraits.hpp:244
Anasazi::MultiVecTraits::GetNumberVecs
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Definition: AnasaziMultiVecTraits.hpp:182
Teuchos_RCP.hpp
Anasazi::MultiVecTraits::Clone
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Definition: AnasaziMultiVecTraits.hpp:104
Generated on Fri Jul 26 2024 09:18:32 for Anasazi by   doxygen 1.8.5