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"
50 #include "AnasaziMatOrthoManager.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 {
104  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
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 
145  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
146  typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
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 {
207  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
208  typedef Anasazi::OperatorTraits<Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > OPT;
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 {
243  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
244  typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
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 {
278  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
279  typedef Anasazi::OperatorTraits<Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > OPT;
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
Anasazi::MultiVecTraits::MvTimesMatAddMv
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Definition: AnasaziMultiVecTraits.hpp:224
AnasaziTraceMinBase.hpp
Abstract base class for trace minimization eigensolvers.
AnasaziMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition: AnasaziOperatorTraits.hpp:84
Anasazi::MultiVecTraits::MvDot
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
Definition: AnasaziMultiVecTraits.hpp:254
Anasazi::OperatorTraits::Apply
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Definition: AnasaziOperatorTraits.hpp:92
Anasazi::MultiVecTraits::Assign
static void Assign(const MV &A, MV &mv)
mv := A
Definition: AnasaziMultiVecTraits.hpp:297
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp:127
Anasazi::MultiVecTraits::CloneCopy
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Definition: AnasaziMultiVecTraits.hpp:143
Anasazi::MultiVecTraits::MvTransMv
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Definition: AnasaziMultiVecTraits.hpp:248
Anasazi::MultiVecTraits::MvScale
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Definition: AnasaziMultiVecTraits.hpp:236
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Anasazi::MultiVecTraits::CloneViewNonConst
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).
Definition: AnasaziMultiVecTraits.hpp:171
Anasazi::MultiVecTraits::MvNorm
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Definition: AnasaziMultiVecTraits.hpp:264
Anasazi::MultiVecTraits::CloneView
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).
Definition: AnasaziMultiVecTraits.hpp:190
Anasazi::MultiVecTraits::SetBlock
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
Definition: AnasaziMultiVecTraits.hpp:276
Anasazi::MultiVecTraits::GetNumberVecs
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Definition: AnasaziMultiVecTraits.hpp:214
Anasazi::MultiVecTraits::MvAddMv
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Definition: AnasaziMultiVecTraits.hpp:231
Anasazi::MultiVecTraits::Clone
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Definition: AnasaziMultiVecTraits.hpp:136
Generated by   doxygen 1.8.13