Anasazi Version of the Day
Loading...
Searching...
No Matches
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"
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
72using Teuchos::RCP;
73
74namespace Anasazi {
75namespace Experimental {
76
77
78
81// Abstract base class for all operators
84template <class Scalar, class MV, class OP>
85class TraceMinOp
86{
87public:
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
101template <class Scalar, class MV, class OP>
102class TraceMinProjOp
103{
105 const Scalar ONE;
106
107public:
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
121private:
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
138template <class Scalar, class MV, class OP>
139class 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
150public:
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
181private:
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
204template <class Scalar, class MV, class OP>
205class TraceMinProjRitzOp : public TraceMinOp<Scalar,MV,OP>
206{
209
210public:
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
224private:
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?
240template <class Scalar, class MV, class OP>
241class TraceMinProjectedPrecOp : public TraceMinOp<Scalar,MV,OP>
242{
245 const Scalar ONE;
246
247public:
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
258private:
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
275template <class Scalar, class MV, class OP>
276class TraceMinProjRitzOpWithPrec : public TraceMinOp<Scalar,MV,OP>
277{
280 const Scalar ONE;
281
282public:
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
295private:
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
314namespace Anasazi
315{
316template <class Scalar, class MV, class OP>
317class 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
333namespace Belos
334{
335template <class Scalar, class MV, class OP>
336class 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
353namespace Anasazi
354{
355template <class Scalar, class MV, class OP>
356class 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
372namespace Anasazi
373{
374template <class Scalar, class MV, class OP>
375class 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
391namespace Anasazi
392{
393template <class Scalar, class MV, class OP>
394class 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
410namespace Anasazi {
411namespace Experimental {
414// TraceMinProjOp implementations
417template <class Scalar, class MV, class OP>
418TraceMinProjOp<Scalar,MV,OP>::TraceMinProjOp(const Teuchos::RCP<const MV> X, const Teuchos::RCP<const OP> B, Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
419ONE(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
459template <class Scalar, class MV, class OP>
460TraceMinProjOp<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) :
461ONE(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
490template <class Scalar, class MV, class OP>
491TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
492{
493 if(orthman_ != Teuchos::null)
494 orthman_->setOp(B_);
495}
496
497
498// Compute Px = x - proj proj'x
499template <class Scalar, class MV, class OP>
500void 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
524template <class Scalar, class MV, class OP>
525void 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
550template <class Scalar, class MV, class OP>
551TraceMinRitzOp<Scalar,MV,OP>::TraceMinRitzOp(const Teuchos::RCP<const OP>& A, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<const OP>& Prec) :
552ONE(Teuchos::ScalarTraits<Scalar>::one()),
553ZERO(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
579template <class Scalar, class MV, class OP>
580void 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
643template <class Scalar, class MV, class OP>
644void 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
668template <class Scalar, class MV, class OP>
669void 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
697template <class Scalar, class MV, class OP>
698TraceMinProjRitzOp<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
714template <class Scalar, class MV, class OP>
715TraceMinProjRitzOp<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
733template <class Scalar, class MV, class OP>
734void 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
752template <class Scalar, class MV, class OP>
753void 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
778template <class Scalar, class MV, class OP>
779void 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
792// TraceMinProjRitzOpWithPrec implementations
795template <class Scalar, class MV, class OP>
796TraceMinProjRitzOpWithPrec<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) :
797ONE(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
822template <class Scalar, class MV, class OP>
823TraceMinProjRitzOpWithPrec<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) :
824ONE(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
849template <class Scalar, class MV, class OP>
850void 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
859template <class Scalar, class MV, class OP>
860void 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
896template <class Scalar, class MV, class OP>
897void 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
911template <class Scalar, class MV, class OP>
912TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(const Teuchos::RCP<const OP> Op, const Teuchos::RCP<const MV> projVecs, Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
913ONE (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
955template <class Scalar, class MV, class OP>
956TraceMinProjectedPrecOp<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) :
957ONE(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
1027template <class Scalar, class MV, class OP>
1028TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
1029{
1030 if(orthman_ != Teuchos::null)
1031 orthman_->setOp(B_);
1032}
1033
1034
1035
1036template <class Scalar, class MV, class OP>
1037void 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
1068template <class Scalar, class MV, class OP>
1069void 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 header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Abstract base class for trace minimization eigensolvers.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
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.