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
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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...
static void Assign(const MV &A, MV &mv)
mv := A
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 MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
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 Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
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 Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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 void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.