176 BlockDavidson(
const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
177 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
178 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
179 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
180 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
181 Teuchos::ParameterList ¶ms
250 void initialize(BlockDavidsonState<ScalarType,MV>& newstate);
286 BlockDavidsonState<ScalarType,MV>
getState()
const;
332 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
340 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
350 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
369 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
372 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
375 const Eigenproblem<ScalarType,MV,OP>&
getProblem()
const;
403 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
406 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const;
422 void setSize(
int blockSize,
int numBlocks);
438 typedef SolverUtils<ScalarType,MV,OP> Utils;
439 typedef MultiVecTraits<ScalarType,MV> MVT;
440 typedef OperatorTraits<ScalarType,MV,OP> OPT;
441 typedef Teuchos::ScalarTraits<ScalarType> SCT;
442 typedef typename SCT::magnitudeType MagnitudeType;
443 const MagnitudeType ONE;
444 const MagnitudeType ZERO;
445 const MagnitudeType NANVAL;
451 bool checkX, checkMX, checkKX;
452 bool checkH, checkMH, checkKH;
455 CheckList() : checkV(
false),
456 checkX(
false),checkMX(
false),checkKX(
false),
457 checkH(
false),checkMH(
false),checkKH(
false),
458 checkR(
false),checkQ(
false),checkKK(
false) {};
463 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
467 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
468 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
469 const Teuchos::RCP<OutputManager<ScalarType> > om_;
470 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
471 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
475 Teuchos::RCP<const OP> Op_;
476 Teuchos::RCP<const OP> MOp_;
477 Teuchos::RCP<const OP> Prec_;
482 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
483 timerSortEval_, timerDS_,
484 timerLocal_, timerCompRes_,
485 timerOrtho_, timerInit_;
489 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
518 Teuchos::RCP<MV> X_, KX_, MX_, R_,
524 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
527 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
534 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
537 bool Rnorms_current_, R2norms_current_;
554 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
555 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
556 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
557 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
558 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
559 Teuchos::ParameterList ¶ms
561 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
562 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
563 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
571#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
572 timerOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation Op*x")),
573 timerMOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation M*x")),
574 timerPrec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation Prec*x")),
575 timerSortEval_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Sorting eigenvalues")),
576 timerDS_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Direct solve")),
577 timerLocal_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Local update")),
578 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Computing residuals")),
579 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Orthogonalization")),
580 timerInit_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Initialization")),
590 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
593 Rnorms_current_(false),
594 R2norms_current_(false)
596 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
597 "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
598 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
599 "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
600 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
601 "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
602 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
603 "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
604 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
605 "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
606 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
607 "Anasazi::BlockDavidson::constructor: problem is not set.");
608 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
609 "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
612 Op_ = problem_->getOperator();
613 TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
614 "Anasazi::BlockDavidson::constructor: problem provides no operator.");
615 MOp_ = problem_->getM();
616 Prec_ = problem_->getPrec();
617 hasM_ = (MOp_ != Teuchos::null);
620 int bs = params.get(
"Block Size", problem_->getNEV());
621 int nb = params.get(
"Num Blocks", 2);
903#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
904 Teuchos::TimeMonitor inittimer( *timerInit_ );
907 std::vector<int> bsind(blockSize_);
908 for (
int i=0; i<blockSize_; ++i) bsind[i] = i;
910 Teuchos::BLAS<int,ScalarType> blas;
934 Teuchos::RCP<MV> lclV;
935 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
937 if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
938 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
939 "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
940 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
941 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
942 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
943 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
944 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
945 "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
947 curDim_ = newstate.curDim;
949 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
951 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
952 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
955 TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
956 "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
959 std::vector<int> nevind(curDim_);
960 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
961 if (newstate.V != V_) {
962 if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
963 newstate.V = MVT::CloneView(*newstate.V,nevind);
965 MVT::SetBlock(*newstate.V,nevind,*V_);
967 lclV = MVT::CloneViewNonConst(*V_,nevind);
970 lclKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
971 if (newstate.KK != KK_) {
972 if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
973 newstate.KK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
975 lclKK->assign(*newstate.KK);
979 for (
int j=0; j<curDim_-1; ++j) {
980 for (
int i=j+1; i<curDim_; ++i) {
981 (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
988 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
989 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
990 "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
992 newstate.X = Teuchos::null;
993 newstate.MX = Teuchos::null;
994 newstate.KX = Teuchos::null;
995 newstate.R = Teuchos::null;
996 newstate.H = Teuchos::null;
997 newstate.T = Teuchos::null;
998 newstate.KK = Teuchos::null;
999 newstate.V = Teuchos::null;
1000 newstate.curDim = 0;
1002 curDim_ = MVT::GetNumberVecs(*ivec);
1004 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1005 if (curDim_ > blockSize_*numBlocks_) {
1008 curDim_ = blockSize_*numBlocks_;
1010 bool userand =
false;
1015 curDim_ = blockSize_;
1025 std::vector<int> dimind(curDim_);
1026 for (
int i=0; i<curDim_; ++i) dimind[i] = i;
1027 lclV = MVT::CloneViewNonConst(*V_,dimind);
1030 MVT::MvRandom(*lclV);
1033 if (MVT::GetNumberVecs(*ivec) > curDim_) {
1034 ivec = MVT::CloneView(*ivec,dimind);
1037 MVT::SetBlock(*ivec,dimind,*lclV);
1039 Teuchos::RCP<MV> tmpVecs;
1040 if (curDim_*2 <= blockSize_*numBlocks_) {
1042 std::vector<int> block2(curDim_);
1043 for (
int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
1044 tmpVecs = MVT::CloneViewNonConst(*V_,block2);
1048 tmpVecs = MVT::Clone(*V_,curDim_);
1053#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1054 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1056 OPT::Apply(*MOp_,*lclV,*tmpVecs);
1057 count_ApplyM_ += curDim_;
1061 if (auxVecs_.size() > 0) {
1062#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1063 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1066 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1067 int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
1069 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1072#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1073 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1076 int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
1078 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1083#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1084 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1086 OPT::Apply(*Op_,*lclV,*tmpVecs);
1087 count_ApplyOp_ += curDim_;
1091 lclKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1092 MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
1095 tmpVecs = Teuchos::null;
1099 if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
1100 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1101 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
1102 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(newstate.T->size()) != curDim_,
1103 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
1105 if (newstate.X != X_) {
1106 MVT::SetBlock(*newstate.X,bsind,*X_);
1109 std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1113 Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
1115#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1116 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1119 Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
1122 "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
1126#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1127 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1130 std::vector<int> order(curDim_);
1133 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);
1136 Utils::permuteVectors(order,S);
1140 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
1142#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1143 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1147 MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
1150 newstate.KX = Teuchos::null;
1151 newstate.MX = Teuchos::null;
1155 lclV = Teuchos::null;
1156 lclKK = Teuchos::null;
1159 if ( newstate.KX != Teuchos::null ) {
1160 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
1161 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1162 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*X_),
1163 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1164 if (newstate.KX != KX_) {
1165 MVT::SetBlock(*newstate.KX,bsind,*KX_);
1171#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1172 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1174 OPT::Apply(*Op_,*X_,*KX_);
1175 count_ApplyOp_ += blockSize_;
1178 newstate.R = Teuchos::null;
1183 if ( newstate.MX != Teuchos::null ) {
1184 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
1185 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1186 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*X_),
1187 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1188 if (newstate.MX != MX_) {
1189 MVT::SetBlock(*newstate.MX,bsind,*MX_);
1195#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1196 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1198 OPT::Apply(*MOp_,*X_,*MX_);
1199 count_ApplyOp_ += blockSize_;
1202 newstate.R = Teuchos::null;
1207 TEUCHOS_TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error,
"Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
1211 if (newstate.R != Teuchos::null) {
1212 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1213 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1214 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1215 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1216 if (newstate.R != R_) {
1217 MVT::SetBlock(*newstate.R,bsind,*R_);
1221#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1226 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1227 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1229 for (
int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1230 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1235 Rnorms_current_ =
false;
1236 R2norms_current_ =
false;
1239 initialized_ =
true;
1241 if (om_->isVerbosity(
Debug ) ) {
1251 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
1255 if (om_->isVerbosity(
Debug)) {
1256 currentStatus( om_->stream(
Debug) );
1281 if (initialized_ ==
false) {
1287 const int searchDim = blockSize_*numBlocks_;
1289 Teuchos::BLAS<int,ScalarType> blas;
1294 Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
1300 while (tester_->checkStatus(
this) !=
Passed && curDim_ < searchDim) {
1303 if (om_->isVerbosity(
Debug)) {
1304 currentStatus( om_->stream(
Debug) );
1313 std::vector<int> curind(blockSize_);
1314 for (
int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1315 H_ = MVT::CloneViewNonConst(*V_,curind);
1319 if (Prec_ != Teuchos::null) {
1320#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1321 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1323 OPT::Apply( *Prec_, *R_, *H_ );
1324 count_ApplyPrec_ += blockSize_;
1327 std::vector<int> bsind(blockSize_);
1328 for (
int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1329 MVT::SetBlock(*R_,bsind,*H_);
1336#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1337 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1339 OPT::Apply( *MOp_, *H_, *MH_);
1340 count_ApplyM_ += blockSize_;
1348 std::vector<int> prevind(curDim_);
1349 for (
int i=0; i<curDim_; ++i) prevind[i] = i;
1350 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
1354#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1355 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1358 Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
1359 against.push_back(Vprev);
1360 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1361 int rank = orthman_->projectAndNormalizeMat(*H_,against,
1365 "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
1372#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1373 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1375 OPT::Apply( *Op_, *H_, *KH_);
1376 count_ApplyOp_ += blockSize_;
1379 if (om_->isVerbosity(
Debug ) ) {
1384 om_->print(
Debug, accuracyCheck(chk,
": after ortho H") );
1391 om_->print(
OrthoDetails, accuracyCheck(chk,
": after ortho H") );
1396 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
1398 nextKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
1399 MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
1401 nextKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1402 MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
1405 nextKK = Teuchos::null;
1406 for (
int i=curDim_; i<curDim_+blockSize_; ++i) {
1407 for (
int j=0; j<i; ++j) {
1408 (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
1413 curDim_ += blockSize_;
1414 H_ = KH_ = MH_ = Teuchos::null;
1415 Vprev = Teuchos::null;
1417 if (om_->isVerbosity(
Debug ) ) {
1420 om_->print(
Debug, accuracyCheck(chk,
": after expanding KK") );
1424 curind.resize(curDim_);
1425 for (
int i=0; i<curDim_; ++i) curind[i] = i;
1426 Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
1430#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1431 Teuchos::TimeMonitor lcltimer(*timerDS_);
1433 int nevlocal = curDim_;
1434 int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
1435 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
1437 TEUCHOS_TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,
"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors.");
1442#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1443 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1446 std::vector<int> order(curDim_);
1449 sm_->sort(theta_, Teuchos::rcp(&order,
false), curDim_);
1452 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
1453 Utils::permuteVectors(order,curS);
1457 Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1461#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1462 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1464 MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
1469#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1470 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1472 OPT::Apply( *Op_, *X_, *KX_);
1473 count_ApplyOp_ += blockSize_;
1477#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1478 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1480 OPT::Apply(*MOp_,*X_,*MX_);
1481 count_ApplyM_ += blockSize_;
1490#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1491 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1494 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1496 for (
int i = 0; i < blockSize_; ++i) {
1499 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1503 Rnorms_current_ =
false;
1504 R2norms_current_ =
false;
1508 if (om_->isVerbosity(
Debug ) ) {
1516 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1524 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );