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
238 virtual void iterate() = 0;
283 bool isInitialized()
const;
300 int getNumIters()
const;
303 void resetNumIters();
312 Teuchos::RCP<const MV> getRitzVectors();
319 std::vector<Value<ScalarType> > getRitzValues();
328 std::vector<int> getRitzIndex();
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();
361 int getCurSubspaceDim()
const;
365 int getMaxSubspaceDim()
const;
376 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest()
const;
390 void setBlockSize(
int blockSize);
394 int getBlockSize()
const;
417 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
420 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs()
const;
428 virtual void currentStatus(std::ostream &os);
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.");
752 V_ =
MVT::Clone(*tmp,numAuxVecs_ + blockSize);
757 BV_ =
MVT::Clone(*tmp,numAuxVecs_ + blockSize);
764 if (hasPrec_ && olsenPrec_) {
766 PBV_ =
MVT::Clone(*tmp,numAuxVecs_ + blockSize);
773 std::vector<int> indV(numAuxVecs_+blockSize);
774 for (
int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
785 if (hasPrec_ && olsenPrec_) {
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;
821 Aeta_ = Teuchos::null;
822 Adelta_ = Teuchos::null;
831 Beta_ = Teuchos::null;
832 Bdelta_ = Teuchos::null;
842 if (hasPrec_ || isSkinny_) {
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;
887 if (hasPrec_ || isSkinny_) {
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;
938 if (hasPrec_ || isSkinny_) {
944 blockSize_ = blockSize;
950 std::vector<int> indX(blockSize_);
951 for (
int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
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) {
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) {
1015 for (
size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
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_);
1029 if (hasPrec_ && olsenPrec_) {
1030 PBV_ =
MVT::Clone(*R_,numAuxVecs_ + blockSize_);
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;
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) {
1106 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
1109 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
1121 if (newstate.
AX != Teuchos::null) {
1123 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
1126 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
1131 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1132 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1135 counterAOp_ += blockSize_;
1138 newstate.
R = Teuchos::null;
1144 if (newstate.
BX != Teuchos::null) {
1146 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
1149 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
1154 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1155 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
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.");
1183 if (initSize > blockSize_) {
1185 initSize = blockSize_;
1186 std::vector<int> ind(blockSize_);
1187 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1193 std::vector<int> ind(initSize);
1194 for (
int i=0; i<initSize; i++) ind[i] = i;
1198 if (blockSize_ > initSize) {
1199 std::vector<int> ind(blockSize_ - initSize);
1200 for (
int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1209 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
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_ );
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_ );
1276 nevLocal_ = blockSize_;
1282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1283 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1288 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1289 Teuchos::TimeMonitor lcltimer( *timerDS_ );
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_);
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_ );
1361 std::vector<int> ind(blockSize_);
1362 for (
int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1368 this->BX_ = this->X_;
1374 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
1377 if (newstate.
R != Teuchos::null) {
1379 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
1381 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
1385 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1386 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1390 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1392 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
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) {
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_) {
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_);
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_) {
1528 os <<
" >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
1530 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
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_);
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_);
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_);
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
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
1648 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
1654 template <
class ScalarType,
class MV,
class OP>
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>
1670 const MV &xi,
const MV &zeta,
1671 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d)
const 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.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
void resetNumIters()
Reset the iteration count.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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.
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.
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.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For RTR, this always returns getBlockSiz...
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 ...
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 Teuchos::ScalarTraits< ScalarType >::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP< const OP > M=Teuchos::null)
Return the maximum coefficient of the matrix scaled by the maximum coefficient of MX...
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...
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.
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...
static int directSolver(int size, const Teuchos::SerialDenseMatrix< int, ScalarType > &KK, Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > MM, Teuchos::SerialDenseMatrix< int, ScalarType > &EV, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &theta, int &nev, int esType=0)
Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK...
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().
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
RTRInitFailure is thrown when the RTR solver is unable to generate an initial iterate in the RTRBase:...
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
Teuchos::ScalarTraits< ScalarType >::magnitudeType rho
The current rho value. This is only valid if the debugging level of verbosity is enabled.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MV > R
The current residual vectors.
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
Types and exceptions used within Anasazi solvers and interfaces.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
RTRState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
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...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
Class which provides internal utilities for the Anasazi solvers.