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