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" 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>
103 void Apply(
const MV& X, MV& Y)
const;
106 void removeIndices(
const std::vector<int>& indicesToRemove);
113 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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;
142 ~TraceMinRitzOp() { };
145 void setRitzShifts(std::vector<Scalar> shifts) {ritzShifts_ = shifts;};
147 Scalar getRitzShift(
const int subscript) {
return ritzShifts_[subscript]; };
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);
173 std::vector<Scalar> ritzShifts_;
174 std::vector<Scalar> tolerances_;
178 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 191 template <
class Scalar,
class MV,
class OP>
192 class TraceMinProjRitzOp :
public TraceMinOp<Scalar,MV,OP>
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);
227 template <
class Scalar,
class MV,
class OP>
228 class TraceMinProjectedPrecOp :
public TraceMinOp<Scalar,MV,OP>
239 ~TraceMinProjectedPrecOp();
241 void Apply(
const MV& X, MV& Y)
const;
243 void removeIndices(
const std::vector<int>& indicesToRemove);
261 #ifdef HAVE_ANASAZI_BELOS 262 template <
class Scalar,
class MV,
class OP>
263 class TraceMinProjRitzOpWithPrec :
public TraceMinOp<Scalar,MV,OP>
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);
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 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>
448 ONE(
Teuchos::ScalarTraits<Scalar>::one())
450 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 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);
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>
539 ONE(
Teuchos::ScalarTraits<Scalar>::one()),
540 ZERO(
Teuchos::ScalarTraits<Scalar>::zero())
547 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 577 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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)
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)
613 OPT::Apply(*B_,*ritzX,*BX);
614 MVT::MvScale(*BX,nonzeroRitzShifts);
615 MVT::MvAddMv(ONE,*ritzY,-ONE,*BX,*ritzY);
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++)
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>
690 projector_ =
Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman) );
697 solver_ =
Teuchos::rcp(
new PseudoBlockMinres< Scalar,MV,TraceMinProjRitzOp<Scalar,MV,OP> >(linSolOp) );
701 template <
class Scalar,
class MV,
class OP>
707 projector_ =
Teuchos::rcp(
new TraceMinProjOp<Scalar,MV,OP>(projVecs, Op_->B_, orthman, auxVecs) );
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);
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++)
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>
784 ONE(
Teuchos::ScalarTraits<Scalar>::one())
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>
811 ONE(
Teuchos::ScalarTraits<Scalar>::one())
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++)
857 Prec_->Apply(X,*rcpX);
860 problem_->setProblem(rcpY,rcpX);
863 solver_->setProblem(problem_);
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());
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>
900 ONE (
Teuchos::ScalarTraits<Scalar>::one())
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>
944 ONE(
Teuchos::ScalarTraits<Scalar>::one())
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();
995 B_ = orthman_->getOp();
996 orthman_->setOp(Op_);
999 const int rank = orthman_->normalizeMat(*locProjVecs,Teuchos::null,helperMV);
1001 projVecs_.push_back(helperMV);
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);
1036 MVT::MvTransMv(ONE, *projVecs_[0], X, *projX);
1038 MVT::MvTimesMatAddMv(-ONE, *projVecs_[0], *projX, ONE, Y);
1043 OPT::Apply(*Op_,X,*MX);
1045 std::vector<Scalar> dotprods(nvecsX);
1046 MVT::MvDot(*projVecs_[0], X, dotprods);
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());
1070 projVecs_[0] = helperMV;
Abstract base class for trace minimization eigensolvers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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...
static RCP< Time > getNewTimer(const std::string &name)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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 ...
ParameterList & setParameters(const ParameterList &source)
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.