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"
84template <
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;
101template <
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_;
138template <
class Scalar,
class MV,
class OP>
139class 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_;
204template <
class Scalar,
class MV,
class OP>
205class 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_;
240template <
class Scalar,
class MV,
class OP>
241class 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
275template <
class Scalar,
class MV,
class OP>
276class 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_;
316template <
class Scalar,
class MV,
class OP>
317class 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
335template <
class Scalar,
class MV,
class OP>
336class 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;};
355template <
class Scalar,
class MV,
class OP>
356class 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;};
374template <
class Scalar,
class MV,
class OP>
375class 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;};
393template <
class Scalar,
class MV,
class OP>
394class 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;};
417template <
class Scalar,
class MV,
class OP>
419ONE(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)
434 int nvec = MVT::GetNumberVecs(*X);
436 if(orthman_ != Teuchos::null)
439 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*X);
440 orthman_->normalizeMat(*helperMV);
441 projVecs_.push_back(helperMV);
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);
455 projVecs_.push_back(MVT::CloneCopy(*X));
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())
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)
479 Teuchos::RCP<MV> helperMV = MVT::CloneCopy(*X);
480 orthman_->normalizeMat(*helperMV);
481 projVecs_.push_back(helperMV);
485 projVecs_.push_back(MVT::CloneCopy(*X));
490template <
class Scalar,
class MV,
class OP>
491TraceMinProjOp<Scalar,MV,OP>::~TraceMinProjOp()
493 if(orthman_ != Teuchos::null)
499template <
class Scalar,
class MV,
class OP>
500void 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);
524template <
class Scalar,
class MV,
class OP>
525void 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;
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())
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_) );
579template <
class Scalar,
class MV,
class OP>
580void 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);
643template <
class Scalar,
class MV,
class OP>
644void 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);
668template <
class Scalar,
class MV,
class OP>
669void 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;
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)
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) );
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)
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) );
733template <
class Scalar,
class MV,
class OP>
734void 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);
752template <
class Scalar,
class MV,
class OP>
753void 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);
778template <
class Scalar,
class MV,
class OP>
779void TraceMinProjRitzOp<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
781 Op_->removeIndices(indicesToRemove);
783 projector_->removeIndices(indicesToRemove);
789#ifdef HAVE_ANASAZI_BELOS
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())
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> >());
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())
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> >());
849template <
class Scalar,
class MV,
class OP>
850void 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);
859template <
class Scalar,
class MV,
class OP>
860void 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);
896template <
class Scalar,
class MV,
class OP>
897void TraceMinProjRitzOpWithPrec<Scalar,MV,OP>::removeIndices(
const std::vector<int>& indicesToRemove)
899 Op_->removeIndices(indicesToRemove);
901 Prec_->removeIndices(indicesToRemove);
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())
918 int nvecs = MVT::GetNumberVecs(*projVecs);
919 Teuchos::RCP<MV> helperMV = MVT::Clone(*projVecs,nvecs);
922 OPT::Apply(*Op_,*projVecs,*helperMV);
924 if(orthman_ != Teuchos::null)
927 B_ = orthman_->getOp();
928 orthman_->setOp(Op_);
930 Teuchos::RCP<MV> locProjVecs = MVT::CloneCopy(*projVecs);
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);
943 MVT::MvDot(*projVecs,*helperMV,dotprods);
945 for(
int i=0; i<nvecs; i++)
946 dotprods[i] = ONE/sqrt(dotprods[i]);
948 MVT::MvScale(*helperMV, dotprods);
951 projVecs_.push_back(helperMV);
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())
963 Teuchos::RCP<MV> locProjVecs;
966 if(auxVecs.size() > 0)
969 nvecs = MVT::GetNumberVecs(*projVecs);
970 for(
int i=0; i<auxVecs.size(); i++)
971 nvecs += MVT::GetNumberVecs(*auxVecs[i]);
974 locProjVecs = MVT::Clone(*projVecs, nvecs);
978 std::vector<int> locind(nvecs);
980 locind.resize(MVT::GetNumberVecs(*projVecs));
981 for (
size_t i = 0; i<locind.size(); i++) {
982 locind[i] = startIndex + i;
984 startIndex += locind.size();
985 MVT::SetBlock(*projVecs,locind,*locProjVecs);
987 for (
size_t i=0; i < static_cast<size_t> (auxVecs.size ()); ++i)
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);
998 nvecs = MVT::GetNumberVecs(*projVecs);
999 locProjVecs = MVT::CloneCopy(*projVecs);
1002 Teuchos::RCP<MV> helperMV = MVT::Clone(*projVecs,nvecs);
1005 OPT::Apply(*Op_,*locProjVecs,*helperMV);
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);
1027template <
class Scalar,
class MV,
class OP>
1028TraceMinProjectedPrecOp<Scalar,MV,OP>::~TraceMinProjectedPrecOp()
1030 if(orthman_ != Teuchos::null)
1031 orthman_->setOp(B_);
1036template <
class Scalar,
class MV,
class OP>
1037void 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);
1068template <
class Scalar,
class MV,
class OP>
1069void 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;
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.