33 #ifndef TRACEMIN_RITZ_OP_HPP 34 #define TRACEMIN_RITZ_OP_HPP 41 #ifdef HAVE_ANASAZI_BELOS 42 #include "BelosMultiVecTraits.hpp" 43 #include "BelosLinearProblem.hpp" 44 #include "BelosPseudoBlockGmresSolMgr.hpp" 45 #include "BelosOperator.hpp" 46 #ifdef HAVE_ANASAZI_TPETRA 47 #include "BelosTpetraAdapter.hpp" 49 #ifdef HAVE_ANASAZI_EPETRA 50 #include "BelosEpetraAdapter.hpp" 54 #include "Teuchos_RCP.hpp" 55 #include "Teuchos_SerialDenseSolver.hpp" 56 #include "Teuchos_ScalarTraitsDecl.hpp" 71 template <
class Scalar,
class MV,
class OP>
75 virtual ~TraceMinOp() { };
76 virtual void Apply(
const MV& X, MV& Y)
const =0;
77 virtual void removeIndices(
const std::vector<int>& indicesToRemove) =0;
88 template <
class Scalar,
class MV,
class OP>
97 TraceMinProjOp(
const Teuchos::RCP<const MV> X,
const Teuchos::RCP<const OP> B, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs);
103 void Apply(
const MV& X, MV& Y)
const;
106 void removeIndices(
const std::vector<int>& indicesToRemove);
109 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
110 Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
111 Teuchos::RCP<const OP> B_;
113 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 114 RCP<Teuchos::Time> ProjTime_;
125 template <
class Scalar,
class MV,
class OP>
126 class TraceMinRitzOp :
public TraceMinOp<Scalar,MV,OP>
128 template <
class Scalar_,
class MV_,
class OP_>
friend class TraceMinProjRitzOp;
129 template <
class Scalar_,
class MV_,
class OP_>
friend class TraceMinProjRitzOpWithPrec;
130 template <
class Scalar_,
class MV_,
class OP_>
friend class TraceMinProjectedPrecOp;
139 TraceMinRitzOp(
const Teuchos::RCP<const OP>& A,
const Teuchos::RCP<const OP>& B = Teuchos::null,
const Teuchos::RCP<const OP>& Prec = Teuchos::null);
142 ~TraceMinRitzOp() { };
145 void setRitzShifts(std::vector<Scalar> shifts) {ritzShifts_ = shifts;};
147 Scalar getRitzShift(
const int subscript) {
return ritzShifts_[subscript]; };
149 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > getPrec() {
return Prec_; };
152 void setInnerTol(
const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
154 void getInnerTol(std::vector<Scalar>& tolerances)
const { tolerances = tolerances_; };
156 void setMaxIts(
const int maxits) { maxits_ = maxits; };
158 int getMaxIts()
const {
return maxits_; };
161 void Apply(
const MV& X, MV& Y)
const;
164 void ApplyInverse(
const MV& X, MV& Y);
166 void removeIndices(
const std::vector<int>& indicesToRemove);
169 Teuchos::RCP<const OP> A_, B_;
170 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Prec_;
173 std::vector<Scalar> ritzShifts_;
174 std::vector<Scalar> tolerances_;
176 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > > solver_;
178 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 179 RCP<Teuchos::Time> PetraMultTime_, AopMultTime_;
191 template <
class Scalar,
class MV,
class OP>
192 class TraceMinProjRitzOp :
public TraceMinOp<Scalar,MV,OP>
199 TraceMinProjRitzOp(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs,
const Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
200 TraceMinProjRitzOp(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs,
const Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs);
203 void Apply(
const MV& X, MV& Y)
const;
207 void ApplyInverse(
const MV& X, MV& Y);
209 void removeIndices(
const std::vector<int>& indicesToRemove);
212 Teuchos::RCP< TraceMinRitzOp<Scalar,MV,OP> > Op_;
213 Teuchos::RCP< TraceMinProjOp<Scalar,MV,OP> > projector_;
215 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> > > solver_;
227 template <
class Scalar,
class MV,
class OP>
228 class TraceMinProjectedPrecOp :
public TraceMinOp<Scalar,MV,OP>
237 TraceMinProjectedPrecOp(
const Teuchos::RCP<const OP> Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs);
239 ~TraceMinProjectedPrecOp();
241 void Apply(
const MV& X, MV& Y)
const;
243 void removeIndices(
const std::vector<int>& indicesToRemove);
246 Teuchos::RCP<const OP> Op_;
247 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
249 Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
250 Teuchos::RCP<const OP> B_;
261 #ifdef HAVE_ANASAZI_BELOS 262 template <
class Scalar,
class MV,
class OP>
263 class TraceMinProjRitzOpWithPrec :
public TraceMinOp<Scalar,MV,OP>
271 TraceMinProjRitzOpWithPrec(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman = Teuchos::null);
272 TraceMinProjRitzOpWithPrec(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs);
274 void Apply(
const MV& X, MV& Y)
const;
278 void ApplyInverse(
const MV& X, MV& Y);
280 void removeIndices(
const std::vector<int>& indicesToRemove);
283 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Op_;
284 Teuchos::RCP<TraceMinProjectedPrecOp<Scalar,MV,OP> > Prec_;
286 Teuchos::RCP< Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> > > solver_;
287 Teuchos::RCP< Belos::LinearProblem<Scalar,MV,TraceMinOp<Scalar,MV,OP> > > problem_;
303 template <
class Scalar,
class MV,
class OP>
304 class OperatorTraits <Scalar, MV,
Experimental::TraceMinOp<Scalar,MV,OP> >
308 Apply (
const Experimental::TraceMinOp<Scalar,MV,OP>& Op,
310 MV& y) {Op.Apply(x,y);};
314 HasApplyTranspose (
const Experimental::TraceMinOp<Scalar,MV,OP>& Op) {
return true;};
319 #ifdef HAVE_ANASAZI_BELOS 322 template <
class Scalar,
class MV,
class OP>
323 class OperatorTraits <Scalar, MV,
Anasazi::Experimental::TraceMinOp<Scalar,MV,OP> >
327 Apply (
const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op,
329 MV& y, Belos::ETrans trans = NOTRANS) {Op.Apply(x,y);};
333 HasApplyTranspose (
const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op) {
return true;};
342 template <
class Scalar,
class MV,
class OP>
343 class OperatorTraits <Scalar, MV,
Experimental::TraceMinRitzOp<Scalar,MV,OP> >
347 Apply (
const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op,
349 MV& y) {Op.Apply(x,y);};
353 HasApplyTranspose (
const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op) {
return true;};
361 template <
class Scalar,
class MV,
class OP>
362 class OperatorTraits <Scalar, MV,
Experimental::TraceMinProjRitzOp<Scalar,MV,OP> >
366 Apply (
const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op,
368 MV& y) {Op.Apply(x,y);};
372 HasApplyTranspose (
const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op) {
return true;};
380 template <
class Scalar,
class MV,
class OP>
381 class OperatorTraits <Scalar, MV,
Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP> >
385 Apply (
const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op,
387 MV& y) {Op.Apply(x,y);};
391 HasApplyTranspose (
const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op) {
return true;};
404 template <
class Scalar,
class MV,
class OP>
406 ONE(Teuchos::ScalarTraits<Scalar>::one())
408 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 409 ProjTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinProjOp::Apply()");
414 if(orthman_ != Teuchos::null && B_ != Teuchos::null)
415 orthman_->setOp(Teuchos::null);
419 if(B_ != Teuchos::null)
423 if(orthman_ != Teuchos::null)
427 orthman_->normalizeMat(*helperMV);
428 projVecs_.push_back(helperMV);
432 std::vector<Scalar> normvec(nvec);
434 for(
int i=0; i<nvec; i++)
435 normvec[i] = ONE/normvec[i];
438 projVecs_.push_back(helperMV);
446 template <
class Scalar,
class MV,
class OP>
447 TraceMinProjOp<Scalar,MV,OP>::TraceMinProjOp(
const Teuchos::RCP<const MV> X,
const Teuchos::RCP<const OP> B, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs) :
448 ONE(Teuchos::ScalarTraits<Scalar>::one())
450 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 451 ProjTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinProjOp::Apply()");
456 if(B_ != Teuchos::null)
457 orthman_->setOp(Teuchos::null);
463 if(B_ != Teuchos::null)
467 orthman_->normalizeMat(*helperMV);
468 projVecs_.push_back(helperMV);
477 template <
class Scalar,
class MV,
class OP>
478 TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
480 if(orthman_ != Teuchos::null)
486 template <
class Scalar,
class MV,
class OP>
487 void TraceMinProjOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 489 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 490 Teuchos::TimeMonitor lcltimer( *ProjTime_ );
493 if(orthman_ != Teuchos::null)
496 orthman_->projectMat(Y,projVecs_);
500 int nvec = MVT::GetNumberVecs(X);
501 std::vector<Scalar> dotProducts(nvec);
502 MVT::MvDot(*projVecs_[0],X,dotProducts);
503 Teuchos::RCP<MV> helper = MVT::CloneCopy(*projVecs_[0]);
504 MVT::MvScale(*helper,dotProducts);
505 MVT::MvAddMv(ONE,X,-ONE,*helper,Y);
511 template <
class Scalar,
class MV,
class OP>
512 void TraceMinProjOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
514 if (orthman_ == Teuchos::null) {
515 const int nprojvecs = projVecs_.size();
516 const int nvecs = MVT::GetNumberVecs(*projVecs_[nprojvecs-1]);
517 const int numRemoving = indicesToRemove.size();
518 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
520 for (
int i=0; i<nvecs; i++) {
524 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
526 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*projVecs_[nprojvecs-1],indicesToLeave);
527 projVecs_[nprojvecs-1] = helperMV;
537 template <
class Scalar,
class MV,
class OP>
538 TraceMinRitzOp<Scalar,MV,OP>::TraceMinRitzOp(
const Teuchos::RCP<const OP>& A,
const Teuchos::RCP<const OP>& B,
const Teuchos::RCP<const OP>& Prec) :
539 ONE(Teuchos::ScalarTraits<Scalar>::one()),
540 ZERO(Teuchos::ScalarTraits<Scalar>::zero())
547 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 548 PetraMultTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinRitzOp: *Petra::Apply()");
549 AopMultTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinRitzOp::Apply()");
553 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
557 if(Prec != Teuchos::null)
558 Prec_ = Teuchos::rcp(
new TraceMinRitzOp<Scalar,MV,OP>(Prec) );
561 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> >(linSolOp,Prec_) );
566 template <
class Scalar,
class MV,
class OP>
567 void TraceMinRitzOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 569 int nvecs = MVT::GetNumberVecs(X);
571 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 572 Teuchos::TimeMonitor outertimer( *AopMultTime_ );
577 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 578 Teuchos::TimeMonitor lcltimer( *PetraMultTime_ );
585 if(ritzShifts_.size() > 0)
588 std::vector<int> nonzeroRitzIndices;
589 nonzeroRitzIndices.reserve(nvecs);
590 for(
int i=0; i<nvecs; i++)
592 if(ritzShifts_[i] != ZERO)
593 nonzeroRitzIndices.push_back(i);
597 int numRitzShifts = nonzeroRitzIndices.size();
598 if(numRitzShifts > 0)
601 Teuchos::RCP<const MV> ritzX = MVT::CloneView(X,nonzeroRitzIndices);
602 Teuchos::RCP<MV> ritzY = MVT::CloneViewNonConst(Y,nonzeroRitzIndices);
605 std::vector<Scalar> nonzeroRitzShifts(numRitzShifts);
606 for(
int i=0; i<numRitzShifts; i++)
607 nonzeroRitzShifts[i] = ritzShifts_[nonzeroRitzIndices[i]];
610 if(B_ != Teuchos::null)
612 Teuchos::RCP<MV> BX = MVT::Clone(*ritzX,numRitzShifts);
613 OPT::Apply(*B_,*ritzX,*BX);
614 MVT::MvScale(*BX,nonzeroRitzShifts);
615 MVT::MvAddMv(ONE,*ritzY,-ONE,*BX,*ritzY);
620 Teuchos::RCP<MV> scaledX = MVT::CloneCopy(*ritzX);
621 MVT::MvScale(*scaledX,nonzeroRitzShifts);
622 MVT::MvAddMv(ONE,*ritzY,-ONE,*scaledX,*ritzY);
630 template <
class Scalar,
class MV,
class OP>
631 void TraceMinRitzOp<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
633 int nvecs = MVT::GetNumberVecs(X);
634 std::vector<int> indices(nvecs);
635 for(
int i=0; i<nvecs; i++)
638 Teuchos::RCP<const MV> rcpX = MVT::CloneView(X,indices);
639 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
642 solver_->setTol(tolerances_);
643 solver_->setMaxIter(maxits_);
646 solver_->setSol(rcpY);
647 solver_->setRHS(rcpX);
655 template <
class Scalar,
class MV,
class OP>
656 void TraceMinRitzOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
658 int nvecs = tolerances_.size();
659 int numRemoving = indicesToRemove.size();
660 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
661 std::vector<Scalar> helperS(nvecs-numRemoving);
663 for(
int i=0; i<nvecs; i++)
666 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
668 for(
int i=0; i<nvecs-numRemoving; i++)
669 helperS[i] = ritzShifts_[indicesToLeave[i]];
670 ritzShifts_ = helperS;
672 for(
int i=0; i<nvecs-numRemoving; i++)
673 helperS[i] = tolerances_[indicesToLeave[i]];
674 tolerances_ = helperS;
684 template <
class Scalar,
class MV,
class OP>
685 TraceMinProjRitzOp<Scalar,MV,OP>::TraceMinProjRitzOp(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs,
const Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman)
690 projector_ = Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman) );
693 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
697 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
701 template <
class Scalar,
class MV,
class OP>
702 TraceMinProjRitzOp<Scalar,MV,OP>::TraceMinProjRitzOp(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs,
const Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array<Teuchos::RCP<const MV> > auxVecs)
707 projector_ = Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman, auxVecs) );
710 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
714 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
720 template <
class Scalar,
class MV,
class OP>
721 void TraceMinProjRitzOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 723 int nvecs = MVT::GetNumberVecs(X);
730 Teuchos::RCP<MV> APX = MVT::Clone(X,nvecs);
731 OPT::Apply(*Op_,X,*APX);
734 projector_->Apply(*APX,Y);
739 template <
class Scalar,
class MV,
class OP>
740 void TraceMinProjRitzOp<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
742 int nvecs = MVT::GetNumberVecs(X);
743 std::vector<int> indices(nvecs);
744 for(
int i=0; i<nvecs; i++)
747 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
748 Teuchos::RCP<MV> PX = MVT::Clone(X,nvecs);
749 projector_->Apply(X,*PX);
752 solver_->setTol(Op_->tolerances_);
753 solver_->setMaxIter(Op_->maxits_);
756 solver_->setSol(rcpY);
765 template <
class Scalar,
class MV,
class OP>
766 void TraceMinProjRitzOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
768 Op_->removeIndices(indicesToRemove);
770 projector_->removeIndices(indicesToRemove);
776 #ifdef HAVE_ANASAZI_BELOS 782 template <
class Scalar,
class MV,
class OP>
783 TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::TraceMinProjRitzOpWithPrec(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
784 ONE(Teuchos::ScalarTraits<Scalar>::one())
789 Teuchos::RCP<TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
793 problem_ = Teuchos::rcp(
new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
796 problem_->setOperator(linSolOp);
800 Prec_ = Teuchos::rcp(
new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman) );
805 solver_ = Teuchos::rcp(
new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
809 template <
class Scalar,
class MV,
class OP>
810 TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::TraceMinProjRitzOpWithPrec(
const Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> >& Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs) :
811 ONE(Teuchos::ScalarTraits<Scalar>::one())
816 Teuchos::RCP<const TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
820 problem_ = Teuchos::rcp(
new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
823 problem_->setOperator(linSolOp);
827 Prec_ = Teuchos::rcp(
new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman, auxVecs) );
832 solver_ = Teuchos::rcp(
new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
836 template <
class Scalar,
class MV,
class OP>
837 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 839 int nvecs = MVT::GetNumberVecs(X);
840 RCP<MV> Ydot = MVT::Clone(Y,nvecs);
841 OPT::Apply(*Op_,X,*Ydot);
842 Prec_->Apply(*Ydot,Y);
846 template <
class Scalar,
class MV,
class OP>
847 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
849 int nvecs = MVT::GetNumberVecs(X);
850 std::vector<int> indices(nvecs);
851 for(
int i=0; i<nvecs; i++)
854 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
855 Teuchos::RCP<MV> rcpX = MVT::Clone(X,nvecs);
857 Prec_->Apply(X,*rcpX);
860 problem_->setProblem(rcpY,rcpX);
863 solver_->setProblem(problem_);
869 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(
new Teuchos::ParameterList());
870 pl->set(
"Convergence Tolerance", Op_->tolerances_[0]);
871 pl->set(
"Block Size", nvecs);
874 pl->set(
"Maximum Iterations", Op_->getMaxIts());
875 pl->set(
"Num Blocks", Op_->getMaxIts());
876 solver_->setParameters(pl);
883 template <
class Scalar,
class MV,
class OP>
884 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
886 Op_->removeIndices(indicesToRemove);
888 Prec_->removeIndices(indicesToRemove);
898 template <
class Scalar,
class MV,
class OP>
899 TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(
const Teuchos::RCP<const OP> Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman) :
900 ONE (Teuchos::ScalarTraits<Scalar>::one())
906 Teuchos::RCP<MV> helperMV =
MVT::Clone(*projVecs,nvecs);
911 if(orthman_ != Teuchos::null)
914 B_ = orthman_->getOp();
915 orthman_->setOp(Op_);
920 const int rank = orthman_->normalizeMat (*locProjVecs, Teuchos::null, helperMV);
923 TEUCHOS_TEST_FOR_EXCEPTION(rank != nvecs, std::runtime_error,
"Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
925 orthman_->setOp(Teuchos::null);
929 std::vector<Scalar> dotprods(nvecs);
932 for(
int i=0; i<nvecs; i++)
933 dotprods[i] = ONE/sqrt(dotprods[i]);
938 projVecs_.push_back(helperMV);
942 template <
class Scalar,
class MV,
class OP>
943 TraceMinProjectedPrecOp<Scalar,MV,OP>::TraceMinProjectedPrecOp(
const Teuchos::RCP<const OP> Op,
const Teuchos::RCP<const MV> projVecs, Teuchos::RCP<
Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman, Teuchos::Array< Teuchos::RCP<const MV> > auxVecs) :
944 ONE(Teuchos::ScalarTraits<Scalar>::one())
950 Teuchos::RCP<MV> locProjVecs;
953 if(auxVecs.size() > 0)
957 for(
int i=0; i<auxVecs.size(); i++)
965 std::vector<int> locind(nvecs);
968 for (
size_t i = 0; i<locind.size(); i++) {
969 locind[i] = startIndex + i;
971 startIndex += locind.size();
974 for (
size_t i=0; i < static_cast<size_t> (auxVecs.size ()); ++i)
977 for(
size_t j=0; j<locind.size(); j++) locind[j] = startIndex + j;
978 startIndex += locind.size();
989 Teuchos::RCP<MV> helperMV =
MVT::Clone(*projVecs,nvecs);
995 B_ = orthman_->getOp();
996 orthman_->setOp(Op_);
999 const int rank = orthman_->normalizeMat(*locProjVecs,Teuchos::null,helperMV);
1001 projVecs_.push_back(helperMV);
1006 TEUCHOS_TEST_FOR_EXCEPTION(
1007 rank != nvecs, std::runtime_error,
1008 "Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
1010 orthman_->setOp(Teuchos::null);
1014 template <
class Scalar,
class MV,
class OP>
1015 TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
1017 if(orthman_ != Teuchos::null)
1018 orthman_->setOp(B_);
1023 template <
class Scalar,
class MV,
class OP>
1024 void TraceMinProjectedPrecOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 1026 int nvecsX = MVT::GetNumberVecs(X);
1028 if(orthman_ != Teuchos::null)
1031 int nvecsP = MVT::GetNumberVecs(*projVecs_[0]);
1032 OPT::Apply(*Op_,X,Y);
1034 Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> > projX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,Scalar>(nvecsP,nvecsX));
1036 MVT::MvTransMv(ONE, *projVecs_[0], X, *projX);
1038 MVT::MvTimesMatAddMv(-ONE, *projVecs_[0], *projX, ONE, Y);
1042 Teuchos::RCP<MV> MX = MVT::Clone(X,nvecsX);
1043 OPT::Apply(*Op_,X,*MX);
1045 std::vector<Scalar> dotprods(nvecsX);
1046 MVT::MvDot(*projVecs_[0], X, dotprods);
1048 Teuchos::RCP<MV> helper = MVT::CloneCopy(*projVecs_[0]);
1049 MVT::MvScale(*helper, dotprods);
1050 MVT::MvAddMv(ONE, *MX, -ONE, *helper, Y);
1055 template <
class Scalar,
class MV,
class OP>
1056 void TraceMinProjectedPrecOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
1058 if(orthman_ == Teuchos::null)
1060 int nvecs = MVT::GetNumberVecs(*projVecs_[0]);
1061 int numRemoving = indicesToRemove.size();
1062 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
1064 for(
int i=0; i<nvecs; i++)
1067 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
1069 Teuchos::RCP<const MV> helperMV = MVT::CloneCopy(*projVecs_[0],indicesToLeave);
1070 projVecs_[0] = helperMV;
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.