52 #ifndef ANASAZI_LOBPCG_HPP 53 #define ANASAZI_LOBPCG_HPP 60 #include "Teuchos_ScalarTraits.hpp" 65 #include "Teuchos_LAPACK.hpp" 66 #include "Teuchos_BLAS.hpp" 67 #include "Teuchos_SerialDenseMatrix.hpp" 68 #include "Teuchos_ParameterList.hpp" 69 #include "Teuchos_TimeMonitor.hpp" 100 template <
class ScalarType,
class MultiVector>
103 Teuchos::RCP<const MultiVector>
V;
105 Teuchos::RCP<const MultiVector>
KV;
107 Teuchos::RCP<const MultiVector>
MV;
110 Teuchos::RCP<const MultiVector>
X;
112 Teuchos::RCP<const MultiVector>
KX;
114 Teuchos::RCP<const MultiVector>
MX;
117 Teuchos::RCP<const MultiVector>
P;
119 Teuchos::RCP<const MultiVector>
KP;
121 Teuchos::RCP<const MultiVector>
MP;
127 Teuchos::RCP<const MultiVector>
H;
129 Teuchos::RCP<const MultiVector>
KH;
131 Teuchos::RCP<const MultiVector>
MH;
134 Teuchos::RCP<const MultiVector>
R;
137 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >
T;
140 V(Teuchos::null),KV(Teuchos::null),MV(Teuchos::null),
141 X(Teuchos::null),KX(Teuchos::null),MX(Teuchos::null),
142 P(Teuchos::null),KP(Teuchos::null),MP(Teuchos::null),
143 H(Teuchos::null),KH(Teuchos::null),MH(Teuchos::null),
144 R(Teuchos::null),T(Teuchos::null) {};
208 template <
class ScalarType,
class MV,
class OP>
223 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
227 Teuchos::ParameterList ¶ms
317 bool isInitialized()
const;
337 int getNumIters()
const;
340 void resetNumIters();
349 Teuchos::RCP<const MV> getRitzVectors();
356 std::vector<Value<ScalarType> > getRitzValues();
365 std::vector<int> getRitzIndex();
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();
402 int getCurSubspaceDim()
const;
407 int getMaxSubspaceDim()
const;
418 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest()
const;
432 void setBlockSize(
int blockSize);
436 int getBlockSize()
const;
450 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
453 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs()
const;
465 void setFullOrtho(
bool fullOrtho);
468 bool getFullOrtho()
const;
479 void currentStatus(std::ostream &os);
494 typedef Teuchos::ScalarTraits<ScalarType> SCT;
495 typedef typename SCT::magnitudeType MagnitudeType;
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;
519 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
520 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
521 const Teuchos::RCP<OutputManager<ScalarType> > om_;
522 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
523 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
527 Teuchos::RCP<const OP> Op_;
528 Teuchos::RCP<const OP> MOp_;
529 Teuchos::RCP<const OP> Prec_;
534 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
536 timerLocalProj_, timerDS_,
537 timerLocalUpdate_, timerCompRes_,
538 timerOrtho_, timerInit_;
543 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
572 Teuchos::RCP<MV> V_, KV_, MV_, R_;
574 Teuchos::RCP<MV> X_, KX_, MX_,
588 Teuchos::RCP<MV> tmpmvec_;
591 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
598 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
601 bool Rnorms_current_, R2norms_current_;
610 template <
class ScalarType,
class MV,
class OP>
613 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
617 Teuchos::ParameterList ¶ms
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)),
650 auxVecs_( Teuchos::Array<Teuchos::RCP<const
MV> >(0) ),
653 Rnorms_current_(false),
654 R2norms_current_(false)
656 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
657 "Anasazi::LOBPCG::constructor: user passed null problem pointer.");
658 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
659 "Anasazi::LOBPCG::constructor: user passed null sort manager pointer.");
660 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
661 "Anasazi::LOBPCG::constructor: user passed null output manager pointer.");
662 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
663 "Anasazi::LOBPCG::constructor: user passed null status test pointer.");
664 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
665 "Anasazi::LOBPCG::constructor: user passed null orthogonalization manager pointer.");
666 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
667 "Anasazi::LOBPCG::constructor: problem is not set.");
668 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
669 "Anasazi::LOBPCG::constructor: problem is not Hermitian; LOBPCG requires Hermitian problem.");
672 Op_ = problem_->getOperator();
673 TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
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 692 Teuchos::TimeMonitor lcltimer( *timerInit_ );
699 Teuchos::RCP<const MV> tmp;
705 if (blockSize_ > 0) {
709 tmp = problem_->getInitVec();
710 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
711 "Anasazi::LOBPCG::setBlockSize(): eigenproblem did not specify initial vectors to clone from.");
714 TEUCHOS_TEST_FOR_EXCEPTION(newBS <= 0 || static_cast<ptrdiff_t>(newBS) >
MVT::GetGlobalLength(*tmp),
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++) {
742 Teuchos::RCP<MV> newV, newMV, newKV, newR;
743 Teuchos::RCP<const MV> src;
769 theta_.resize(3*newBS);
770 Rnorms_.resize(newBS);
771 R2norms_.resize(newBS);
790 for (
int i=0; i<newBS; i++) {
791 newind[i] += 2*newBS;
792 oldind[i] += 2*blockSize_;
819 initialized_ =
false;
838 theta_.resize(3*newBS,NANVAL);
839 Rnorms_.resize(newBS,NANVAL);
840 R2norms_.resize(newBS,NANVAL);
855 tmpmvec_ = Teuchos::null;
870 template <
class ScalarType,
class MV,
class OP>
873 std::vector<int> ind(blockSize_);
875 for (
int i=0; i<blockSize_; i++) {
887 for (
int i=0; i<blockSize_; i++) {
888 ind[i] += blockSize_;
899 for (
int i=0; i<blockSize_; i++) {
900 ind[i] += blockSize_;
915 template <
class ScalarType,
class MV,
class OP>
917 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
923 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); 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 967 Teuchos::TimeMonitor inittimer( *timerInit_ );
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.");
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.");
1018 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1019 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
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.");
1041 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1042 Teuchos::TimeMonitor lcltimer( *timerOp_ );
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;
1063 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1064 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1065 "Anasazi::LOBPCG::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1068 if (initSize > blockSize_) {
1070 initSize = blockSize_;
1071 std::vector<int> ind(blockSize_);
1072 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1078 std::vector<int> ind(initSize);
1079 for (
int i=0; i<initSize; i++) ind[i] = i;
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 1094 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1097 count_ApplyM_ += blockSize_;
1101 if (numAuxVecs_ > 0) {
1102 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1103 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1105 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
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 1112 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
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 1122 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1125 count_ApplyOp_ += blockSize_;
1133 theta_.resize(3*blockSize_,NANVAL);
1134 if (newstate.
T != Teuchos::null) {
1135 TEUCHOS_TEST_FOR_EXCEPTION( (
signed int)(newstate.
T->size()) < blockSize_,
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_;
1144 Teuchos::SerialDenseMatrix<int,ScalarType> KK(blockSize_,blockSize_),
1145 MM(blockSize_,blockSize_),
1146 S(blockSize_,blockSize_);
1148 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1149 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1155 nevLocal_ = blockSize_;
1160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1161 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1165 "Anasazi::LOBPCG::initialize(): Initial Ritz analysis did not produce enough Ritz pairs to initialize algorithm.");
1171 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1172 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1175 std::vector<int> order(blockSize_);
1178 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1186 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1187 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
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." );
1214 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1215 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1219 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1220 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
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." );
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." );
1248 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1249 Teuchos::TimeMonitor lcltimer( *timerOp_ );
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." );
1265 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1266 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
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) {
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_;
1364 Teuchos::SerialDenseMatrix<int,ScalarType> KK( threeBlocks, threeBlocks ),
1365 MM( threeBlocks, threeBlocks ),
1366 S( threeBlocks, threeBlocks );
1368 while (tester_->checkStatus(
this) !=
Passed) {
1371 if (om_->isVerbosity(
Debug)) {
1382 if (Prec_ != Teuchos::null) {
1383 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1384 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1387 count_ApplyPrec_ += blockSize_;
1395 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1396 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1399 count_ApplyM_ += blockSize_;
1404 Teuchos::Array<Teuchos::RCP<const MV> > Q;
1415 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1416 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1418 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC =
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 1442 Teuchos::TimeMonitor lcltimer( *timerOp_ );
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);
1539 Teuchos::SerialDenseMatrix<int,ScalarType>
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 1548 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1554 if (nevLocal_ == threeBlocks) {
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;
1580 Teuchos::SerialDenseMatrix<int,ScalarType> lclKK(Teuchos::View,KK,nevLocal_,nevLocal_),
1581 lclMM(Teuchos::View,MM,nevLocal_,nevLocal_),
1582 lclS(Teuchos::View, S,nevLocal_,nevLocal_);
1584 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1585 Teuchos::TimeMonitor lcltimer( *timerDS_ );
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." );
1616 Teuchos::LAPACK<int,ScalarType> lapack;
1617 Teuchos::BLAS<int,ScalarType> blas;
1619 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1620 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1623 std::vector<int> order(nevLocal_);
1626 sm_->sort(theta_, Teuchos::rcpFromRef(order), nevLocal_);
1641 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > CX, CP;
1652 Teuchos::SerialDenseMatrix<int,ScalarType> C(nevLocal_,twoBlocks),
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);
1689 teuchosret = MMC.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,C,ZERO);
1690 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1691 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1693 teuchosret = CMMC.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,C,MMC,ZERO);
1694 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1695 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1701 lapack.POTRF(
'U',twoBlocks,CMMC.values(),CMMC.stride(),&info);
1705 Teuchos::SerialDenseMatrix<int,ScalarType> K22(Teuchos::View,KK,blockSize_,blockSize_,1*blockSize_,1*blockSize_);
1706 Teuchos::SerialDenseMatrix<int,ScalarType> K22_trans( K22, Teuchos::CONJ_TRANS );
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!";
1732 blas.TRSM(Teuchos::RIGHT_SIDE,Teuchos::UPPER_TRI,Teuchos::NO_TRANS,Teuchos::NON_UNIT_DIAG,
1733 nevLocal_,twoBlocks,ONE,CMMC.values(),CMMC.stride(),C.values(),C.stride());
1735 CX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,0) );
1737 CP = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,oneBlock) );
1740 if (om_->isVerbosity(
Debug ) ) {
1741 Teuchos::SerialDenseMatrix<int,ScalarType> tmp1(nevLocal_,oneBlock),
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;
1753 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CX,ZERO);
1754 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1755 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1757 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1758 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1759 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1761 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1762 tmp = tmp2.normFrobenius();
1763 os <<
" >> Error in CX^H MM CX == I : " << tmp << std::endl;
1767 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1768 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1769 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1771 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CP,tmp1,ZERO);
1772 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1773 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1775 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1776 tmp = tmp2.normFrobenius();
1777 os <<
" >> Error in CP^H MM CP == I : " << tmp << std::endl;
1781 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1782 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1784 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1785 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1787 tmp = tmp2.normFrobenius();
1788 os <<
" >> Error in CX^H MM CP == 0 : " << tmp << std::endl;
1791 om_->print(
Debug,os.str());
1806 CX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_ ,oneBlock,0 ,0) );
1807 CP = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_-oneBlock,oneBlock,oneBlock,0) );
1816 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1817 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1847 cXHP = Teuchos::null;
1853 cK_XHP = Teuchos::null;
1860 cM_XHP = Teuchos::null;
1865 cM_XHP = Teuchos::null;
1871 cXHP = Teuchos::null;
1874 cHP = Teuchos::null;
1878 cK_XHP = Teuchos::null;
1881 cK_HP = Teuchos::null;
1886 cM_XHP = Teuchos::null;
1889 cM_HP = Teuchos::null;
1893 cM_XHP = Teuchos::null;
1894 cM_HP = Teuchos::null;
1908 TEUCHOS_TEST_FOR_EXCEPTION( cXHP != Teuchos::null || cK_XHP != Teuchos::null || cM_XHP != 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 1922 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1925 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1926 for (
int i = 0; i < blockSize_; i++) {
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) {
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_) {
2040 os <<
" >> Error in MX == M*X : " << tmp << endl;
2042 if (chk.checkKX && initialized_) {
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_) {
2062 os <<
" >> Error in MP == M*P : " << tmp << endl;
2064 if (chk.checkKP && hasP_ && initialized_) {
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_) {
2088 os <<
" >> Error in MH == M*H : " << tmp << endl;
2092 if (chk.checkR && initialized_) {
2093 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
2095 tmp = xTx.normFrobenius();
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>
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>
2261 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
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 static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
Declaration of basic traits for the multivector type.
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.
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
Teuchos::RCP< const MultiVector > P
The current search direction.
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...
An exception class parent to all Anasazi exceptions.
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.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Anasazi's templated, static class providing utilities for the solvers.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
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.
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...
Traits class which defines basic operations on multivectors.
int getNumIters() const
Get the current iteration count.
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.
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.
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).
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().
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...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For LOBPCG, this always returns 3*getBlo...
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.
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 .
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
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...
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...
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.
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
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...
LOBPCGState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
bool hasP()
Indicates whether the search direction given by getState() is valid.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
Class which provides internal utilities for the Anasazi solvers.