33 #ifndef ANASAZI_BLOCK_DAVIDSON_HPP 34 #define ANASAZI_BLOCK_DAVIDSON_HPP 41 #include "Teuchos_ScalarTraits.hpp" 46 #include "Teuchos_LAPACK.hpp" 47 #include "Teuchos_BLAS.hpp" 48 #include "Teuchos_SerialDenseMatrix.hpp" 49 #include "Teuchos_ParameterList.hpp" 50 #include "Teuchos_TimeMonitor.hpp" 76 template <
class ScalarType,
class MV>
87 Teuchos::RCP<const MV>
V;
89 Teuchos::RCP<const MV>
X;
91 Teuchos::RCP<const MV>
KX;
93 Teuchos::RCP<const MV>
MX;
95 Teuchos::RCP<const MV>
R;
100 Teuchos::RCP<const MV>
H;
102 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >
T;
108 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
KK;
110 X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
111 R(Teuchos::null), H(Teuchos::null),
112 T(Teuchos::null), KK(Teuchos::null) {}
150 template <
class ScalarType,
class MV,
class OP>
164 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
168 Teuchos::ParameterList ¶ms
259 bool isInitialized()
const;
282 int getNumIters()
const;
285 void resetNumIters();
294 Teuchos::RCP<const MV> getRitzVectors();
301 std::vector<Value<ScalarType> > getRitzValues();
311 std::vector<int> getRitzIndex();
319 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
327 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
337 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
344 int getCurSubspaceDim()
const;
347 int getMaxSubspaceDim()
const;
359 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest()
const;
373 void setBlockSize(
int blockSize);
376 int getBlockSize()
const;
390 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
393 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs()
const;
409 void setSize(
int blockSize,
int numBlocks);
417 void currentStatus(std::ostream &os);
428 typedef Teuchos::ScalarTraits<ScalarType> SCT;
429 typedef typename SCT::magnitudeType MagnitudeType;
430 const MagnitudeType ONE;
431 const MagnitudeType ZERO;
432 const MagnitudeType NANVAL;
438 bool checkX, checkMX, checkKX;
439 bool checkH, checkMH, checkKH;
442 CheckList() : checkV(
false),
443 checkX(
false),checkMX(
false),checkKX(
false),
444 checkH(
false),checkMH(
false),checkKH(
false),
445 checkR(
false),checkQ(
false),checkKK(
false) {};
450 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
454 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
455 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
456 const Teuchos::RCP<OutputManager<ScalarType> > om_;
457 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
458 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
462 Teuchos::RCP<const OP> Op_;
463 Teuchos::RCP<const OP> MOp_;
464 Teuchos::RCP<const OP> Prec_;
469 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
470 timerSortEval_, timerDS_,
471 timerLocal_, timerCompRes_,
472 timerOrtho_, timerInit_;
476 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
505 Teuchos::RCP<MV> X_, KX_, MX_, R_,
511 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
514 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
521 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
524 bool Rnorms_current_, R2norms_current_;
539 template <
class ScalarType,
class MV,
class OP>
542 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
546 Teuchos::ParameterList ¶ms
548 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
549 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
550 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
558 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
559 timerOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation Op*x")),
560 timerMOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation M*x")),
561 timerPrec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Operation Prec*x")),
562 timerSortEval_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Sorting eigenvalues")),
563 timerDS_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Direct solve")),
564 timerLocal_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Local update")),
565 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Computing residuals")),
566 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Orthogonalization")),
567 timerInit_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidson::Initialization")),
577 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
580 Rnorms_current_(false),
581 R2norms_current_(false)
583 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
584 "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
585 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
586 "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
587 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
588 "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
589 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
590 "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
591 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
592 "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
593 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
594 "Anasazi::BlockDavidson::constructor: problem is not set.");
595 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
596 "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
599 Op_ = problem_->getOperator();
600 TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
601 "Anasazi::BlockDavidson::constructor: problem provides no operator.");
602 MOp_ = problem_->getM();
603 Prec_ = problem_->getPrec();
604 hasM_ = (MOp_ != Teuchos::null);
607 int bs = params.get(
"Block Size", problem_->getNEV());
608 int nb = params.get(
"Num Blocks", 2);
615 template <
class ScalarType,
class MV,
class OP>
622 template <
class ScalarType,
class MV,
class OP>
631 template <
class ScalarType,
class MV,
class OP>
640 template <
class ScalarType,
class MV,
class OP>
648 template <
class ScalarType,
class MV,
class OP>
656 template <
class ScalarType,
class MV,
class OP>
658 return blockSize_*numBlocks_;
663 template <
class ScalarType,
class MV,
class OP>
665 if (!initialized_)
return 0;
672 template <
class ScalarType,
class MV,
class OP>
673 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
681 template <
class ScalarType,
class MV,
class OP>
683 std::vector<int> ret(curDim_,0);
690 template <
class ScalarType,
class MV,
class OP>
692 std::vector<Value<ScalarType> > ret(curDim_);
693 for (
int i=0; i<curDim_; ++i) {
694 ret[i].realpart = theta_[i];
695 ret[i].imagpart = ZERO;
703 template <
class ScalarType,
class MV,
class OP>
711 template <
class ScalarType,
class MV,
class OP>
719 template <
class ScalarType,
class MV,
class OP>
727 template <
class ScalarType,
class MV,
class OP>
738 state.
MX = Teuchos::null;
744 state.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
746 state.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(0));
754 template <
class ScalarType,
class MV,
class OP>
760 template <
class ScalarType,
class MV,
class OP>
764 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 765 Teuchos::TimeMonitor initimer( *timerInit_ );
771 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument,
"Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
772 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument,
"Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
773 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
778 blockSize_ = blockSize;
779 numBlocks_ = numBlocks;
781 Teuchos::RCP<const MV> tmp;
786 if (X_ != Teuchos::null) {
790 tmp = problem_->getInitVec();
791 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
792 "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
795 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) >
MVT::GetGlobalLength(*tmp),std::invalid_argument,
796 "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
803 Rnorms_.resize(blockSize_,NANVAL);
804 R2norms_.resize(blockSize_,NANVAL);
815 om_->print(
Debug,
" >> Allocating X_\n");
817 om_->print(
Debug,
" >> Allocating KX_\n");
820 om_->print(
Debug,
" >> Allocating MX_\n");
826 om_->print(
Debug,
" >> Allocating R_\n");
832 int newsd = blockSize_*numBlocks_;
833 theta_.resize(blockSize_*numBlocks_,NANVAL);
834 om_->print(
Debug,
" >> Allocating V_\n");
836 KK_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
838 om_->print(
Debug,
" >> done allocating.\n");
840 initialized_ =
false;
847 template <
class ScalarType,
class MV,
class OP>
849 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
854 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
859 if (numAuxVecs_ > 0 && initialized_) {
860 initialized_ =
false;
863 if (om_->isVerbosity(
Debug ) ) {
866 om_->print(
Debug, accuracyCheck(chk,
": in setAuxVecs()") );
884 template <
class ScalarType,
class MV,
class OP>
890 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 891 Teuchos::TimeMonitor inittimer( *timerInit_ );
894 std::vector<int> bsind(blockSize_);
895 for (
int i=0; i<blockSize_; ++i) bsind[i] = i;
897 Teuchos::BLAS<int,ScalarType> blas;
921 Teuchos::RCP<MV> lclV;
922 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
924 if (newstate.
V != Teuchos::null && newstate.
KK != Teuchos::null) {
926 "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
927 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim < blockSize_, std::invalid_argument,
928 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
929 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim > blockSize_*numBlocks_, std::invalid_argument,
930 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
932 "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
934 curDim_ = newstate.
curDim;
936 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
938 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.
curDim, std::invalid_argument,
939 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
942 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
KK->numRows() < curDim_ || newstate.
KK->numCols() < curDim_, std::invalid_argument,
943 "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
946 std::vector<int> nevind(curDim_);
947 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
948 if (newstate.
V != V_) {
957 lclKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
958 if (newstate.
KK != KK_) {
959 if (newstate.
KK->numRows() > curDim_ || newstate.
KK->numCols() > curDim_) {
960 newstate.
KK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.
KK,curDim_,curDim_) );
962 lclKK->assign(*newstate.
KK);
966 for (
int j=0; j<curDim_-1; ++j) {
967 for (
int i=j+1; i<curDim_; ++i) {
968 (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
975 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
976 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
977 "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
979 newstate.
X = Teuchos::null;
980 newstate.
MX = Teuchos::null;
981 newstate.
KX = Teuchos::null;
982 newstate.
R = Teuchos::null;
983 newstate.
H = Teuchos::null;
984 newstate.
T = Teuchos::null;
985 newstate.
KK = Teuchos::null;
986 newstate.
V = Teuchos::null;
991 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
992 if (curDim_ > blockSize_*numBlocks_) {
995 curDim_ = blockSize_*numBlocks_;
997 bool userand =
false;
1002 curDim_ = blockSize_;
1012 std::vector<int> dimind(curDim_);
1013 for (
int i=0; i<curDim_; ++i) dimind[i] = i;
1026 Teuchos::RCP<MV> tmpVecs;
1027 if (curDim_*2 <= blockSize_*numBlocks_) {
1029 std::vector<int> block2(curDim_);
1030 for (
int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
1040 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1041 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1044 count_ApplyM_ += curDim_;
1048 if (auxVecs_.size() > 0) {
1049 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1050 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1053 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1054 int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
1056 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1059 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1060 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1063 int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
1065 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1070 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1071 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1074 count_ApplyOp_ += curDim_;
1078 lclKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1082 tmpVecs = Teuchos::null;
1086 if (newstate.
X != Teuchos::null && newstate.
T != Teuchos::null) {
1088 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
1089 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(newstate.
T->size()) != curDim_,
1090 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
1092 if (newstate.
X != X_) {
1096 std::copy(newstate.
T->begin(),newstate.
T->end(),theta_.begin());
1100 Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
1102 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1103 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1109 "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
1113 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1114 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1117 std::vector<int> order(curDim_);
1120 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);
1127 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
1129 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1130 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1137 newstate.
KX = Teuchos::null;
1138 newstate.
MX = Teuchos::null;
1142 lclV = Teuchos::null;
1143 lclKK = Teuchos::null;
1146 if ( newstate.
KX != Teuchos::null ) {
1148 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1150 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1151 if (newstate.
KX != KX_) {
1158 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1159 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1162 count_ApplyOp_ += blockSize_;
1165 newstate.
R = Teuchos::null;
1170 if ( newstate.
MX != Teuchos::null ) {
1172 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1174 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1175 if (newstate.
MX != MX_) {
1182 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1183 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1186 count_ApplyOp_ += blockSize_;
1189 newstate.
R = Teuchos::null;
1194 TEUCHOS_TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error,
"Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
1198 if (newstate.
R != Teuchos::null) {
1200 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1202 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1203 if (newstate.
R != R_) {
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1209 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1214 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1216 for (
int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1222 Rnorms_current_ =
false;
1223 R2norms_current_ =
false;
1226 initialized_ =
true;
1228 if (om_->isVerbosity(
Debug ) ) {
1238 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
1242 if (om_->isVerbosity(
Debug)) {
1253 template <
class ScalarType,
class MV,
class OP>
1263 template <
class ScalarType,
class MV,
class OP>
1268 if (initialized_ ==
false) {
1274 const int searchDim = blockSize_*numBlocks_;
1276 Teuchos::BLAS<int,ScalarType> blas;
1281 Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
1287 while (tester_->checkStatus(
this) !=
Passed && curDim_ < searchDim) {
1290 if (om_->isVerbosity(
Debug)) {
1300 std::vector<int> curind(blockSize_);
1301 for (
int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1306 if (Prec_ != Teuchos::null) {
1307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1308 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1311 count_ApplyPrec_ += blockSize_;
1314 std::vector<int> bsind(blockSize_);
1315 for (
int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1323 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1324 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1327 count_ApplyM_ += blockSize_;
1335 std::vector<int> prevind(curDim_);
1336 for (
int i=0; i<curDim_; ++i) prevind[i] = i;
1341 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1342 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1345 Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
1346 against.push_back(Vprev);
1347 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1348 int rank = orthman_->projectAndNormalizeMat(*H_,against,
1352 "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
1359 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1360 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1363 count_ApplyOp_ += blockSize_;
1366 if (om_->isVerbosity(
Debug ) ) {
1371 om_->print(
Debug, accuracyCheck(chk,
": after ortho H") );
1378 om_->print(
OrthoDetails, accuracyCheck(chk,
": after ortho H") );
1383 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
1385 nextKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
1388 nextKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1392 nextKK = Teuchos::null;
1393 for (
int i=curDim_; i<curDim_+blockSize_; ++i) {
1394 for (
int j=0; j<i; ++j) {
1395 (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
1400 curDim_ += blockSize_;
1401 H_ = KH_ = MH_ = Teuchos::null;
1402 Vprev = Teuchos::null;
1404 if (om_->isVerbosity(
Debug ) ) {
1407 om_->print(
Debug, accuracyCheck(chk,
": after expanding KK") );
1411 curind.resize(curDim_);
1412 for (
int i=0; i<curDim_; ++i) curind[i] = i;
1417 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1418 Teuchos::TimeMonitor lcltimer(*timerDS_);
1420 int nevlocal = curDim_;
1422 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
1424 TEUCHOS_TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,
"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors.");
1429 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1430 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1433 std::vector<int> order(curDim_);
1436 sm_->sort(theta_, Teuchos::rcp(&order,
false), curDim_);
1439 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
1444 Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1448 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1449 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1456 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1457 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1460 count_ApplyOp_ += blockSize_;
1464 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1465 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1468 count_ApplyM_ += blockSize_;
1477 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1478 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1482 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1483 for (
int i = 0; i < blockSize_; ++i) {
1490 Rnorms_current_ =
false;
1491 R2norms_current_ =
false;
1495 if (om_->isVerbosity(
Debug ) ) {
1503 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1511 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );
1521 template <
class ScalarType,
class MV,
class OP>
1522 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1524 if (Rnorms_current_ ==
false) {
1526 orthman_->norm(*R_,Rnorms_);
1527 Rnorms_current_ =
true;
1535 template <
class ScalarType,
class MV,
class OP>
1536 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1538 if (R2norms_current_ ==
false) {
1541 R2norms_current_ =
true;
1548 template <
class ScalarType,
class MV,
class OP>
1550 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1551 "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
1557 template <
class ScalarType,
class MV,
class OP>
1589 template <
class ScalarType,
class MV,
class OP>
1594 std::stringstream os;
1596 os.setf(std::ios::scientific, std::ios::floatfield);
1598 os <<
" Debugging checks: iteration " << iter_ << where << endl;
1601 std::vector<int> lclind(curDim_);
1602 for (
int i=0; i<curDim_; ++i) lclind[i] = i;
1603 Teuchos::RCP<const MV> lclV;
1607 if (chk.checkV && initialized_) {
1608 MagnitudeType err = orthman_->orthonormError(*lclV);
1609 os <<
" >> Error in V^H M V == I : " << err << endl;
1610 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1611 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
1612 os <<
" >> Error in V^H M Q[" << i <<
"] == 0 : " << err << endl;
1614 Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
1615 Teuchos::RCP<MV> lclKV =
MVT::Clone(*V_,curDim_);
1618 Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
1621 for (
int j=0; j<curDim_; ++j) {
1622 for (
int i=j+1; i<curDim_; ++i) {
1623 curKK(i,j) = curKK(j,i);
1626 os <<
" >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
1630 if (chk.checkX && initialized_) {
1631 MagnitudeType err = orthman_->orthonormError(*X_);
1632 os <<
" >> Error in X^H M X == I : " << err << endl;
1633 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1634 err = orthman_->orthogError(*X_,*auxVecs_[i]);
1635 os <<
" >> Error in X^H M Q[" << i <<
"] == 0 : " << err << endl;
1638 if (chk.checkMX && hasM_ && initialized_) {
1640 os <<
" >> Error in MX == M*X : " << err << endl;
1642 if (chk.checkKX && initialized_) {
1644 os <<
" >> Error in KX == K*X : " << err << endl;
1648 if (chk.checkH && initialized_) {
1649 MagnitudeType err = orthman_->orthonormError(*H_);
1650 os <<
" >> Error in H^H M H == I : " << err << endl;
1651 err = orthman_->orthogError(*H_,*lclV);
1652 os <<
" >> Error in H^H M V == 0 : " << err << endl;
1653 err = orthman_->orthogError(*H_,*X_);
1654 os <<
" >> Error in H^H M X == 0 : " << err << endl;
1655 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1656 err = orthman_->orthogError(*H_,*auxVecs_[i]);
1657 os <<
" >> Error in H^H M Q[" << i <<
"] == 0 : " << err << endl;
1660 if (chk.checkKH && initialized_) {
1662 os <<
" >> Error in KH == K*H : " << err << endl;
1664 if (chk.checkMH && hasM_ && initialized_) {
1666 os <<
" >> Error in MH == M*H : " << err << endl;
1670 if (chk.checkR && initialized_) {
1671 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1673 MagnitudeType err = xTx.normFrobenius();
1674 os <<
" >> Error in X^H R == 0 : " << err << endl;
1678 if (chk.checkKK && initialized_) {
1679 Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
1680 for (
int j=0; j<curDim_; ++j) {
1681 for (
int i=0; i<curDim_; ++i) {
1682 SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
1685 os <<
" >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
1690 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1691 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
1692 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << i <<
"] == I : " << err << endl;
1693 for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
1694 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1695 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << j <<
"] == 0 : " << err << endl;
1708 template <
class ScalarType,
class MV,
class OP>
1714 os.setf(std::ios::scientific, std::ios::floatfield);
1717 os <<
"================================================================================" << endl;
1719 os <<
" BlockDavidson Solver Status" << endl;
1721 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
1722 os <<
"The number of iterations performed is " <<iter_<<endl;
1723 os <<
"The block size is " << blockSize_<<endl;
1724 os <<
"The number of blocks is " << numBlocks_<<endl;
1725 os <<
"The current basis size is " << curDim_<<endl;
1726 os <<
"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
1727 os <<
"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1728 os <<
"The number of operations M*x is "<<count_ApplyM_<<endl;
1729 os <<
"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
1731 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1735 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1736 os << std::setw(20) <<
"Eigenvalue" 1737 << std::setw(20) <<
"Residual(M)" 1738 << std::setw(20) <<
"Residual(2)" 1740 os <<
"--------------------------------------------------------------------------------"<<endl;
1741 for (
int i=0; i<blockSize_; ++i) {
1742 os << std::setw(20) << theta_[i];
1743 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
1744 else os << std::setw(20) <<
"not current";
1745 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
1746 else os << std::setw(20) <<
"not current";
1750 os <<
"================================================================================" << endl;
virtual ~BlockDavidson()
Anasazi::BlockDavidson destructor.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Teuchos::RCP< const MV > H
The current preconditioned residual vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
int getNumIters() const
Get the current iteration count.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
This class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Declaration of basic traits for the multivector type.
BlockDavidsonOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the...
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to the given output stream...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
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...
BlockDavidsonInitFailure is thrown when the BlockDavidson solver is unable to generate an initial ite...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
An exception class parent to all Anasazi exceptions.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
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...
Anasazi's templated, static class providing utilities for the solvers.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For BlockDavidson, this always returns n...
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::RCP< const MV > V
The basis for the Krylov space.
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.
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.
int getBlockSize() const
Get the blocksize used by the iterative solver.
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).
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 ...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
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).
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
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 curDim
The current dimension of the solver.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Structure to contain pointers to BlockDavidson state variables.
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...
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Types and exceptions used within Anasazi solvers and interfaces.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
void iterate()
This method performs BlockDavidson iterations until the status test indicates the need to stop or an ...
BlockDavidsonState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
BlockDavidson(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)
BlockDavidson constructor with eigenproblem, solver utilities, and parameter list of solver options...
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms, computing the norms if they are not up-to-date with the current res...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void setBlockSize(int blockSize)
Set the blocksize.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
void resetNumIters()
Reset the iteration count.
Class which provides internal utilities for the Anasazi solvers.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms, computing the norms if they are not up-to-date with the current resid...
Teuchos::RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.