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