33 #ifndef ANASAZI_RTRBASE_HPP 34 #define ANASAZI_RTRBASE_HPP 41 #include "Teuchos_ScalarTraits.hpp" 46 #include "Teuchos_LAPACK.hpp" 47 #include "Teuchos_BLAS.hpp" 48 #include "Teuchos_SerialDenseMatrix.hpp" 49 #include "Teuchos_ParameterList.hpp" 50 #include "Teuchos_TimeMonitor.hpp" 110 template <
class ScalarType,
class MV>
113 Teuchos::RCP<const MV>
X;
115 Teuchos::RCP<const MV>
AX;
117 Teuchos::RCP<const MV>
BX;
119 Teuchos::RCP<const MV>
R;
121 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >
T;
125 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
rho;
126 RTRState() :
X(Teuchos::null),
AX(Teuchos::null),
BX(Teuchos::null),
127 R(Teuchos::null),
T(Teuchos::null),
rho(0) {};
188 template <
class ScalarType,
class MV,
class OP>
201 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
205 Teuchos::ParameterList ¶ms,
206 const std::string &solverLabel,
bool skinnySolver
335 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
343 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
350 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
376 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
417 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
420 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const;
439 typedef Teuchos::ScalarTraits<ScalarType> SCT;
440 typedef typename SCT::magnitudeType MagnitudeType;
441 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
442 const MagnitudeType ONE;
443 const MagnitudeType ZERO;
444 const MagnitudeType NANVAL;
445 typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
446 typedef typename std::vector<ScalarType>::iterator vecSTiter;
451 bool checkX, checkAX, checkBX;
452 bool checkEta, checkAEta, checkBEta;
453 bool checkR, checkQ, checkBR;
454 bool checkZ, checkPBX;
455 CheckList() : checkX(false),checkAX(false),checkBX(false),
456 checkEta(false),checkAEta(false),checkBEta(false),
457 checkR(false),checkQ(false),checkBR(false),
458 checkZ(false), checkPBX(false) {};
463 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
465 virtual void solveTRSubproblem() = 0;
467 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(
const MV &xi)
const;
468 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(
const MV &xi,
const MV &zeta)
const;
469 void ginnersep(
const MV &xi, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d)
const;
470 void ginnersep(
const MV &xi,
const MV &zeta, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d)
const;
474 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
475 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
476 const Teuchos::RCP<OutputManager<ScalarType> > om_;
477 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
478 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_;
482 Teuchos::RCP<const OP> AOp_;
483 Teuchos::RCP<const OP> BOp_;
484 Teuchos::RCP<const OP> Prec_;
485 bool hasBOp_, hasPrec_, olsenPrec_;
489 Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
491 timerLocalProj_, timerDS_,
492 timerLocalUpdate_, timerCompRes_,
493 timerOrtho_, timerInit_;
498 int counterAOp_, counterBOp_, counterPrec_;
561 Teuchos::RCP<MV> V_, BV_, PBV_,
564 delta_, Adelta_, Bdelta_, Hdelta_,
566 Teuchos::RCP<const MV> X_, BX_;
569 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
576 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
579 bool Rnorms_current_, R2norms_current_;
582 MagnitudeType conv_kappa_, conv_theta_;
594 template <
class ScalarType,
class MV,
class OP>
597 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
601 Teuchos::ParameterList ¶ms,
602 const std::string &solverLabel,
bool skinnySolver
604 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
605 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
606 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
613 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614 timerAOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation A*x")),
615 timerBOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation B*x")),
616 timerPrec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation Prec*x")),
617 timerSort_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Sorting eigenvalues")),
618 timerLocalProj_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Local projection")),
619 timerDS_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Direct solve")),
620 timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Local update")),
621 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Computing residuals")),
622 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Orthogonalization")),
623 timerInit_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Initialization")),
632 isSkinny_(skinnySolver),
634 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
637 Rnorms_current_(false),
638 R2norms_current_(false),
644 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
645 "Anasazi::RTRBase::constructor: user passed null problem pointer.");
646 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
647 "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
648 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
649 "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
650 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
651 "Anasazi::RTRBase::constructor: user passed null status test pointer.");
652 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
653 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
654 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
655 "Anasazi::RTRBase::constructor: problem is not set.");
656 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
657 "Anasazi::RTRBase::constructor: problem is not Hermitian.");
660 AOp_ = problem_->getOperator();
661 TEUCHOS_TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
662 "Anasazi::RTRBase::constructor: problem provides no A matrix.");
663 BOp_ = problem_->getM();
664 Prec_ = problem_->getPrec();
665 hasBOp_ = (BOp_ != Teuchos::null);
666 hasPrec_ = (Prec_ != Teuchos::null);
667 olsenPrec_ = params.get<
bool>(
"Olsen Prec",
true);
669 TEUCHOS_TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
670 "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
673 int bs = params.get(
"Block Size", problem_->getNEV());
677 leftMost_ = params.get(
"Leftmost",leftMost_);
679 conv_kappa_ = params.get(
"Kappa Convergence",conv_kappa_);
680 TEUCHOS_TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
681 "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
682 conv_theta_ = params.get(
"Theta Convergence",conv_theta_);
683 TEUCHOS_TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
684 "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
690 template <
class ScalarType,
class MV,
class OP>
694 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 695 Teuchos::TimeMonitor lcltimer( *timerInit_ );
702 Teuchos::RCP<const MV> tmp;
708 if (blockSize_ > 0) {
712 tmp = problem_->getInitVec();
713 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
714 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
717 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetGlobalLength(*tmp), std::invalid_argument,
718 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
721 if (blockSize == blockSize_) {
740 if (blockSize > blockSize_)
744 Teuchos::RCP<const MV> Q;
745 std::vector<int> indQ(numAuxVecs_);
746 for (
int i=0; i<numAuxVecs_; i++) indQ[i] = i;
748 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
749 "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
751 if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
752 V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
753 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
756 if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
757 BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
758 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
764 if (hasPrec_ && olsenPrec_) {
765 if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
766 PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
767 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
773 std::vector<int> indV(numAuxVecs_+blockSize);
774 for (
int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
776 V_ = MVT::CloneCopy(*V_,indV);
779 BV_ = MVT::CloneCopy(*BV_,indV);
785 if (hasPrec_ && olsenPrec_) {
786 PBV_ = MVT::CloneCopy(*PBV_,indV);
790 if (blockSize < blockSize_) {
792 blockSize_ = blockSize;
794 theta_.resize(blockSize_);
795 ritz2norms_.resize(blockSize_);
796 Rnorms_.resize(blockSize_);
797 R2norms_.resize(blockSize_);
801 std::vector<int> ind(blockSize_);
802 for (
int i=0; i<blockSize_; i++) ind[i] = i;
810 R_ = MVT::CloneCopy(*R_ ,ind);
811 eta_ = MVT::CloneCopy(*eta_ ,ind);
812 delta_ = MVT::CloneCopy(*delta_ ,ind);
813 Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
815 AX_ = MVT::CloneCopy(*AX_ ,ind);
816 Aeta_ = MVT::CloneCopy(*Aeta_ ,ind);
817 Adelta_ = MVT::CloneCopy(*Adelta_,ind);
821 Aeta_ = Teuchos::null;
822 Adelta_ = Teuchos::null;
827 Beta_ = MVT::CloneCopy(*Beta_,ind);
828 Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
831 Beta_ = Teuchos::null;
832 Bdelta_ = Teuchos::null;
842 if (hasPrec_ || isSkinny_) {
843 Z_ = MVT::Clone(*V_,blockSize_);
855 eta_ = Teuchos::null;
856 Aeta_ = Teuchos::null;
857 delta_ = Teuchos::null;
858 Adelta_ = Teuchos::null;
859 Hdelta_ = Teuchos::null;
860 Beta_ = Teuchos::null;
861 Bdelta_ = Teuchos::null;
864 R_ = MVT::Clone(*tmp,blockSize_);
865 eta_ = MVT::Clone(*tmp,blockSize_);
866 delta_ = MVT::Clone(*tmp,blockSize_);
867 Hdelta_ = MVT::Clone(*tmp,blockSize_);
869 AX_ = MVT::Clone(*tmp,blockSize_);
870 Aeta_ = MVT::Clone(*tmp,blockSize_);
871 Adelta_ = MVT::Clone(*tmp,blockSize_);
876 Beta_ = MVT::Clone(*tmp,blockSize_);
877 Bdelta_ = MVT::Clone(*tmp,blockSize_);
887 if (hasPrec_ || isSkinny_) {
888 Z_ = MVT::Clone(*tmp,blockSize_);
897 initialized_ =
false;
900 theta_.resize(blockSize,NANVAL);
901 ritz2norms_.resize(blockSize,NANVAL);
902 Rnorms_.resize(blockSize,NANVAL);
903 R2norms_.resize(blockSize,NANVAL);
908 eta_ = Teuchos::null;
909 Aeta_ = Teuchos::null;
910 delta_ = Teuchos::null;
911 Adelta_ = Teuchos::null;
912 Hdelta_ = Teuchos::null;
913 Beta_ = Teuchos::null;
914 Bdelta_ = Teuchos::null;
918 R_ = MVT::Clone(*tmp,blockSize);
919 eta_ = MVT::Clone(*tmp,blockSize);
920 delta_ = MVT::Clone(*tmp,blockSize);
921 Hdelta_ = MVT::Clone(*tmp,blockSize);
923 AX_ = MVT::Clone(*tmp,blockSize);
924 Aeta_ = MVT::Clone(*tmp,blockSize);
925 Adelta_ = MVT::Clone(*tmp,blockSize);
930 Beta_ = MVT::Clone(*tmp,blockSize);
931 Bdelta_ = MVT::Clone(*tmp,blockSize);
938 if (hasPrec_ || isSkinny_) {
939 Z_ = MVT::Clone(*tmp,blockSize);
944 blockSize_ = blockSize;
950 std::vector<int> indX(blockSize_);
951 for (
int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
952 X_ = MVT::CloneView(*V_,indX);
954 BX_ = MVT::CloneView(*BV_,indX);
965 template <
class ScalarType,
class MV,
class OP>
967 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
968 "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
975 template <
class ScalarType,
class MV,
class OP>
983 template <
class ScalarType,
class MV,
class OP>
985 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
989 auxVecs_.reserve(auxvecs.size());
992 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
993 numAuxVecs_ += MVT::GetNumberVecs(**v);
997 if (numAuxVecs_ > 0 && initialized_) {
998 initialized_ =
false;
1003 BX_ = Teuchos::null;
1006 BV_ = Teuchos::null;
1007 PBV_ = Teuchos::null;
1010 if (numAuxVecs_ > 0) {
1011 V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1013 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1014 std::vector<int> ind(MVT::GetNumberVecs(**v));
1015 for (
size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
1016 MVT::SetBlock(**v,ind,*V_);
1017 auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
1019 TEUCHOS_TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
1020 "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
1023 BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1024 OPT::Apply(*BOp_,*V_,*BV_);
1029 if (hasPrec_ && olsenPrec_) {
1030 PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1031 OPT::Apply(*Prec_,*BV_,*V_);
1036 if (om_->isVerbosity(
Debug ) ) {
1040 om_->print(
Debug, accuracyCheck(chk,
"in setAuxVecs()") );
1057 template <
class ScalarType,
class MV,
class OP>
1066 BX_ = Teuchos::null;
1067 Teuchos::RCP<MV> X, BX;
1069 std::vector<int> ind(blockSize_);
1070 for (
int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1071 X = MVT::CloneViewNonConst(*V_,ind);
1073 BX = MVT::CloneViewNonConst(*BV_,ind);
1080 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1081 Teuchos::TimeMonitor inittimer( *timerInit_ );
1084 std::vector<int> bsind(blockSize_);
1085 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
1104 if (newstate.
X != Teuchos::null) {
1105 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X),
1106 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
1108 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
X) < blockSize_,
1109 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
1112 MVT::SetBlock(*newstate.
X,bsind,*X);
1121 if (newstate.
AX != Teuchos::null) {
1122 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
AX) != MVT::GetGlobalLength(*X),
1123 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
1125 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
AX) < blockSize_,
1126 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
1127 MVT::SetBlock(*newstate.
AX,bsind,*AX_);
1131 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1132 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1134 OPT::Apply(*AOp_,*X,*AX_);
1135 counterAOp_ += blockSize_;
1138 newstate.
R = Teuchos::null;
1144 if (newstate.
BX != Teuchos::null) {
1145 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
BX) != MVT::GetGlobalLength(*X),
1146 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
1148 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
BX) < blockSize_,
1149 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
1150 MVT::SetBlock(*newstate.
BX,bsind,*BX);
1154 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1155 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1157 OPT::Apply(*BOp_,*X,*BX);
1158 counterBOp_ += blockSize_;
1161 newstate.
R = Teuchos::null;
1166 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error,
"Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1174 newstate.
R = Teuchos::null;
1175 newstate.
T = Teuchos::null;
1178 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1179 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1180 "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1182 int initSize = MVT::GetNumberVecs(*ivec);
1183 if (initSize > blockSize_) {
1185 initSize = blockSize_;
1186 std::vector<int> ind(blockSize_);
1187 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1188 ivec = MVT::CloneView(*ivec,ind);
1193 std::vector<int> ind(initSize);
1194 for (
int i=0; i<initSize; i++) ind[i] = i;
1195 MVT::SetBlock(*ivec,ind,*X);
1198 if (blockSize_ > initSize) {
1199 std::vector<int> ind(blockSize_ - initSize);
1200 for (
int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1201 Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind);
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1209 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1211 OPT::Apply(*BOp_,*X,*BX);
1212 counterBOp_ += blockSize_;
1216 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error,
"Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1220 if (numAuxVecs_ > 0) {
1221 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1222 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1224 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1225 int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
1227 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1230 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1231 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1233 int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
1235 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1243 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1244 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1246 OPT::Apply(*AOp_,*X,*AX_);
1247 counterAOp_ += blockSize_;
1254 if (newstate.
T != Teuchos::null) {
1255 TEUCHOS_TEST_FOR_EXCEPTION( (
signed int)(newstate.
T->size()) < blockSize_,
1256 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1257 for (
int i=0; i<blockSize_; i++) {
1258 theta_[i] = (*newstate.
T)[i];
1263 Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
1264 BB(blockSize_,blockSize_),
1265 S(blockSize_,blockSize_);
1268 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1269 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1271 MVT::MvTransMv(ONE,*X,*AX_,AA);
1273 MVT::MvTransMv(ONE,*X,*BX,BB);
1276 nevLocal_ = blockSize_;
1282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1283 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1285 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
1288 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1289 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1291 ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
1294 "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
1295 TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,
RTRInitFailure,
"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
1300 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1301 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1304 std::vector<int> order(blockSize_);
1307 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1310 Utils::permuteVectors(order,S);
1315 Teuchos::BLAS<int,ScalarType> blas;
1316 Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
1320 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
1321 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1326 for (
int i=0; i<blockSize_; i++) {
1327 blas.SCAL(blockSize_,theta_[i],RR[i],1);
1329 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
1330 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1331 for (
int i=0; i<blockSize_; i++) {
1332 ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
1340 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1341 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1344 MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );
1345 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
1347 MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );
1348 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
1351 MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );
1352 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
1361 std::vector<int> ind(blockSize_);
1362 for (
int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1363 this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
1365 this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
1368 this->BX_ = this->X_;
1374 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
1377 if (newstate.
R != Teuchos::null) {
1378 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) < blockSize_,
1379 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
1380 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
1381 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
1382 MVT::SetBlock(*newstate.
R,bsind,*R_);
1385 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1386 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1389 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
1390 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1392 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1393 MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
1397 Rnorms_current_ =
false;
1398 R2norms_current_ =
false;
1403 AX_ = Teuchos::null;
1407 initialized_ =
true;
1409 if (om_->isVerbosity(
Debug ) ) {
1417 om_->print(
Debug, accuracyCheck(chk,
"after initialize()") );
1421 template <
class ScalarType,
class MV,
class OP>
1433 template <
class ScalarType,
class MV,
class OP>
1434 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1436 if (Rnorms_current_ ==
false) {
1438 orthman_->norm(*R_,Rnorms_);
1439 Rnorms_current_ =
true;
1447 template <
class ScalarType,
class MV,
class OP>
1448 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1450 if (R2norms_current_ ==
false) {
1452 MVT::MvNorm(*R_,R2norms_);
1453 R2norms_current_ =
true;
1485 template <
class ScalarType,
class MV,
class OP>
1488 using std::setprecision;
1489 using std::scientific;
1492 std::stringstream os;
1495 os <<
" Debugging checks: " << where << endl;
1498 if (chk.checkX && initialized_) {
1499 tmp = orthman_->orthonormError(*X_);
1500 os <<
" >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl;
1501 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1502 tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
1503 os <<
" >> Error in X^H B Q[" << i <<
"] == 0 : " << scientific << setprecision(10) << tmp << endl;
1506 if (chk.checkAX && !isSkinny_ && initialized_) {
1507 tmp = Utils::errorEquality(*X_, *AX_, AOp_);
1508 os <<
" >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl;
1510 if (chk.checkBX && hasBOp_ && initialized_) {
1511 Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
1512 OPT::Apply(*BOp_,*this->X_,*BX);
1513 tmp = Utils::errorEquality(*BX, *BX_);
1514 os <<
" >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl;
1518 if (chk.checkEta && initialized_) {
1519 tmp = orthman_->orthogError(*X_,*eta_);
1520 os <<
" >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl;
1521 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1522 tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
1523 os <<
" >> Error in Eta^H B Q[" << i <<
"]==0 : " << scientific << setprecision(10) << tmp << endl;
1526 if (chk.checkAEta && !isSkinny_ && initialized_) {
1527 tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
1528 os <<
" >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
1530 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
1531 tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
1532 os <<
" >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl;
1536 if (chk.checkR && initialized_) {
1537 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1538 MVT::MvTransMv(ONE,*X_,*R_,xTx);
1539 tmp = xTx.normFrobenius();
1540 os <<
" >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl;
1545 if (chk.checkBR && hasBOp_ && initialized_) {
1546 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1547 MVT::MvTransMv(ONE,*BX_,*R_,xTx);
1548 tmp = xTx.normFrobenius();
1549 os <<
" >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl;
1553 if (chk.checkZ && initialized_) {
1554 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1555 MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
1556 tmp = xTx.normFrobenius();
1557 os <<
" >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl;
1562 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1563 tmp = orthman_->orthonormError(*auxVecs_[i]);
1564 os <<
" >> Error in Q[" << i <<
"]^H B Q[" << i <<
"]==I: " << scientific << setprecision(10) << tmp << endl;
1565 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1566 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1567 os <<
" >> Error in Q[" << i <<
"]^H B Q[" << j <<
"]==0: " << scientific << setprecision(10) << tmp << endl;
1578 template <
class ScalarType,
class MV,
class OP>
1582 using std::setprecision;
1583 using std::scientific;
1588 os <<
"================================================================================" << endl;
1590 os <<
" RTRBase Solver Status" << endl;
1592 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
1593 os <<
"The number of iterations performed is " << iter_ << endl;
1594 os <<
"The current block size is " << blockSize_ << endl;
1595 os <<
"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1596 os <<
"The number of operations A*x is " << counterAOp_ << endl;
1597 os <<
"The number of operations B*x is " << counterBOp_ << endl;
1598 os <<
"The number of operations Prec*x is " << counterPrec_ << endl;
1599 os <<
"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
1600 os <<
"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
1604 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1605 os << setw(20) <<
"Eigenvalue" 1606 << setw(20) <<
"Residual(B)" 1607 << setw(20) <<
"Residual(2)" 1609 os <<
"--------------------------------------------------------------------------------"<<endl;
1610 for (
int i=0; i<blockSize_; i++) {
1611 os << scientific << setprecision(10) << setw(20) << theta_[i];
1612 if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
1613 else os << scientific << setprecision(10) << setw(20) <<
"not current";
1614 if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
1615 else os << scientific << setprecision(10) << setw(20) <<
"not current";
1619 os <<
"================================================================================" << endl;
1626 template <
class ScalarType,
class MV,
class OP>
1627 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1630 std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
1632 MagnitudeType ret = 0;
1633 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1642 template <
class ScalarType,
class MV,
class OP>
1643 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1644 RTRBase<ScalarType,MV,OP>::ginner(
const MV &xi,
const MV &zeta)
const 1646 std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
1647 MVT::MvDot(xi,zeta,d);
1648 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
1654 template <
class ScalarType,
class MV,
class OP>
1655 void RTRBase<ScalarType,MV,OP>::ginnersep(
1657 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d)
const 1660 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1668 template <
class ScalarType,
class MV,
class OP>
1669 void RTRBase<ScalarType,MV,OP>::ginnersep(
1670 const MV &xi,
const MV &zeta,
1671 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d)
const 1673 std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
1674 MVT::MvDot(xi,zeta,dC);
1675 vecMTiter iR=d.begin();
1676 vecSTiter iS=dC.begin();
1677 for (; iR != d.end(); iR++, iS++) {
1678 (*iR) = SCT::real(*iS);
1682 template <
class ScalarType,
class MV,
class OP>
1687 template <
class ScalarType,
class MV,
class OP>
1692 template <
class ScalarType,
class MV,
class OP>
1697 template <
class ScalarType,
class MV,
class OP>
1702 template <
class ScalarType,
class MV,
class OP>
1705 if (!initialized_)
return 0;
1709 template <
class ScalarType,
class MV,
class OP>
1710 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1713 std::vector<MagnitudeType> ret = ritz2norms_;
1714 ret.resize(nevLocal_);
1718 template <
class ScalarType,
class MV,
class OP>
1719 std::vector<Value<ScalarType> >
1722 std::vector<Value<ScalarType> > ret(nevLocal_);
1723 for (
int i=0; i<nevLocal_; i++) {
1724 ret[i].realpart = theta_[i];
1725 ret[i].imagpart = ZERO;
1730 template <
class ScalarType,
class MV,
class OP>
1731 Teuchos::RCP<const MV>
1737 template <
class ScalarType,
class MV,
class OP>
1743 template <
class ScalarType,
class MV,
class OP>
1749 template <
class ScalarType,
class MV,
class OP>
1759 state.
BX = Teuchos::null;
1763 state.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta_));
1767 template <
class ScalarType,
class MV,
class OP>
1770 return initialized_;
1773 template <
class ScalarType,
class MV,
class OP>
1776 std::vector<int> ret(nevLocal_,0);
1783 #endif // ANASAZI_RTRBASE_HPP RTROrthoFailure is thrown when an orthogonalization attempt fails.
void resetNumIters()
Reset the iteration count.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
Structure to contain pointers to RTR state variables.
virtual void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
Teuchos::RCP< const MV > AX
The image of the current eigenvectors under A, or Teuchos::null is we implement a skinny solver...
This class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > BX
The image of the current eigenvectors under B, or Teuchos::null if B was not specified.
Declaration of basic traits for the multivector type.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the Ritz residuals.
RTRState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers. The class provides the interfaces shared by the IRTR solvers (e.g., getState() and initialize()) as well as the shared implementations (e.g., inner products).
Virtual base class which defines basic traits for the operator type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
virtual void iterate()=0
This method performs RTR iterations until the status test indicates the need to stop or an error occu...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For RTR, this always returns getBlockSiz...
Teuchos::RCP< const MV > X
The current eigenvectors.
An exception class parent to all Anasazi exceptions.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
virtual ~RTRBase()
RTRBase destructor
Anasazi's templated, static class providing utilities for the solvers.
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
RTRInitFailure is thrown when the RTR solver is unable to generate an initial iterate in the RTRBase:...
Teuchos::ScalarTraits< ScalarType >::magnitudeType rho
The current rho value. This is only valid if the debugging level of verbosity is enabled.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MV > R
The current residual vectors.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
Types and exceptions used within Anasazi solvers and interfaces.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
RTRBase(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms, const std::string &solverLabel, bool skinnySolver)
RTRBase constructor with eigenproblem, solver utilities, and parameter list of solver options...
Class which provides internal utilities for the Anasazi solvers.