33 #ifndef ANASAZI_RTRBASE_HPP 34 #define ANASAZI_RTRBASE_HPP 110 template <
class ScalarType,
class MV>
188 template <
class ScalarType,
class MV,
class OP>
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();
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;
485 bool hasBOp_, hasPrec_, olsenPrec_;
491 timerLocalProj_, timerDS_,
492 timerLocalUpdate_, timerCompRes_,
493 timerOrtho_, timerInit_;
498 int counterAOp_, counterBOp_, counterPrec_;
564 delta_, Adelta_, Bdelta_, Hdelta_,
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>
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),
637 Rnorms_current_(false),
638 R2norms_current_(false),
645 "Anasazi::RTRBase::constructor: user passed null problem pointer.");
647 "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
649 "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
651 "Anasazi::RTRBase::constructor: user passed null status test pointer.");
653 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
655 "Anasazi::RTRBase::constructor: problem is not set.");
657 "Anasazi::RTRBase::constructor: problem is not Hermitian.");
660 AOp_ = problem_->getOperator();
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);
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_);
681 "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
682 conv_theta_ = params.
get(
"Theta Convergence",conv_theta_);
684 "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
690 template <
class ScalarType,
class MV,
class OP>
694 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 708 if (blockSize_ > 0) {
712 tmp = problem_->getInitVec();
714 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
718 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
721 if (blockSize == blockSize_) {
740 if (blockSize > blockSize_)
745 std::vector<int> indQ(numAuxVecs_);
746 for (
int i=0; i<numAuxVecs_; i++) indQ[i] = i;
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>
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>
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));
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;
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 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.");
1112 MVT::SetBlock(*newstate.
X,bsind,*X);
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.");
1127 MVT::SetBlock(*newstate.
AX,bsind,*AX_);
1131 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1134 OPT::Apply(*AOp_,*X,*AX_);
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.");
1150 MVT::SetBlock(*newstate.
BX,bsind,*BX);
1154 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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;
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;
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 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 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 1246 OPT::Apply(*AOp_,*X,*AX_);
1247 counterAOp_ += blockSize_;
1254 if (newstate.
T != Teuchos::null) {
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];
1264 BB(blockSize_,blockSize_),
1265 S(blockSize_,blockSize_);
1268 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1271 MVT::MvTransMv(ONE,*X,*AX_,AA);
1273 MVT::MvTransMv(ONE,*X,*BX,BB);
1276 nevLocal_ = blockSize_;
1282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1285 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
1288 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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);
1300 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1304 std::vector<int> order(blockSize_);
1307 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1310 Utils::permuteVectors(order,S);
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);
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 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) {
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." );
1382 MVT::SetBlock(*newstate.
R,bsind,*R_);
1385 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1389 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
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_) {
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_) {
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_) {
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_) {
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>
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>
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(
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,
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>
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.
T & get(const std::string &name, T def_value)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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...
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
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.
void resize(size_type new_size, const value_type &x=value_type())
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...
void SCAL(const OrdinalType n, const ScalarType alpha, ScalarType *x, const OrdinalType incx) const
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...
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
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.