52 #ifndef ANASAZI_LOBPCG_HPP 53 #define ANASAZI_LOBPCG_HPP 100 template <
class ScalarType,
class MultiVector>
208 template <
class ScalarType,
class MV,
class OP>
373 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
381 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
391 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
496 const MagnitudeType ONE;
497 const MagnitudeType ZERO;
498 const MagnitudeType NANVAL;
503 bool checkX, checkMX, checkKX;
504 bool checkH, checkMH;
505 bool checkP, checkMP, checkKP;
507 CheckList() : checkX(false),checkMX(false),checkKX(false),
508 checkH(false),checkMH(false),
509 checkP(false),checkMP(false),checkKP(false),
510 checkR(false),checkQ(false) {};
515 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
536 timerLocalProj_, timerDS_,
537 timerLocalUpdate_, timerCompRes_,
538 timerOrtho_, timerInit_;
543 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
598 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
601 bool Rnorms_current_, R2norms_current_;
610 template <
class ScalarType,
class MV,
class OP>
619 ONE(
Teuchos::ScalarTraits<MagnitudeType>::one()),
620 ZERO(
Teuchos::ScalarTraits<MagnitudeType>::zero()),
621 NANVAL(
Teuchos::ScalarTraits<MagnitudeType>::nan()),
629 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
630 timerOp_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Operation Op*x")),
631 timerMOp_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Operation M*x")),
632 timerPrec_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Operation Prec*x")),
633 timerSort_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Sorting eigenvalues")),
634 timerLocalProj_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Local projection")),
635 timerDS_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Direct solve")),
636 timerLocalUpdate_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Local update")),
637 timerCompRes_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Computing residuals")),
638 timerOrtho_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Orthogonalization")),
639 timerInit_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCG::Initialization")),
646 fullOrtho_(params.get(
"Full Ortho", true)),
653 Rnorms_current_(false),
654 R2norms_current_(false)
657 "Anasazi::LOBPCG::constructor: user passed null problem pointer.");
659 "Anasazi::LOBPCG::constructor: user passed null sort manager pointer.");
661 "Anasazi::LOBPCG::constructor: user passed null output manager pointer.");
663 "Anasazi::LOBPCG::constructor: user passed null status test pointer.");
665 "Anasazi::LOBPCG::constructor: user passed null orthogonalization manager pointer.");
667 "Anasazi::LOBPCG::constructor: problem is not set.");
669 "Anasazi::LOBPCG::constructor: problem is not Hermitian; LOBPCG requires Hermitian problem.");
672 Op_ = problem_->getOperator();
674 "Anasazi::LOBPCG::constructor: problem provides no operator.");
675 MOp_ = problem_->getM();
676 Prec_ = problem_->getPrec();
677 hasM_ = (MOp_ != Teuchos::null);
680 int bs = params.
get(
"Block Size", problem_->getNEV());
687 template <
class ScalarType,
class MV,
class OP>
691 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 705 if (blockSize_ > 0) {
709 tmp = problem_->getInitVec();
711 "Anasazi::LOBPCG::setBlockSize(): eigenproblem did not specify initial vectors to clone from.");
715 std::invalid_argument,
"Anasazi::LOBPCG::setBlockSize(): block size must be strictly positive.");
716 if (newBS == blockSize_) {
720 else if (newBS < blockSize_ && initialized_) {
736 std::vector<int> newind(newBS), oldind(newBS);
737 for (
int i=0; i<newBS; i++) {
745 newR = MVT::Clone(*tmp,newBS);
746 newV = MVT::Clone(*tmp,newBS*3);
747 newKV = MVT::Clone(*tmp,newBS*3);
749 newMV = MVT::Clone(*tmp,newBS*3);
769 theta_.resize(3*newBS);
770 Rnorms_.resize(newBS);
771 R2norms_.resize(newBS);
774 src = MVT::CloneView(*R_,newind);
775 MVT::SetBlock(*src,newind,*newR);
781 src = MVT::CloneView(*V_,oldind);
782 MVT::SetBlock(*src,newind,*newV);
783 src = MVT::CloneView(*KV_,oldind);
784 MVT::SetBlock(*src,newind,*newKV);
786 src = MVT::CloneView(*MV_,oldind);
787 MVT::SetBlock(*src,newind,*newMV);
790 for (
int i=0; i<newBS; i++) {
791 newind[i] += 2*newBS;
792 oldind[i] += 2*blockSize_;
794 src = MVT::CloneView(*V_,oldind);
795 MVT::SetBlock(*src,newind,*newV);
796 src = MVT::CloneView(*KV_,oldind);
797 MVT::SetBlock(*src,newind,*newKV);
799 src = MVT::CloneView(*MV_,oldind);
800 MVT::SetBlock(*src,newind,*newMV);
819 initialized_ =
false;
838 theta_.resize(3*newBS,NANVAL);
839 Rnorms_.resize(newBS,NANVAL);
840 R2norms_.resize(newBS,NANVAL);
843 R_ = MVT::Clone(*tmp,newBS);
844 V_ = MVT::Clone(*tmp,newBS*3);
845 KV_ = MVT::Clone(*tmp,newBS*3);
847 MV_ = MVT::Clone(*tmp,newBS*3);
855 tmpmvec_ = Teuchos::null;
857 tmpmvec_ = MVT::Clone(*tmp,newBS);
870 template <
class ScalarType,
class MV,
class OP>
873 std::vector<int> ind(blockSize_);
875 for (
int i=0; i<blockSize_; i++) {
878 X_ = MVT::CloneViewNonConst(*V_,ind);
879 KX_ = MVT::CloneViewNonConst(*KV_,ind);
881 MX_ = MVT::CloneViewNonConst(*MV_,ind);
887 for (
int i=0; i<blockSize_; i++) {
888 ind[i] += blockSize_;
890 H_ = MVT::CloneViewNonConst(*V_,ind);
891 KH_ = MVT::CloneViewNonConst(*KV_,ind);
893 MH_ = MVT::CloneViewNonConst(*MV_,ind);
899 for (
int i=0; i<blockSize_; i++) {
900 ind[i] += blockSize_;
902 P_ = MVT::CloneViewNonConst(*V_,ind);
903 KP_ = MVT::CloneViewNonConst(*KV_,ind);
905 MP_ = MVT::CloneViewNonConst(*MV_,ind);
915 template <
class ScalarType,
class MV,
class OP>
923 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
924 numAuxVecs_ += MVT::GetNumberVecs(**i);
928 if (numAuxVecs_ > 0 && initialized_) {
929 initialized_ =
false;
933 if (om_->isVerbosity(
Debug ) ) {
937 om_->print(
Debug, accuracyCheck(chk,
": in setAuxVecs()") );
960 template <
class ScalarType,
class MV,
class OP>
966 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 970 std::vector<int> bsind(blockSize_);
971 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
995 if (newstate.
X != Teuchos::null) {
997 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.X not correct." );
1000 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.X must have at least block size vectors.");
1003 MVT::SetBlock(*newstate.
X,bsind,*X_);
1007 if (newstate.
MX != Teuchos::null) {
1009 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MX not correct." );
1012 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.MX must have at least block size vectors.");
1013 MVT::SetBlock(*newstate.
MX,bsind,*MX_);
1018 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1021 OPT::Apply(*MOp_,*X_,*MX_);
1022 count_ApplyM_ += blockSize_;
1025 newstate.
R = Teuchos::null;
1030 if (newstate.
KX != Teuchos::null) {
1032 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KX not correct." );
1035 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.KX must have at least block size vectors.");
1036 MVT::SetBlock(*newstate.
KX,bsind,*KX_);
1041 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1044 OPT::Apply(*Op_,*X_,*KX_);
1045 count_ApplyOp_ += blockSize_;
1048 newstate.
R = Teuchos::null;
1056 newstate.
P = Teuchos::null;
1057 newstate.
KP = Teuchos::null;
1058 newstate.
MP = Teuchos::null;
1059 newstate.
R = Teuchos::null;
1060 newstate.
T = Teuchos::null;
1065 "Anasazi::LOBPCG::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1067 int initSize = MVT::GetNumberVecs(*ivec);
1068 if (initSize > blockSize_) {
1070 initSize = blockSize_;
1071 std::vector<int> ind(blockSize_);
1072 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1073 ivec = MVT::CloneView(*ivec,ind);
1078 std::vector<int> ind(initSize);
1079 for (
int i=0; i<initSize; i++) ind[i] = i;
1080 MVT::SetBlock(*ivec,ind,*X_);
1083 if (blockSize_ > initSize) {
1084 std::vector<int> ind(blockSize_ - initSize);
1085 for (
int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1093 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1096 OPT::Apply(*MOp_,*X_,*MX_);
1097 count_ApplyM_ += blockSize_;
1101 if (numAuxVecs_ > 0) {
1102 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1106 int rank = orthman_->projectAndNormalizeMat(*X_,auxVecs_,dummy,Teuchos::null,MX_);
1108 "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1111 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1114 int rank = orthman_->normalizeMat(*X_,Teuchos::null,MX_);
1116 "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1121 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1124 OPT::Apply(*Op_,*X_,*KX_);
1125 count_ApplyOp_ += blockSize_;
1133 theta_.resize(3*blockSize_,NANVAL);
1134 if (newstate.
T != Teuchos::null) {
1136 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1137 for (
int i=0; i<blockSize_; i++) {
1138 theta_[i] = (*newstate.
T)[i];
1140 nevLocal_ = blockSize_;
1145 MM(blockSize_,blockSize_),
1146 S(blockSize_,blockSize_);
1148 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1152 MVT::MvTransMv(ONE,*X_,*KX_,KK);
1154 MVT::MvTransMv(ONE,*X_,*MX_,MM);
1155 nevLocal_ = blockSize_;
1160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1163 Utils::directSolver(blockSize_, KK, Teuchos::rcpFromRef(MM), S, theta_, nevLocal_, 1);
1165 "Anasazi::LOBPCG::initialize(): Initial Ritz analysis did not produce enough Ritz pairs to initialize algorithm.");
1171 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1175 std::vector<int> order(blockSize_);
1178 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1181 Utils::permuteVectors(order,S);
1186 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1190 MVT::MvAddMv( ONE, *X_, ZERO, *X_, *R_ );
1191 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X_ );
1193 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1194 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *KX_ );
1197 MVT::MvAddMv( ONE, *MX_, ZERO, *MX_, *R_ );
1198 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *MX_ );
1206 if (newstate.
R != Teuchos::null) {
1208 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.R not correct." );
1210 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.R must have blockSize number of vectors." );
1211 MVT::SetBlock(*newstate.
R,bsind,*R_);
1214 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1218 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1220 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1221 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1225 Rnorms_current_ =
false;
1226 R2norms_current_ =
false;
1229 if (newstate.
P != Teuchos::null) {
1231 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.P must have blockSize number of vectors." );
1233 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.P not correct." );
1237 MVT::SetBlock(*newstate.
P,bsind,*P_);
1240 if (newstate.
KP != Teuchos::null) {
1242 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.KP must have blockSize number of vectors." );
1244 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KP not correct." );
1245 MVT::SetBlock(*newstate.
KP,bsind,*KP_);
1248 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1251 OPT::Apply(*Op_,*P_,*KP_);
1252 count_ApplyOp_ += blockSize_;
1257 if (newstate.
MP != Teuchos::null) {
1259 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.MP must have blockSize number of vectors." );
1261 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MP not correct." );
1262 MVT::SetBlock(*newstate.
MP,bsind,*MP_);
1265 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1268 OPT::Apply(*MOp_,*P_,*MP_);
1269 count_ApplyM_ += blockSize_;
1278 initialized_ =
true;
1280 if (om_->isVerbosity(
Debug ) ) {
1291 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
1296 template <
class ScalarType,
class MV,
class OP>
1306 template <
class ScalarType,
class MV,
class OP>
1309 if ( fullOrtho_ ==
true || initialized_ ==
false || fullOrtho == fullOrtho_ ) {
1311 fullOrtho_ = fullOrtho;
1323 if (fullOrtho_ && tmpmvec_ == Teuchos::null) {
1325 tmpmvec_ = MVT::Clone(*X_,blockSize_);
1327 else if (fullOrtho_==
false) {
1329 tmpmvec_ = Teuchos::null;
1336 template <
class ScalarType,
class MV,
class OP>
1342 if (initialized_ ==
false) {
1348 const int oneBlock = blockSize_;
1349 const int twoBlocks = 2*blockSize_;
1350 const int threeBlocks = 3*blockSize_;
1352 std::vector<int> indblock1(blockSize_), indblock2(blockSize_), indblock3(blockSize_);
1353 for (
int i=0; i<blockSize_; i++) {
1355 indblock2[i] = i + blockSize_;
1356 indblock3[i] = i + 2*blockSize_;
1365 MM( threeBlocks, threeBlocks ),
1366 S( threeBlocks, threeBlocks );
1368 while (tester_->checkStatus(
this) !=
Passed) {
1371 if (om_->isVerbosity(
Debug)) {
1372 currentStatus( om_->stream(
Debug) );
1382 if (Prec_ != Teuchos::null) {
1383 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1386 OPT::Apply( *Prec_, *R_, *H_ );
1387 count_ApplyPrec_ += blockSize_;
1390 MVT::MvAddMv(ONE,*R_,ZERO,*R_,*H_);
1395 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1398 OPT::Apply( *MOp_, *H_, *MH_);
1399 count_ApplyM_ += blockSize_;
1415 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1419 Teuchos::tuple<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >(Teuchos::null);
1420 int rank = orthman_->projectAndNormalizeMat(*H_,Q,dummyC,Teuchos::null,MH_);
1423 "Anasazi::LOBPCG::iterate(): unable to compute orthonormal basis for H.");
1426 if (om_->isVerbosity(
Debug ) ) {
1430 om_->print(
Debug, accuracyCheck(chk,
": after ortho H") );
1436 om_->print(
OrthoDetails, accuracyCheck(chk,
": after ortho H") );
1441 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1444 OPT::Apply( *Op_, *H_, *KH_);
1445 count_ApplyOp_ += blockSize_;
1449 nevLocal_ = threeBlocks;
1452 nevLocal_ = twoBlocks;
1474 KX_ = Teuchos::null;
1475 MX_ = Teuchos::null;
1477 KH_ = Teuchos::null;
1478 MH_ = Teuchos::null;
1480 KP_ = Teuchos::null;
1481 MP_ = Teuchos::null;
1482 Teuchos::RCP<const MV> cX, cH, cXHP, cHP, cK_XHP, cK_HP, cM_XHP, cM_HP, cP, cK_P, cM_P;
1484 cX = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock1);
1485 cH = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock2);
1487 std::vector<int> indXHP(nevLocal_);
1488 for (
int i=0; i<nevLocal_; i++) {
1491 cXHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indXHP);
1492 cK_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indXHP);
1494 cM_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indXHP);
1500 std::vector<int> indHP(nevLocal_-blockSize_);
1501 for (
int i=blockSize_; i<nevLocal_; i++) {
1502 indHP[i-blockSize_] = i;
1504 cHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indHP);
1505 cK_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indHP);
1507 cM_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indHP);
1513 if (nevLocal_ == threeBlocks) {
1514 cP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock3);
1515 cK_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indblock3);
1517 cM_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indblock3);
1540 K1(
Teuchos::View,KK,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1541 K2(
Teuchos::View,KK,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1542 K3(
Teuchos::View,KK,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_),
1543 M1(
Teuchos::View,MM,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1544 M2(
Teuchos::View,MM,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1545 M3(
Teuchos::View,MM,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_);
1547 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1550 MVT::MvTransMv( ONE, *cX, *cK_XHP, K1 );
1551 MVT::MvTransMv( ONE, *cX, *cM_XHP, M1 );
1552 MVT::MvTransMv( ONE, *cH, *cK_HP , K2 );
1553 MVT::MvTransMv( ONE, *cH, *cM_HP , M2 );
1554 if (nevLocal_ == threeBlocks) {
1555 MVT::MvTransMv( ONE, *cP, *cK_P, K3 );
1556 MVT::MvTransMv( ONE, *cP, *cM_P, M3 );
1566 cK_P = Teuchos::null;
1567 cM_P = Teuchos::null;
1568 if (fullOrtho_ ==
true) {
1569 cHP = Teuchos::null;
1570 cK_HP = Teuchos::null;
1571 cM_HP = Teuchos::null;
1584 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1587 int localSize = nevLocal_;
1588 Utils::directSolver(localSize, lclKK, Teuchos::rcpFromRef(lclMM), lclS, theta_, nevLocal_, 0);
1597 if (nevLocal_ != localSize) {
1600 cXHP = Teuchos::null;
1601 cK_XHP = Teuchos::null;
1602 cM_XHP = Teuchos::null;
1603 cHP = Teuchos::null;
1604 cK_HP = Teuchos::null;
1605 cM_HP = Teuchos::null;
1609 "Anasazi::LOBPCG::iterate(): indefiniteness detected in projected mass matrix." );
1619 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1623 std::vector<int> order(nevLocal_);
1626 sm_->sort(theta_, Teuchos::rcpFromRef(order), nevLocal_);
1629 Utils::permuteVectors(order,lclS);
1653 MMC(nevLocal_,twoBlocks),
1654 CMMC(twoBlocks ,twoBlocks);
1657 for (
int j=0; j<oneBlock; j++) {
1658 for (
int i=0; i<oneBlock; i++) {
1662 C(i,j+oneBlock) = ZERO;
1666 for (
int j=0; j<oneBlock; j++) {
1667 for (
int i=oneBlock; i<twoBlocks; i++) {
1671 C(i,j+oneBlock) = lclS(i,j);
1675 if (nevLocal_ == threeBlocks) {
1676 for (
int j=0; j<oneBlock; j++) {
1677 for (
int i=twoBlocks; i<threeBlocks; i++) {
1681 C(i,j+oneBlock) = lclS(i,j);
1691 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1695 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1708 MagnitudeType symNorm = K22_trans.
normOne();
1709 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
1711 cXHP = Teuchos::null;
1712 cHP = Teuchos::null;
1713 cK_XHP = Teuchos::null;
1714 cK_HP = Teuchos::null;
1715 cM_XHP = Teuchos::null;
1716 cM_HP = Teuchos::null;
1719 std::string errMsg =
"Anasazi::LOBPCG::iterate(): Cholesky factorization failed during full orthogonalization.";
1720 if ( symNorm < tol )
1726 errMsg +=
" Check the operator A (or K), it may not be Hermitian!";
1733 nevLocal_,twoBlocks,ONE,CMMC.
values(),CMMC.
stride(),C.values(),C.stride());
1740 if (om_->isVerbosity(
Debug ) ) {
1742 tmp2(oneBlock,oneBlock);
1745 std::stringstream os;
1747 os.setf(std::ios::scientific, std::ios::floatfield);
1749 os <<
" Checking Full Ortho: iteration " << iter_ << std::endl;
1755 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1759 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1761 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1763 os <<
" >> Error in CX^H MM CX == I : " << tmp << std::endl;
1769 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1773 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1775 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1777 os <<
" >> Error in CP^H MM CP == I : " << tmp << std::endl;
1782 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1785 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1788 os <<
" >> Error in CX^H MM CP == 0 : " << tmp << std::endl;
1791 om_->print(
Debug,os.str());
1816 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1845 MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*tmpmvec_);
1846 MVT::MvTimesMatAddMv(ONE,*cXHP,*CP,ZERO,*R_);
1847 cXHP = Teuchos::null;
1848 MVT::SetBlock(*tmpmvec_,indblock1,*V_);
1849 MVT::SetBlock(*R_ ,indblock3,*V_);
1851 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*tmpmvec_);
1852 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CP,ZERO,*R_);
1853 cK_XHP = Teuchos::null;
1854 MVT::SetBlock(*tmpmvec_,indblock1,*KV_);
1855 MVT::SetBlock(*R_ ,indblock3,*KV_);
1858 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*tmpmvec_);
1859 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CP,ZERO,*R_);
1860 cM_XHP = Teuchos::null;
1861 MVT::SetBlock(*tmpmvec_,indblock1,*MV_);
1862 MVT::SetBlock(*R_ ,indblock3,*MV_);
1865 cM_XHP = Teuchos::null;
1870 MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*R_);
1871 cXHP = Teuchos::null;
1872 MVT::SetBlock(*R_,indblock1,*V_);
1873 MVT::MvTimesMatAddMv(ONE,*cHP,*CP,ZERO,*R_);
1874 cHP = Teuchos::null;
1875 MVT::SetBlock(*R_,indblock3,*V_);
1877 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*R_);
1878 cK_XHP = Teuchos::null;
1879 MVT::SetBlock(*R_,indblock1,*KV_);
1880 MVT::MvTimesMatAddMv(ONE,*cK_HP,*CP,ZERO,*R_);
1881 cK_HP = Teuchos::null;
1882 MVT::SetBlock(*R_,indblock3,*KV_);
1885 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*R_);
1886 cM_XHP = Teuchos::null;
1887 MVT::SetBlock(*R_,indblock1,*MV_);
1888 MVT::MvTimesMatAddMv(ONE,*cM_HP,*CP,ZERO,*R_);
1889 cM_HP = Teuchos::null;
1890 MVT::SetBlock(*R_,indblock3,*MV_);
1893 cM_XHP = Teuchos::null;
1894 cM_HP = Teuchos::null;
1909 || cHP != Teuchos::null || cK_HP != Teuchos::null || cM_HP != Teuchos::null
1910 || cP != Teuchos::null || cK_P != Teuchos::null || cM_P != Teuchos::null,
1912 "Anasazi::BlockKrylovSchur::iterate(): const views were not all cleared! Something went wrong!" );
1921 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1924 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1926 for (
int i = 0; i < blockSize_; i++) {
1929 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1933 Rnorms_current_ =
false;
1934 R2norms_current_ =
false;
1937 if (om_->isVerbosity(
Debug ) ) {
1947 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1954 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );
1962 template <
class ScalarType,
class MV,
class OP>
1963 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1965 if (Rnorms_current_ ==
false) {
1967 orthman_->norm(*R_,Rnorms_);
1968 Rnorms_current_ =
true;
1976 template <
class ScalarType,
class MV,
class OP>
1977 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1979 if (R2norms_current_ ==
false) {
1981 MVT::MvNorm(*R_,R2norms_);
1982 R2norms_current_ =
true;
2017 template <
class ScalarType,
class MV,
class OP>
2022 std::stringstream os;
2024 os.setf(std::ios::scientific, std::ios::floatfield);
2027 os <<
" Debugging checks: iteration " << iter_ << where << endl;
2030 if (chk.checkX && initialized_) {
2031 tmp = orthman_->orthonormError(*X_);
2032 os <<
" >> Error in X^H M X == I : " << tmp << endl;
2033 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2034 tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
2035 os <<
" >> Error in X^H M Q[" << i <<
"] == 0 : " << tmp << endl;
2038 if (chk.checkMX && hasM_ && initialized_) {
2039 tmp = Utils::errorEquality(*X_, *MX_, MOp_);
2040 os <<
" >> Error in MX == M*X : " << tmp << endl;
2042 if (chk.checkKX && initialized_) {
2043 tmp = Utils::errorEquality(*X_, *KX_, Op_);
2044 os <<
" >> Error in KX == K*X : " << tmp << endl;
2048 if (chk.checkP && hasP_ && initialized_) {
2050 tmp = orthman_->orthonormError(*P_);
2051 os <<
" >> Error in P^H M P == I : " << tmp << endl;
2052 tmp = orthman_->orthogError(*P_,*X_);
2053 os <<
" >> Error in P^H M X == 0 : " << tmp << endl;
2055 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2056 tmp = orthman_->orthogError(*P_,*auxVecs_[i]);
2057 os <<
" >> Error in P^H M Q[" << i <<
"] == 0 : " << tmp << endl;
2060 if (chk.checkMP && hasM_ && hasP_ && initialized_) {
2061 tmp = Utils::errorEquality(*P_, *MP_, MOp_);
2062 os <<
" >> Error in MP == M*P : " << tmp << endl;
2064 if (chk.checkKP && hasP_ && initialized_) {
2065 tmp = Utils::errorEquality(*P_, *KP_, Op_);
2066 os <<
" >> Error in KP == K*P : " << tmp << endl;
2070 if (chk.checkH && initialized_) {
2072 tmp = orthman_->orthonormError(*H_);
2073 os <<
" >> Error in H^H M H == I : " << tmp << endl;
2074 tmp = orthman_->orthogError(*H_,*X_);
2075 os <<
" >> Error in H^H M X == 0 : " << tmp << endl;
2077 tmp = orthman_->orthogError(*H_,*P_);
2078 os <<
" >> Error in H^H M P == 0 : " << tmp << endl;
2081 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2082 tmp = orthman_->orthogError(*H_,*auxVecs_[i]);
2083 os <<
" >> Error in H^H M Q[" << i <<
"] == 0 : " << tmp << endl;
2086 if (chk.checkMH && hasM_ && initialized_) {
2087 tmp = Utils::errorEquality(*H_, *MH_, MOp_);
2088 os <<
" >> Error in MH == M*H : " << tmp << endl;
2092 if (chk.checkR && initialized_) {
2094 MVT::MvTransMv(ONE,*X_,*R_,xTx);
2095 tmp = xTx.normFrobenius();
2096 MVT::MvTransMv(ONE,*R_,*R_,xTx);
2097 double normR = xTx.normFrobenius();
2098 os <<
" >> RelError in X^H R == 0: " << tmp/normR << endl;
2103 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2104 tmp = orthman_->orthonormError(*auxVecs_[i]);
2105 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << i <<
"] == I : " << tmp << endl;
2106 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
2107 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
2108 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << j <<
"] == 0 : " << tmp << endl;
2121 template <
class ScalarType,
class MV,
class OP>
2127 os.setf(std::ios::scientific, std::ios::floatfield);
2130 os <<
"================================================================================" << endl;
2132 os <<
" LOBPCG Solver Status" << endl;
2134 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
2135 os <<
"The number of iterations performed is " << iter_ << endl;
2136 os <<
"The current block size is " << blockSize_ << endl;
2137 os <<
"The number of auxiliary vectors is " << numAuxVecs_ << endl;
2138 os <<
"The number of operations Op*x is " << count_ApplyOp_ << endl;
2139 os <<
"The number of operations M*x is " << count_ApplyM_ << endl;
2140 os <<
"The number of operations Prec*x is " << count_ApplyPrec_ << endl;
2142 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2146 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
2147 os << std::setw(20) <<
"Eigenvalue" 2148 << std::setw(20) <<
"Residual(M)" 2149 << std::setw(20) <<
"Residual(2)" 2151 os <<
"--------------------------------------------------------------------------------"<<endl;
2152 for (
int i=0; i<blockSize_; i++) {
2153 os << std::setw(20) << theta_[i];
2154 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
2155 else os << std::setw(20) <<
"not current";
2156 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
2157 else os << std::setw(20) <<
"not current";
2161 os <<
"================================================================================" << endl;
2167 template <
class ScalarType,
class MV,
class OP>
2169 return initialized_;
2175 template <
class ScalarType,
class MV,
class OP>
2182 template <
class ScalarType,
class MV,
class OP>
2190 template <
class ScalarType,
class MV,
class OP>
2197 template <
class ScalarType,
class MV,
class OP>
2204 template <
class ScalarType,
class MV,
class OP>
2212 template <
class ScalarType,
class MV,
class OP>
2214 return 3*blockSize_;
2219 template <
class ScalarType,
class MV,
class OP>
2221 if (!initialized_)
return 0;
2228 template <
class ScalarType,
class MV,
class OP>
2229 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2232 return this->getRes2Norms();
2238 template <
class ScalarType,
class MV,
class OP>
2240 std::vector<int> ret(nevLocal_,0);
2247 template <
class ScalarType,
class MV,
class OP>
2249 std::vector<Value<ScalarType> > ret(nevLocal_);
2250 for (
int i=0; i<nevLocal_; i++) {
2251 ret[i].realpart = theta_[i];
2252 ret[i].imagpart = ZERO;
2259 template <
class ScalarType,
class MV,
class OP>
2262 "Anasazi::LOBPCG::setStatusTest() was passed a null StatusTest.");
2268 template <
class ScalarType,
class MV,
class OP>
2275 template <
class ScalarType,
class MV,
class OP>
2283 template <
class ScalarType,
class MV,
class OP>
2291 template <
class ScalarType,
class MV,
class OP>
2299 template <
class ScalarType,
class MV,
class OP>
2311 state.
T =
Teuchos::rcp(
new std::vector<MagnitudeType>(theta_));
2319 state.
MX = Teuchos::null;
2320 state.
MP = Teuchos::null;
2321 state.
MH = Teuchos::null;
2328 #endif // ANASAZI_LOBPCG_HPP ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
ScalarType * values() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
LOBPCG(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< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
LOBPCG constructor with eigenproblem, solver utilities, and parameter list of solver options...
This class defines the interface required by an eigensolver and status test class to compute solution...
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
Declaration of basic traits for the multivector type.
T & get(const std::string &name, T def_value)
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb) const
LOBPCGState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration...
Teuchos::RCP< const MultiVector > V
The current test basis.
virtual ~LOBPCG()
LOBPCG destructor
Virtual base class which defines basic traits for the operator type.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For LOBPCG, this always returns 3*getBlo...
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
Teuchos::RCP< const MultiVector > P
The current search direction.
An exception class parent to all Anasazi exceptions.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
Teuchos::RCP< const MultiVector > MV
The image of the current test basis under M, or Teuchos::null if M was not specified.
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...
Teuchos::RCP< const MultiVector > R
The current residual vectors.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Anasazi's templated, static class providing utilities for the solvers.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
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...
Traits class which defines basic operations on multivectors.
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.
Teuchos::RCP< const MultiVector > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Teuchos::RCP< const MultiVector > KX
The image of the current eigenvectors under K.
bool getFullOrtho() const
Determine if the LOBPCG iteration is using full orthogonality.
Teuchos::RCP< const MultiVector > KH
The image of the current preconditioned residual vectors under K.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Teuchos::RCP< const MultiVector > MH
The image of the current preconditioned residual vectors under M, or Teuchos::null if M was not speci...
Teuchos::RCP< const MultiVector > KP
The image of the current search direction under K.
ScalarTraits< ScalarType >::magnitudeType normOne() const
void push_back(const value_type &x)
int getNumIters() const
Get the current iteration count.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
void iterate()
This method performs LOBPCG iterations until the status test indicates the need to stop or an error o...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
LOBPCGRitzFailure is thrown when the LOBPCG solver is unable to continue a call to LOBPCG::iterate() ...
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
Teuchos::RCP< const MultiVector > KV
The image of the current test basis under K.
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 resetNumIters()
Reset the iteration count.
OrdinalType stride() const
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
LOBPCGOrthoFailure is thrown when an orthogonalization attempt fails.
LOBPCGInitFailure is thrown when the LOBPCG solver is unable to generate an initial iterate in the LO...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void setFullOrtho(bool fullOrtho)
Instruct the LOBPCG iteration to use full orthogonality.
Structure to contain pointers to Anasazi state variables.
Teuchos::RCP< const MultiVector > X
The current eigenvectors.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
bool hasP()
Indicates whether the search direction given by getState() is valid.
void POTRF(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Class which provides internal utilities for the Anasazi solvers.