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
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();
359 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
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);
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>
625 setSize(blockSize,numBlocks_);
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>
675 return this->getRes2Norms();
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");
816 X_ = MVT::Clone(*tmp,blockSize_);
817 om_->print(
Debug,
" >> Allocating KX_\n");
818 KX_ = MVT::Clone(*tmp,blockSize_);
820 om_->print(
Debug,
" >> Allocating MX_\n");
821 MX_ = MVT::Clone(*tmp,blockSize_);
826 om_->print(
Debug,
" >> Allocating R_\n");
827 R_ = MVT::Clone(*tmp,blockSize_);
832 int newsd = blockSize_*numBlocks_;
833 theta_.resize(blockSize_*numBlocks_,NANVAL);
834 om_->print(
Debug,
" >> Allocating V_\n");
835 V_ = MVT::Clone(*tmp,newsd);
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) {
855 numAuxVecs_ += MVT::GetNumberVecs(**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) {
925 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
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().");
931 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
V) < newstate.
curDim, std::invalid_argument,
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_) {
949 if (curDim_ < MVT::GetNumberVecs(*newstate.
V)) {
950 newstate.
V = MVT::CloneView(*newstate.
V,nevind);
952 MVT::SetBlock(*newstate.
V,nevind,*V_);
954 lclV = MVT::CloneViewNonConst(*V_,nevind);
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;
989 curDim_ = MVT::GetNumberVecs(*ivec);
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;
1014 lclV = MVT::CloneViewNonConst(*V_,dimind);
1017 MVT::MvRandom(*lclV);
1020 if (MVT::GetNumberVecs(*ivec) > curDim_) {
1021 ivec = MVT::CloneView(*ivec,dimind);
1024 MVT::SetBlock(*ivec,dimind,*lclV);
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;
1031 tmpVecs = MVT::CloneViewNonConst(*V_,block2);
1035 tmpVecs = MVT::Clone(*V_,curDim_);
1040 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1041 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1043 OPT::Apply(*MOp_,*lclV,*tmpVecs);
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_ );
1073 OPT::Apply(*Op_,*lclV,*tmpVecs);
1074 count_ApplyOp_ += curDim_;
1078 lclKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1079 MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
1082 tmpVecs = Teuchos::null;
1086 if (newstate.
X != Teuchos::null && newstate.
T != Teuchos::null) {
1087 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
X) != blockSize_ || MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X_),
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_) {
1093 MVT::SetBlock(*newstate.
X,bsind,*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_ );
1106 Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
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_);
1123 Utils::permuteVectors(order,S);
1127 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
1129 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1130 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1134 MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
1137 newstate.
KX = Teuchos::null;
1138 newstate.
MX = Teuchos::null;
1142 lclV = Teuchos::null;
1143 lclKK = Teuchos::null;
1146 if ( newstate.
KX != Teuchos::null ) {
1147 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
KX) != blockSize_,
1148 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1149 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.
KX) != MVT::GetGlobalLength(*X_),
1150 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1151 if (newstate.
KX != KX_) {
1152 MVT::SetBlock(*newstate.
KX,bsind,*KX_);
1158 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1159 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1161 OPT::Apply(*Op_,*X_,*KX_);
1162 count_ApplyOp_ += blockSize_;
1165 newstate.
R = Teuchos::null;
1170 if ( newstate.
MX != Teuchos::null ) {
1171 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
MX) != blockSize_,
1172 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1173 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.
MX) != MVT::GetGlobalLength(*X_),
1174 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1175 if (newstate.
MX != MX_) {
1176 MVT::SetBlock(*newstate.
MX,bsind,*MX_);
1182 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1183 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1185 OPT::Apply(*MOp_,*X_,*MX_);
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) {
1199 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
R) != blockSize_,
1200 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1201 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*X_),
1202 std::invalid_argument,
"Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1203 if (newstate.
R != R_) {
1204 MVT::SetBlock(*newstate.
R,bsind,*R_);
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1209 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1213 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1214 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1216 for (
int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1217 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
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)) {
1243 currentStatus( om_->stream(
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)) {
1291 currentStatus( om_->stream(
Debug) );
1300 std::vector<int> curind(blockSize_);
1301 for (
int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1302 H_ = MVT::CloneViewNonConst(*V_,curind);
1306 if (Prec_ != Teuchos::null) {
1307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1308 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1310 OPT::Apply( *Prec_, *R_, *H_ );
1311 count_ApplyPrec_ += blockSize_;
1314 std::vector<int> bsind(blockSize_);
1315 for (
int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1316 MVT::SetBlock(*R_,bsind,*H_);
1323 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1324 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1326 OPT::Apply( *MOp_, *H_, *MH_);
1327 count_ApplyM_ += blockSize_;
1335 std::vector<int> prevind(curDim_);
1336 for (
int i=0; i<curDim_; ++i) prevind[i] = i;
1337 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
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_ );
1362 OPT::Apply( *Op_, *H_, *KH_);
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_) );
1386 MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
1388 nextKK = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1389 MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
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;
1413 Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
1417 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1418 Teuchos::TimeMonitor lcltimer(*timerDS_);
1420 int nevlocal = curDim_;
1421 int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
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_);
1440 Utils::permuteVectors(order,curS);
1444 Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1448 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1449 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1451 MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
1456 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1457 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1459 OPT::Apply( *Op_, *X_, *KX_);
1460 count_ApplyOp_ += blockSize_;
1464 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1465 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1467 OPT::Apply(*MOp_,*X_,*MX_);
1468 count_ApplyM_ += blockSize_;
1477 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1478 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1481 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1482 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1483 for (
int i = 0; i < blockSize_; ++i) {
1486 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
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) {
1540 MVT::MvNorm(*R_,R2norms_);
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;
1605 lclV = MVT::CloneView(*V_,lclind);
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_);
1616 OPT::Apply(*Op_,*lclV,*lclKV);
1617 MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
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_) {
1639 MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
1640 os <<
" >> Error in MX == M*X : " << err << endl;
1642 if (chk.checkKX && initialized_) {
1643 MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
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_) {
1661 MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
1662 os <<
" >> Error in KH == K*H : " << err << endl;
1664 if (chk.checkMH && hasM_ && initialized_) {
1665 MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
1666 os <<
" >> Error in MH == M*H : " << err << endl;
1670 if (chk.checkR && initialized_) {
1671 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1672 MVT::MvTransMv(ONE,*X_,*R_,xTx);
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.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const MV > H
The current preconditioned residual vectors.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
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...
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.
BlockDavidsonState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
BlockDavidsonInitFailure is thrown when the BlockDavidson solver is unable to generate an initial ite...
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...
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...
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.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
int curDim
The current dimension of the solver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Structure to contain pointers to BlockDavidson state variables.
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().
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
void iterate()
This method performs BlockDavidson iterations until the status test indicates the need to stop or an ...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
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.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
int getBlockSize() const
Get the blocksize used by the iterative solver.
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.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For BlockDavidson, this always returns n...