235 LOBPCG(
const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
236 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
237 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
238 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
239 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
240 Teuchos::ParameterList ¶ms
304 void initialize(LOBPCGState<ScalarType,MV>& newstate);
342 LOBPCGState<ScalarType,MV>
getState()
const;
386 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
394 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
404 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
428 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
431 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
434 const Eigenproblem<ScalarType,MV,OP>&
getProblem()
const;
463 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
466 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const;
504 typedef SolverUtils<ScalarType,MV,OP> Utils;
505 typedef MultiVecTraits<ScalarType,MV> MVT;
506 typedef OperatorTraits<ScalarType,MV,OP> OPT;
507 typedef Teuchos::ScalarTraits<ScalarType> SCT;
508 typedef typename SCT::magnitudeType MagnitudeType;
509 const MagnitudeType ONE;
510 const MagnitudeType ZERO;
511 const MagnitudeType NANVAL;
516 bool checkX, checkMX, checkKX;
517 bool checkH, checkMH;
518 bool checkP, checkMP, checkKP;
520 CheckList() : checkX(false),checkMX(false),checkKX(false),
521 checkH(false),checkMH(false),
522 checkP(false),checkMP(false),checkKP(false),
523 checkR(false),checkQ(false) {};
528 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
532 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
533 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
534 const Teuchos::RCP<OutputManager<ScalarType> > om_;
535 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
536 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
540 Teuchos::RCP<const OP> Op_;
541 Teuchos::RCP<const OP> MOp_;
542 Teuchos::RCP<const OP> Prec_;
547 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
549 timerLocalProj_, timerDS_,
550 timerLocalUpdate_, timerCompRes_,
551 timerOrtho_, timerInit_;
556 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
585 Teuchos::RCP<MV> V_, KV_, MV_, R_;
587 Teuchos::RCP<MV> X_, KX_, MX_,
601 Teuchos::RCP<MV> tmpmvec_;
604 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
611 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
614 bool Rnorms_current_, R2norms_current_;
979#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
980 Teuchos::TimeMonitor inittimer( *timerInit_ );
983 std::vector<int> bsind(blockSize_);
984 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
1008 if (newstate.X != Teuchos::null) {
1009 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1010 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.X not correct." );
1012 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
1013 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.X must have at least block size vectors.");
1016 MVT::SetBlock(*newstate.X,bsind,*X_);
1020 if (newstate.MX != Teuchos::null) {
1021 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1022 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MX not correct." );
1024 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MX) < blockSize_,
1025 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.MX must have at least block size vectors.");
1026 MVT::SetBlock(*newstate.MX,bsind,*MX_);
1031#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1032 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1034 OPT::Apply(*MOp_,*X_,*MX_);
1035 count_ApplyM_ += blockSize_;
1038 newstate.R = Teuchos::null;
1043 if (newstate.KX != Teuchos::null) {
1044 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1045 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KX not correct." );
1047 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KX) < blockSize_,
1048 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.KX must have at least block size vectors.");
1049 MVT::SetBlock(*newstate.KX,bsind,*KX_);
1054#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1055 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1057 OPT::Apply(*Op_,*X_,*KX_);
1058 count_ApplyOp_ += blockSize_;
1061 newstate.R = Teuchos::null;
1069 newstate.P = Teuchos::null;
1070 newstate.KP = Teuchos::null;
1071 newstate.MP = Teuchos::null;
1072 newstate.R = Teuchos::null;
1073 newstate.T = Teuchos::null;
1076 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1077 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1078 "Anasazi::LOBPCG::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1080 int initSize = MVT::GetNumberVecs(*ivec);
1081 if (initSize > blockSize_) {
1083 initSize = blockSize_;
1084 std::vector<int> ind(blockSize_);
1085 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1086 ivec = MVT::CloneView(*ivec,ind);
1091 std::vector<int> ind(initSize);
1092 for (
int i=0; i<initSize; i++) ind[i] = i;
1093 MVT::SetBlock(*ivec,ind,*X_);
1096 if (blockSize_ > initSize) {
1097 std::vector<int> ind(blockSize_ - initSize);
1098 for (
int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1099 Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X_,ind);
1106#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1107 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1109 OPT::Apply(*MOp_,*X_,*MX_);
1110 count_ApplyM_ += blockSize_;
1114 if (numAuxVecs_ > 0) {
1115#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1116 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1118 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
1119 int rank = orthman_->projectAndNormalizeMat(*X_,auxVecs_,dummy,Teuchos::null,MX_);
1121 "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1124#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1125 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1127 int rank = orthman_->normalizeMat(*X_,Teuchos::null,MX_);
1129 "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1134#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1135 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1137 OPT::Apply(*Op_,*X_,*KX_);
1138 count_ApplyOp_ += blockSize_;
1146 theta_.resize(3*blockSize_,NANVAL);
1147 if (newstate.T != Teuchos::null) {
1148 TEUCHOS_TEST_FOR_EXCEPTION( (
signed int)(newstate.T->size()) < blockSize_,
1149 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1150 for (
int i=0; i<blockSize_; i++) {
1151 theta_[i] = (*newstate.T)[i];
1153 nevLocal_ = blockSize_;
1157 Teuchos::SerialDenseMatrix<int,ScalarType> KK(blockSize_,blockSize_),
1158 MM(blockSize_,blockSize_),
1159 S(blockSize_,blockSize_);
1161#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1162 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1165 MVT::MvTransMv(ONE,*X_,*KX_,KK);
1167 MVT::MvTransMv(ONE,*X_,*MX_,MM);
1168 nevLocal_ = blockSize_;
1173#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1174 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1176 Utils::directSolver(blockSize_, KK, Teuchos::rcpFromRef(MM), S, theta_, nevLocal_, 1);
1178 "Anasazi::LOBPCG::initialize(): Initial Ritz analysis did not produce enough Ritz pairs to initialize algorithm.");
1184#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1185 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1188 std::vector<int> order(blockSize_);
1191 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1194 Utils::permuteVectors(order,S);
1199#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1200 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1203 MVT::MvAddMv( ONE, *X_, ZERO, *X_, *R_ );
1204 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X_ );
1206 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1207 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *KX_ );
1210 MVT::MvAddMv( ONE, *MX_, ZERO, *MX_, *R_ );
1211 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *MX_ );
1219 if (newstate.R != Teuchos::null) {
1220 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
1221 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.R not correct." );
1222 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
1223 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.R must have blockSize number of vectors." );
1224 MVT::SetBlock(*newstate.R,bsind,*R_);
1227#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1228 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1231 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1232 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1233 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1234 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1238 Rnorms_current_ =
false;
1239 R2norms_current_ =
false;
1242 if (newstate.P != Teuchos::null) {
1243 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.P) < blockSize_ ,
1244 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.P must have blockSize number of vectors." );
1245 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.P) != MVT::GetGlobalLength(*P_),
1246 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.P not correct." );
1250 MVT::SetBlock(*newstate.P,bsind,*P_);
1253 if (newstate.KP != Teuchos::null) {
1254 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KP) < blockSize_,
1255 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.KP must have blockSize number of vectors." );
1256 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KP) != MVT::GetGlobalLength(*KP_),
1257 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KP not correct." );
1258 MVT::SetBlock(*newstate.KP,bsind,*KP_);
1261#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1262 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1264 OPT::Apply(*Op_,*P_,*KP_);
1265 count_ApplyOp_ += blockSize_;
1270 if (newstate.MP != Teuchos::null) {
1271 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MP) < blockSize_,
1272 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): newstate.MP must have blockSize number of vectors." );
1273 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MP) != MVT::GetGlobalLength(*MP_),
1274 std::invalid_argument,
"Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MP not correct." );
1275 MVT::SetBlock(*newstate.MP,bsind,*MP_);
1278#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1279 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1281 OPT::Apply(*MOp_,*P_,*MP_);
1282 count_ApplyM_ += blockSize_;
1291 initialized_ =
true;
1293 if (om_->isVerbosity(
Debug ) ) {
1304 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
1355 if (initialized_ ==
false) {
1361 const int oneBlock = blockSize_;
1362 const int twoBlocks = 2*blockSize_;
1363 const int threeBlocks = 3*blockSize_;
1365 std::vector<int> indblock1(blockSize_), indblock2(blockSize_), indblock3(blockSize_);
1366 for (
int i=0; i<blockSize_; i++) {
1368 indblock2[i] = i + blockSize_;
1369 indblock3[i] = i + 2*blockSize_;
1377 Teuchos::SerialDenseMatrix<int,ScalarType> KK( threeBlocks, threeBlocks ),
1378 MM( threeBlocks, threeBlocks ),
1379 S( threeBlocks, threeBlocks );
1381 while (tester_->checkStatus(
this) !=
Passed) {
1384 if (om_->isVerbosity(
Debug)) {
1385 currentStatus( om_->stream(
Debug) );
1395 if (Prec_ != Teuchos::null) {
1396#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1397 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1399 OPT::Apply( *Prec_, *R_, *H_ );
1400 count_ApplyPrec_ += blockSize_;
1403 MVT::MvAddMv(ONE,*R_,ZERO,*R_,*H_);
1408#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1409 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1411 OPT::Apply( *MOp_, *H_, *MH_);
1412 count_ApplyM_ += blockSize_;
1417 Teuchos::Array<Teuchos::RCP<const MV> > Q;
1428#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1429 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1431 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC =
1432 Teuchos::tuple<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >(Teuchos::null);
1433 int rank = orthman_->projectAndNormalizeMat(*H_,Q,dummyC,Teuchos::null,MH_);
1436 "Anasazi::LOBPCG::iterate(): unable to compute orthonormal basis for H.");
1439 if (om_->isVerbosity(
Debug ) ) {
1443 om_->print(
Debug, accuracyCheck(chk,
": after ortho H") );
1449 om_->print(
OrthoDetails, accuracyCheck(chk,
": after ortho H") );
1454#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1455 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1457 OPT::Apply( *Op_, *H_, *KH_);
1458 count_ApplyOp_ += blockSize_;
1462 nevLocal_ = threeBlocks;
1465 nevLocal_ = twoBlocks;
1487 KX_ = Teuchos::null;
1488 MX_ = Teuchos::null;
1490 KH_ = Teuchos::null;
1491 MH_ = Teuchos::null;
1493 KP_ = Teuchos::null;
1494 MP_ = Teuchos::null;
1495 Teuchos::RCP<const MV> cX, cH, cXHP, cHP, cK_XHP, cK_HP, cM_XHP, cM_HP, cP, cK_P, cM_P;
1497 cX = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock1);
1498 cH = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock2);
1500 std::vector<int> indXHP(nevLocal_);
1501 for (
int i=0; i<nevLocal_; i++) {
1504 cXHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indXHP);
1505 cK_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indXHP);
1507 cM_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indXHP);
1513 std::vector<int> indHP(nevLocal_-blockSize_);
1514 for (
int i=blockSize_; i<nevLocal_; i++) {
1515 indHP[i-blockSize_] = i;
1517 cHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indHP);
1518 cK_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indHP);
1520 cM_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indHP);
1526 if (nevLocal_ == threeBlocks) {
1527 cP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock3);
1528 cK_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indblock3);
1530 cM_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indblock3);
1552 Teuchos::SerialDenseMatrix<int,ScalarType>
1553 K1(Teuchos::View,KK,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1554 K2(Teuchos::View,KK,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1555 K3(Teuchos::View,KK,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_),
1556 M1(Teuchos::View,MM,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1557 M2(Teuchos::View,MM,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1558 M3(Teuchos::View,MM,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_);
1560#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1561 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1563 MVT::MvTransMv( ONE, *cX, *cK_XHP, K1 );
1564 MVT::MvTransMv( ONE, *cX, *cM_XHP, M1 );
1565 MVT::MvTransMv( ONE, *cH, *cK_HP , K2 );
1566 MVT::MvTransMv( ONE, *cH, *cM_HP , M2 );
1567 if (nevLocal_ == threeBlocks) {
1568 MVT::MvTransMv( ONE, *cP, *cK_P, K3 );
1569 MVT::MvTransMv( ONE, *cP, *cM_P, M3 );
1579 cK_P = Teuchos::null;
1580 cM_P = Teuchos::null;
1581 if (fullOrtho_ ==
true) {
1582 cHP = Teuchos::null;
1583 cK_HP = Teuchos::null;
1584 cM_HP = Teuchos::null;
1593 Teuchos::SerialDenseMatrix<int,ScalarType> lclKK(Teuchos::View,KK,nevLocal_,nevLocal_),
1594 lclMM(Teuchos::View,MM,nevLocal_,nevLocal_),
1595 lclS(Teuchos::View, S,nevLocal_,nevLocal_);
1597#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1598 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1600 int localSize = nevLocal_;
1601 Utils::directSolver(localSize, lclKK, Teuchos::rcpFromRef(lclMM), lclS, theta_, nevLocal_, 0);
1610 if (nevLocal_ != localSize) {
1613 cXHP = Teuchos::null;
1614 cK_XHP = Teuchos::null;
1615 cM_XHP = Teuchos::null;
1616 cHP = Teuchos::null;
1617 cK_HP = Teuchos::null;
1618 cM_HP = Teuchos::null;
1622 "Anasazi::LOBPCG::iterate(): indefiniteness detected in projected mass matrix." );
1629 Teuchos::LAPACK<int,ScalarType> lapack;
1630 Teuchos::BLAS<int,ScalarType> blas;
1632#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1633 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1636 std::vector<int> order(nevLocal_);
1639 sm_->sort(theta_, Teuchos::rcpFromRef(order), nevLocal_);
1642 Utils::permuteVectors(order,lclS);
1654 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > CX, CP;
1665 Teuchos::SerialDenseMatrix<int,ScalarType> C(nevLocal_,twoBlocks),
1666 MMC(nevLocal_,twoBlocks),
1667 CMMC(twoBlocks ,twoBlocks);
1670 for (
int j=0; j<oneBlock; j++) {
1671 for (
int i=0; i<oneBlock; i++) {
1675 C(i,j+oneBlock) = ZERO;
1679 for (
int j=0; j<oneBlock; j++) {
1680 for (
int i=oneBlock; i<twoBlocks; i++) {
1684 C(i,j+oneBlock) = lclS(i,j);
1688 if (nevLocal_ == threeBlocks) {
1689 for (
int j=0; j<oneBlock; j++) {
1690 for (
int i=twoBlocks; i<threeBlocks; i++) {
1694 C(i,j+oneBlock) = lclS(i,j);
1702 teuchosret = MMC.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,C,ZERO);
1703 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1704 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1706 teuchosret = CMMC.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,C,MMC,ZERO);
1707 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1708 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1714 lapack.POTRF(
'U',twoBlocks,CMMC.values(),CMMC.stride(),&info);
1718 Teuchos::SerialDenseMatrix<int,ScalarType> K22(Teuchos::View,KK,blockSize_,blockSize_,1*blockSize_,1*blockSize_);
1719 Teuchos::SerialDenseMatrix<int,ScalarType> K22_trans( K22, Teuchos::CONJ_TRANS );
1721 MagnitudeType symNorm = K22_trans.normOne();
1722 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
1724 cXHP = Teuchos::null;
1725 cHP = Teuchos::null;
1726 cK_XHP = Teuchos::null;
1727 cK_HP = Teuchos::null;
1728 cM_XHP = Teuchos::null;
1729 cM_HP = Teuchos::null;
1732 std::string errMsg =
"Anasazi::LOBPCG::iterate(): Cholesky factorization failed during full orthogonalization.";
1733 if ( symNorm < tol )
1739 errMsg +=
" Check the operator A (or K), it may not be Hermitian!";
1745 blas.TRSM(Teuchos::RIGHT_SIDE,Teuchos::UPPER_TRI,Teuchos::NO_TRANS,Teuchos::NON_UNIT_DIAG,
1746 nevLocal_,twoBlocks,ONE,CMMC.values(),CMMC.stride(),C.values(),C.stride());
1748 CX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,0) );
1750 CP = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,oneBlock) );
1753 if (om_->isVerbosity(
Debug ) ) {
1754 Teuchos::SerialDenseMatrix<int,ScalarType> tmp1(nevLocal_,oneBlock),
1755 tmp2(oneBlock,oneBlock);
1758 std::stringstream os;
1760 os.setf(std::ios::scientific, std::ios::floatfield);
1762 os <<
" Checking Full Ortho: iteration " << iter_ << std::endl;
1766 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CX,ZERO);
1767 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1768 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1770 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1771 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1772 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1774 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1775 tmp = tmp2.normFrobenius();
1776 os <<
" >> Error in CX^H MM CX == I : " << tmp << std::endl;
1780 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1781 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1782 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1784 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CP,tmp1,ZERO);
1785 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1786 "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1788 for (
int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1789 tmp = tmp2.normFrobenius();
1790 os <<
" >> Error in CP^H MM CP == I : " << tmp << std::endl;
1794 teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1795 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1797 teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1798 TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1800 tmp = tmp2.normFrobenius();
1801 os <<
" >> Error in CX^H MM CP == 0 : " << tmp << std::endl;
1804 om_->print(
Debug,os.str());
1819 CX = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_ ,oneBlock,0 ,0) );
1820 CP = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_-oneBlock,oneBlock,oneBlock,0) );
1829#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1830 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1858 MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*tmpmvec_);
1859 MVT::MvTimesMatAddMv(ONE,*cXHP,*CP,ZERO,*R_);
1860 cXHP = Teuchos::null;
1861 MVT::SetBlock(*tmpmvec_,indblock1,*V_);
1862 MVT::SetBlock(*R_ ,indblock3,*V_);
1864 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*tmpmvec_);
1865 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CP,ZERO,*R_);
1866 cK_XHP = Teuchos::null;
1867 MVT::SetBlock(*tmpmvec_,indblock1,*KV_);
1868 MVT::SetBlock(*R_ ,indblock3,*KV_);
1871 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*tmpmvec_);
1872 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CP,ZERO,*R_);
1873 cM_XHP = Teuchos::null;
1874 MVT::SetBlock(*tmpmvec_,indblock1,*MV_);
1875 MVT::SetBlock(*R_ ,indblock3,*MV_);
1878 cM_XHP = Teuchos::null;
1883 MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*R_);
1884 cXHP = Teuchos::null;
1885 MVT::SetBlock(*R_,indblock1,*V_);
1886 MVT::MvTimesMatAddMv(ONE,*cHP,*CP,ZERO,*R_);
1887 cHP = Teuchos::null;
1888 MVT::SetBlock(*R_,indblock3,*V_);
1890 MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*R_);
1891 cK_XHP = Teuchos::null;
1892 MVT::SetBlock(*R_,indblock1,*KV_);
1893 MVT::MvTimesMatAddMv(ONE,*cK_HP,*CP,ZERO,*R_);
1894 cK_HP = Teuchos::null;
1895 MVT::SetBlock(*R_,indblock3,*KV_);
1898 MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*R_);
1899 cM_XHP = Teuchos::null;
1900 MVT::SetBlock(*R_,indblock1,*MV_);
1901 MVT::MvTimesMatAddMv(ONE,*cM_HP,*CP,ZERO,*R_);
1902 cM_HP = Teuchos::null;
1903 MVT::SetBlock(*R_,indblock3,*MV_);
1906 cM_XHP = Teuchos::null;
1907 cM_HP = Teuchos::null;
1921 TEUCHOS_TEST_FOR_EXCEPTION( cXHP != Teuchos::null || cK_XHP != Teuchos::null || cM_XHP != Teuchos::null
1922 || cHP != Teuchos::null || cK_HP != Teuchos::null || cM_HP != Teuchos::null
1923 || cP != Teuchos::null || cK_P != Teuchos::null || cM_P != Teuchos::null,
1925 "Anasazi::BlockKrylovSchur::iterate(): const views were not all cleared! Something went wrong!" );
1934#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1935 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1937 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1938 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1939 for (
int i = 0; i < blockSize_; i++) {
1942 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1946 Rnorms_current_ =
false;
1947 R2norms_current_ =
false;
1950 if (om_->isVerbosity(
Debug ) ) {
1960 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1967 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );