46 #ifndef TRACEMIN_RITZ_OP_HPP 47 #define TRACEMIN_RITZ_OP_HPP 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" 62 #ifdef HAVE_ANASAZI_EPETRA 63 #include "BelosEpetraAdapter.hpp" 67 #include "Teuchos_RCP.hpp" 68 #include "Teuchos_SerialDenseSolver.hpp" 69 #include "Teuchos_ScalarTraitsDecl.hpp" 84 template <
class Scalar,
class MV,
class OP>
88 virtual ~TraceMinOp() { };
89 virtual void Apply(
const MV& X, MV& Y)
const =0;
90 virtual void removeIndices(
const std::vector<int>& indicesToRemove) =0;
101 template <
class Scalar,
class MV,
class OP>
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);
116 void Apply(
const MV& X, MV& Y)
const;
119 void removeIndices(
const std::vector<int>& indicesToRemove);
122 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
123 Teuchos::RCP<Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
124 Teuchos::RCP<const OP> B_;
126 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 127 RCP<Teuchos::Time> ProjTime_;
138 template <
class Scalar,
class MV,
class OP>
139 class TraceMinRitzOp :
public TraceMinOp<Scalar,MV,OP>
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;
152 TraceMinRitzOp(
const Teuchos::RCP<const OP>& A,
const Teuchos::RCP<const OP>& B = Teuchos::null,
const Teuchos::RCP<const OP>& Prec = Teuchos::null);
155 ~TraceMinRitzOp() { };
158 void setRitzShifts(std::vector<Scalar> shifts) {ritzShifts_ = shifts;};
160 Scalar getRitzShift(
const int subscript) {
return ritzShifts_[subscript]; };
162 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > getPrec() {
return Prec_; };
165 void setInnerTol(
const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
167 void getInnerTol(std::vector<Scalar>& tolerances)
const { tolerances = tolerances_; };
169 void setMaxIts(
const int maxits) { maxits_ = maxits; };
171 int getMaxIts()
const {
return maxits_; };
174 void Apply(
const MV& X, MV& Y)
const;
177 void ApplyInverse(
const MV& X, MV& Y);
179 void removeIndices(
const std::vector<int>& indicesToRemove);
182 Teuchos::RCP<const OP> A_, B_;
183 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Prec_;
186 std::vector<Scalar> ritzShifts_;
187 std::vector<Scalar> tolerances_;
189 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> > > solver_;
191 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 192 RCP<Teuchos::Time> PetraMultTime_, AopMultTime_;
204 template <
class Scalar,
class MV,
class OP>
205 class TraceMinProjRitzOp :
public TraceMinOp<Scalar,MV,OP>
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);
216 void Apply(
const MV& X, MV& Y)
const;
220 void ApplyInverse(
const MV& X, MV& Y);
222 void removeIndices(
const std::vector<int>& indicesToRemove);
225 Teuchos::RCP< TraceMinRitzOp<Scalar,MV,OP> > Op_;
226 Teuchos::RCP< TraceMinProjOp<Scalar,MV,OP> > projector_;
228 Teuchos::RCP< PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> > > solver_;
240 template <
class Scalar,
class MV,
class OP>
241 class TraceMinProjectedPrecOp :
public TraceMinOp<Scalar,MV,OP>
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);
252 ~TraceMinProjectedPrecOp();
254 void Apply(
const MV& X, MV& Y)
const;
256 void removeIndices(
const std::vector<int>& indicesToRemove);
259 Teuchos::RCP<const OP> Op_;
260 Teuchos::Array< Teuchos::RCP<const MV> > projVecs_;
262 Teuchos::RCP< Anasazi::MatOrthoManager<Scalar,MV,OP> > orthman_;
263 Teuchos::RCP<const OP> B_;
274 #ifdef HAVE_ANASAZI_BELOS 275 template <
class Scalar,
class MV,
class OP>
276 class TraceMinProjRitzOpWithPrec :
public TraceMinOp<Scalar,MV,OP>
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);
287 void Apply(
const MV& X, MV& Y)
const;
291 void ApplyInverse(
const MV& X, MV& Y);
293 void removeIndices(
const std::vector<int>& indicesToRemove);
296 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > Op_;
297 Teuchos::RCP<TraceMinProjectedPrecOp<Scalar,MV,OP> > Prec_;
299 Teuchos::RCP< Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> > > solver_;
300 Teuchos::RCP< Belos::LinearProblem<Scalar,MV,TraceMinOp<Scalar,MV,OP> > > problem_;
316 template <
class Scalar,
class MV,
class OP>
317 class OperatorTraits <Scalar, MV,
Experimental::TraceMinOp<Scalar,MV,OP> >
321 Apply (
const Experimental::TraceMinOp<Scalar,MV,OP>& Op,
323 MV& y) {Op.Apply(x,y);};
327 HasApplyTranspose (
const Experimental::TraceMinOp<Scalar,MV,OP>& Op) {
return true;};
332 #ifdef HAVE_ANASAZI_BELOS 335 template <
class Scalar,
class MV,
class OP>
336 class OperatorTraits <Scalar, MV,
Anasazi::Experimental::TraceMinOp<Scalar,MV,OP> >
340 Apply (
const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op,
342 MV& y, Belos::ETrans trans = NOTRANS) {Op.Apply(x,y);};
346 HasApplyTranspose (
const Anasazi::Experimental::TraceMinOp<Scalar,MV,OP>& Op) {
return true;};
355 template <
class Scalar,
class MV,
class OP>
356 class OperatorTraits <Scalar, MV,
Experimental::TraceMinRitzOp<Scalar,MV,OP> >
360 Apply (
const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op,
362 MV& y) {Op.Apply(x,y);};
366 HasApplyTranspose (
const Experimental::TraceMinRitzOp<Scalar,MV,OP>& Op) {
return true;};
374 template <
class Scalar,
class MV,
class OP>
375 class OperatorTraits <Scalar, MV,
Experimental::TraceMinProjRitzOp<Scalar,MV,OP> >
379 Apply (
const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op,
381 MV& y) {Op.Apply(x,y);};
385 HasApplyTranspose (
const Experimental::TraceMinProjRitzOp<Scalar,MV,OP>& Op) {
return true;};
393 template <
class Scalar,
class MV,
class OP>
394 class OperatorTraits <Scalar, MV,
Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP> >
398 Apply (
const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op,
400 MV& y) {Op.Apply(x,y);};
404 HasApplyTranspose (
const Experimental::TraceMinProjectedPrecOp<Scalar,MV,OP>& Op) {
return true;};
417 template <
class Scalar,
class MV,
class OP>
419 ONE(Teuchos::ScalarTraits<Scalar>::one())
421 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 422 ProjTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinProjOp::Apply()");
427 if(orthman_ != Teuchos::null && B_ != Teuchos::null)
428 orthman_->setOp(Teuchos::null);
432 if(B_ != Teuchos::null)
436 if(orthman_ != Teuchos::null)
440 orthman_->normalizeMat(*helperMV);
441 projVecs_.push_back(helperMV);
445 std::vector<Scalar> normvec(nvec);
447 for(
int i=0; i<nvec; i++)
448 normvec[i] = ONE/normvec[i];
451 projVecs_.push_back(helperMV);
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())
463 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 464 ProjTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinProjOp::Apply()");
469 if(B_ != Teuchos::null)
470 orthman_->setOp(Teuchos::null);
476 if(B_ != Teuchos::null)
480 orthman_->normalizeMat(*helperMV);
481 projVecs_.push_back(helperMV);
490 template <
class Scalar,
class MV,
class OP>
491 TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
493 if(orthman_ != Teuchos::null)
499 template <
class Scalar,
class MV,
class OP>
500 void TraceMinProjOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 502 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 503 Teuchos::TimeMonitor lcltimer( *ProjTime_ );
506 if(orthman_ != Teuchos::null)
509 orthman_->projectMat(Y,projVecs_);
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);
524 template <
class Scalar,
class MV,
class OP>
525 void TraceMinProjOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
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);
533 for (
int i=0; i<nvecs; i++) {
537 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
539 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*projVecs_[nprojvecs-1],indicesToLeave);
540 projVecs_[nprojvecs-1] = helperMV;
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())
560 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 561 PetraMultTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinRitzOp: *Petra::Apply()");
562 AopMultTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinRitzOp::Apply()");
566 Teuchos::RCP<TraceMinRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
570 if(Prec != Teuchos::null)
571 Prec_ = Teuchos::rcp(
new TraceMinRitzOp<Scalar,MV,OP>(Prec) );
574 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinRitzOp<Scalar,MV,OP> >(linSolOp,Prec_) );
579 template <
class Scalar,
class MV,
class OP>
580 void TraceMinRitzOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 582 int nvecs = MVT::GetNumberVecs(X);
584 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 585 Teuchos::TimeMonitor outertimer( *AopMultTime_ );
590 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 591 Teuchos::TimeMonitor lcltimer( *PetraMultTime_ );
598 if(ritzShifts_.size() > 0)
601 std::vector<int> nonzeroRitzIndices;
602 nonzeroRitzIndices.reserve(nvecs);
603 for(
int i=0; i<nvecs; i++)
605 if(ritzShifts_[i] != ZERO)
606 nonzeroRitzIndices.push_back(i);
610 int numRitzShifts = nonzeroRitzIndices.size();
611 if(numRitzShifts > 0)
614 Teuchos::RCP<const MV> ritzX = MVT::CloneView(X,nonzeroRitzIndices);
615 Teuchos::RCP<MV> ritzY = MVT::CloneViewNonConst(Y,nonzeroRitzIndices);
618 std::vector<Scalar> nonzeroRitzShifts(numRitzShifts);
619 for(
int i=0; i<numRitzShifts; i++)
620 nonzeroRitzShifts[i] = ritzShifts_[nonzeroRitzIndices[i]];
623 if(B_ != Teuchos::null)
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);
633 Teuchos::RCP<MV> scaledX = MVT::CloneCopy(*ritzX);
634 MVT::MvScale(*scaledX,nonzeroRitzShifts);
635 MVT::MvAddMv(ONE,*ritzY,-ONE,*scaledX,*ritzY);
643 template <
class Scalar,
class MV,
class OP>
644 void TraceMinRitzOp<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
646 int nvecs = MVT::GetNumberVecs(X);
647 std::vector<int> indices(nvecs);
648 for(
int i=0; i<nvecs; i++)
651 Teuchos::RCP<const MV> rcpX = MVT::CloneView(X,indices);
652 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
655 solver_->setTol(tolerances_);
656 solver_->setMaxIter(maxits_);
659 solver_->setSol(rcpY);
660 solver_->setRHS(rcpX);
668 template <
class Scalar,
class MV,
class OP>
669 void TraceMinRitzOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
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);
676 for(
int i=0; i<nvecs; i++)
679 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
681 for(
int i=0; i<nvecs-numRemoving; i++)
682 helperS[i] = ritzShifts_[indicesToLeave[i]];
683 ritzShifts_ = helperS;
685 for(
int i=0; i<nvecs-numRemoving; i++)
686 helperS[i] = tolerances_[indicesToLeave[i]];
687 tolerances_ = helperS;
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)
703 projector_ = Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman) );
706 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
710 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
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)
720 projector_ = Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman, auxVecs) );
723 Teuchos::RCP<TraceMinProjRitzOp<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
727 solver_ = Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
733 template <
class Scalar,
class MV,
class OP>
734 void TraceMinProjRitzOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 736 int nvecs = MVT::GetNumberVecs(X);
743 Teuchos::RCP<MV> APX = MVT::Clone(X,nvecs);
744 OPT::Apply(*Op_,X,*APX);
747 projector_->Apply(*APX,Y);
752 template <
class Scalar,
class MV,
class OP>
753 void TraceMinProjRitzOp<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
755 int nvecs = MVT::GetNumberVecs(X);
756 std::vector<int> indices(nvecs);
757 for(
int i=0; i<nvecs; i++)
760 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
761 Teuchos::RCP<MV> PX = MVT::Clone(X,nvecs);
762 projector_->Apply(X,*PX);
765 solver_->setTol(Op_->tolerances_);
766 solver_->setMaxIter(Op_->maxits_);
769 solver_->setSol(rcpY);
778 template <
class Scalar,
class MV,
class OP>
779 void TraceMinProjRitzOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
781 Op_->removeIndices(indicesToRemove);
783 projector_->removeIndices(indicesToRemove);
789 #ifdef HAVE_ANASAZI_BELOS 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())
802 Teuchos::RCP<TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
806 problem_ = Teuchos::rcp(
new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
809 problem_->setOperator(linSolOp);
813 Prec_ = Teuchos::rcp(
new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman) );
818 solver_ = Teuchos::rcp(
new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
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())
829 Teuchos::RCP<const TraceMinProjRitzOpWithPrec<Scalar,MV,OP> > linSolOp = Teuchos::rcp(
this);
833 problem_ = Teuchos::rcp(
new Belos::LinearProblem< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
836 problem_->setOperator(linSolOp);
840 Prec_ = Teuchos::rcp(
new TraceMinProjectedPrecOp<Scalar,MV,OP>(Op_->Prec_->A_, projVecs, orthman, auxVecs) );
845 solver_ = Teuchos::rcp(
new Belos::PseudoBlockGmresSolMgr< Scalar,MV,TraceMinOp<Scalar,MV,OP> >());
849 template <
class Scalar,
class MV,
class OP>
850 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 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);
859 template <
class Scalar,
class MV,
class OP>
860 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::ApplyInverse(
const MV& X, MV& Y)
862 int nvecs = MVT::GetNumberVecs(X);
863 std::vector<int> indices(nvecs);
864 for(
int i=0; i<nvecs; i++)
867 Teuchos::RCP<MV> rcpY = MVT::CloneViewNonConst(Y,indices);
868 Teuchos::RCP<MV> rcpX = MVT::Clone(X,nvecs);
870 Prec_->Apply(X,*rcpX);
873 problem_->setProblem(rcpY,rcpX);
876 solver_->setProblem(problem_);
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);
887 pl->set(
"Maximum Iterations", Op_->getMaxIts());
888 pl->set(
"Num Blocks", Op_->getMaxIts());
889 solver_->setParameters(pl);
896 template <
class Scalar,
class MV,
class OP>
897 void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
899 Op_->removeIndices(indicesToRemove);
901 Prec_->removeIndices(indicesToRemove);
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())
919 Teuchos::RCP<MV> helperMV =
MVT::Clone(*projVecs,nvecs);
924 if(orthman_ != Teuchos::null)
927 B_ = orthman_->getOp();
928 orthman_->setOp(Op_);
933 const int rank = orthman_->normalizeMat (*locProjVecs, Teuchos::null, helperMV);
936 TEUCHOS_TEST_FOR_EXCEPTION(rank != nvecs, std::runtime_error,
"Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
938 orthman_->setOp(Teuchos::null);
942 std::vector<Scalar> dotprods(nvecs);
945 for(
int i=0; i<nvecs; i++)
946 dotprods[i] = ONE/sqrt(dotprods[i]);
951 projVecs_.push_back(helperMV);
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())
963 Teuchos::RCP<MV> locProjVecs;
966 if(auxVecs.size() > 0)
970 for(
int i=0; i<auxVecs.size(); i++)
978 std::vector<int> locind(nvecs);
981 for (
size_t i = 0; i<locind.size(); i++) {
982 locind[i] = startIndex + i;
984 startIndex += locind.size();
987 for (
size_t i=0; i < static_cast<size_t> (auxVecs.size ()); ++i)
990 for(
size_t j=0; j<locind.size(); j++) locind[j] = startIndex + j;
991 startIndex += locind.size();
1002 Teuchos::RCP<MV> helperMV =
MVT::Clone(*projVecs,nvecs);
1008 B_ = orthman_->getOp();
1009 orthman_->setOp(Op_);
1012 const int rank = orthman_->normalizeMat(*locProjVecs,Teuchos::null,helperMV);
1014 projVecs_.push_back(helperMV);
1019 TEUCHOS_TEST_FOR_EXCEPTION(
1020 rank != nvecs, std::runtime_error,
1021 "Belos::TraceMinProjectedPrecOp: constructor: Loss of orthogonality detected.");
1023 orthman_->setOp(Teuchos::null);
1027 template <
class Scalar,
class MV,
class OP>
1028 TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
1030 if(orthman_ != Teuchos::null)
1031 orthman_->setOp(B_);
1036 template <
class Scalar,
class MV,
class OP>
1037 void TraceMinProjectedPrecOp<Scalar,MV,OP>::Apply(
const MV& X, MV& Y)
const 1039 int nvecsX = MVT::GetNumberVecs(X);
1041 if(orthman_ != Teuchos::null)
1044 int nvecsP = MVT::GetNumberVecs(*projVecs_[0]);
1045 OPT::Apply(*Op_,X,Y);
1047 Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> > projX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,Scalar>(nvecsP,nvecsX));
1049 MVT::MvTransMv(ONE, *projVecs_[0], X, *projX);
1051 MVT::MvTimesMatAddMv(-ONE, *projVecs_[0], *projX, ONE, Y);
1055 Teuchos::RCP<MV> MX = MVT::Clone(X,nvecsX);
1056 OPT::Apply(*Op_,X,*MX);
1058 std::vector<Scalar> dotprods(nvecsX);
1059 MVT::MvDot(*projVecs_[0], X, dotprods);
1061 Teuchos::RCP<MV> helper = MVT::CloneCopy(*projVecs_[0]);
1062 MVT::MvScale(*helper, dotprods);
1063 MVT::MvAddMv(ONE, *MX, -ONE, *helper, Y);
1068 template <
class Scalar,
class MV,
class OP>
1069 void TraceMinProjectedPrecOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
1071 if(orthman_ == Teuchos::null)
1073 int nvecs = MVT::GetNumberVecs(*projVecs_[0]);
1074 int numRemoving = indicesToRemove.size();
1075 std::vector<int> helper(nvecs), indicesToLeave(nvecs-numRemoving);
1077 for(
int i=0; i<nvecs; i++)
1080 std::set_difference(helper.begin(), helper.end(), indicesToRemove.begin(), indicesToRemove.end(), indicesToLeave.begin());
1082 Teuchos::RCP<const MV> helperMV = MVT::CloneCopy(*projVecs_[0],indicesToLeave);
1083 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.