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_);
501 std::vector<Scalar> dotProducts(nvec);
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();
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 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);
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);
630 template <
class Scalar,
class MV,
class OP>
631 void TraceMinRitzOp<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
634 std::vector<int> indices(nvecs);
635 for(
int i=0; i<nvecs; i++)
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 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)
743 std::vector<int> indices(nvecs);
744 for(
int i=0; i<nvecs; i++)
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 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)
850 std::vector<int> indices(nvecs);
851 for(
int i=0; i<nvecs; i++)
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 1028 if(orthman_ != Teuchos::null)
1034 Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> > projX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,Scalar>(nvecsP,nvecsX));
1045 std::vector<Scalar> dotprods(nvecsX);
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)
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;
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Abstract base class for trace minimization eigensolvers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Virtual base class which defines basic traits for the operator type.
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
static void Assign(const MV &A, MV &mv)
mv := A
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.