33 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_HPP 34 #define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP 41 #include "Teuchos_ScalarTraits.hpp" 42 #include "AnasaziHelperTraits.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" 75 template <
class ScalarType,
class MulVec>
83 Teuchos::RCP<const MulVec>
V;
89 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
H;
91 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
S;
93 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
Q;
95 H(Teuchos::null),
S(Teuchos::null),
133 template <
class ScalarType,
class MV,
class OP>
151 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
155 Teuchos::ParameterList ¶ms
273 std::vector<Value<ScalarType> > ret = ritzValues_;
274 ret.resize(ritzIndex_.size());
289 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms() {
290 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
298 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms() {
299 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
307 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms() {
308 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
309 ret.resize(ritzIndex_.size());
322 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
333 void setSize(
int blockSize,
int numBlocks);
359 if (!initialized_)
return 0;
364 int getMaxSubspaceDim()
const {
return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
379 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
382 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const {
return auxVecs_;}
428 typedef Teuchos::ScalarTraits<ScalarType> SCT;
429 typedef typename SCT::magnitudeType MagnitudeType;
430 typedef typename std::vector<ScalarType>::iterator STiter;
431 typedef typename std::vector<MagnitudeType>::iterator MTiter;
432 const MagnitudeType MT_ONE;
433 const MagnitudeType MT_ZERO;
434 const MagnitudeType NANVAL;
435 const ScalarType ST_ONE;
436 const ScalarType ST_ZERO;
444 CheckList() : checkV(false), checkArn(false), checkAux(false) {};
449 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
450 void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
451 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
452 std::vector<int>& order );
456 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
457 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
458 const Teuchos::RCP<OutputManager<ScalarType> > om_;
459 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
460 const Teuchos::RCP<OrthoManager<ScalarType,MV> > orthman_;
464 Teuchos::RCP<const OP> Op_;
468 Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
469 timerCompSF_, timerSortSF_,
470 timerCompRitzVec_, timerOrtho_;
504 Teuchos::RCP<MV> ritzVectors_, V_;
510 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
515 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > schurH_;
516 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Q_;
519 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
526 bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
529 std::vector<Value<ScalarType> > ritzValues_;
530 std::vector<MagnitudeType> ritzResiduals_;
533 std::vector<int> ritzIndex_;
534 std::vector<int> ritzOrder_;
543 template <
class ScalarType,
class MV,
class OP>
546 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
550 Teuchos::ParameterList ¶ms
552 MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
553 MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
554 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
555 ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
556 ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
564 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
565 timerOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Operation Op*x")),
566 timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Sorting Ritz values")),
567 timerCompSF_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Computing Schur form")),
568 timerSortSF_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Sorting Schur form")),
569 timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
570 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Orthogonalization")),
580 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
583 ritzVecsCurrent_(false),
584 ritzValsCurrent_(false),
585 schurCurrent_(false),
588 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
589 "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
590 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
591 "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
592 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
593 "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
594 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
595 "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
596 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
597 "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
598 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
599 "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
600 TEUCHOS_TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
601 "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
602 TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
603 "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
604 TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
605 "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
606 TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
607 "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
610 Op_ = problem_->getOperator();
613 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Step Size"), std::invalid_argument,
614 "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
615 int ss = params.get(
"Step Size",numBlocks_);
619 int bs = params.get(
"Block Size", 1);
620 int nb = params.get(
"Num Blocks", 3*problem_->getNEV());
625 int numRitzVecs = params.get(
"Number of Ritz Vectors", 0);
629 numRitzPrint_ = params.get(
"Print Number of Ritz Values", bs);
636 template <
class ScalarType,
class MV,
class OP>
639 setSize(blockSize,numBlocks_);
645 template <
class ScalarType,
class MV,
class OP>
648 TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
649 stepSize_ = stepSize;
654 template <
class ScalarType,
class MV,
class OP>
660 TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
663 if (numRitzVecs != numRitzVecs_) {
665 ritzVectors_ = Teuchos::null;
666 ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
668 ritzVectors_ = Teuchos::null;
670 numRitzVecs_ = numRitzVecs;
671 ritzVecsCurrent_ =
false;
677 template <
class ScalarType,
class MV,
class OP>
683 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
684 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
685 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
690 blockSize_ = blockSize;
691 numBlocks_ = numBlocks;
693 Teuchos::RCP<const MV> tmp;
699 if (problem_->getInitVec() != Teuchos::null) {
700 tmp = problem_->getInitVec();
704 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
705 "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
713 if (problem_->isHermitian()) {
714 newsd = blockSize_*numBlocks_;
716 newsd = blockSize_*numBlocks_+1;
719 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(newsd) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
720 "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
722 ritzValues_.resize(newsd);
723 ritzResiduals_.resize(newsd,MT_ONE);
724 ritzOrder_.resize(newsd);
726 V_ = MVT::Clone(*tmp,newsd+blockSize_);
727 H_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
728 Q_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
730 initialized_ =
false;
737 template <
class ScalarType,
class MV,
class OP>
739 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
744 if (om_->isVerbosity(
Debug ) ) {
748 om_->print(
Debug, accuracyCheck(chk,
": in setAuxVecs()") );
752 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
753 numAuxVecs_ += MVT::GetNumberVecs(**i);
757 if (numAuxVecs_ > 0 && initialized_) {
758 initialized_ =
false;
771 template <
class ScalarType,
class MV,
class OP>
776 std::vector<int> bsind(blockSize_);
777 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
785 std::string errstr(
"Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
789 if (newstate.
V != Teuchos::null && newstate.
H != Teuchos::null) {
793 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
V) != MVT::GetGlobalLength(*V_),
794 std::invalid_argument, errstr );
795 if (newstate.
V != V_) {
796 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
V) < blockSize_,
797 std::invalid_argument, errstr );
798 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
V) > getMaxSubspaceDim(),
799 std::invalid_argument, errstr );
801 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim > getMaxSubspaceDim(),
802 std::invalid_argument, errstr );
804 curDim_ = newstate.
curDim;
805 int lclDim = MVT::GetNumberVecs(*newstate.
V);
808 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
H->numRows() < curDim_ || newstate.
H->numCols() < curDim_, std::invalid_argument, errstr);
810 if (curDim_ == 0 && lclDim > blockSize_) {
811 om_->stream(
Warnings) <<
"Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
812 <<
"The block size however is only " << blockSize_ << std::endl
813 <<
"The last " << lclDim - blockSize_ <<
" vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
818 if (newstate.
V != V_) {
819 std::vector<int> nevind(lclDim);
820 for (
int i=0; i<lclDim; i++) nevind[i] = i;
821 MVT::SetBlock(*newstate.
V,nevind,*V_);
825 if (newstate.
H != H_) {
826 H_->putScalar( ST_ZERO );
827 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.
H,curDim_+blockSize_,curDim_);
828 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
829 lclH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
833 lclH = Teuchos::null;
840 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
841 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
842 "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
844 int lclDim = MVT::GetNumberVecs(*ivec);
845 bool userand =
false;
846 if (lclDim < blockSize_) {
854 std::vector<int> dimind2(lclDim);
855 for (
int i=0; i<lclDim; i++) { dimind2[i] = i; }
858 Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,dimind2);
861 MVT::SetBlock(*ivec,dimind2,*newV1);
864 dimind2.resize(blockSize_-lclDim);
865 for (
int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
868 Teuchos::RCP<MV> newV2 = MVT::CloneViewNonConst(*V_,dimind2);
869 MVT::MvRandom(*newV2);
873 Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,bsind);
876 Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
879 MVT::SetBlock(*ivecV,bsind,*newV1);
883 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*V_,bsind);
886 if (auxVecs_.size() > 0) {
887 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 888 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
891 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
892 int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
894 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
897 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 898 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
901 int rank = orthman_->normalize(*newV);
903 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
910 newV = Teuchos::null;
914 ritzVecsCurrent_ =
false;
915 ritzValsCurrent_ =
false;
916 schurCurrent_ =
false;
921 if (om_->isVerbosity(
Debug ) ) {
927 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
931 if (om_->isVerbosity(
Debug)) {
932 currentStatus( om_->stream(
Debug) );
942 template <
class ScalarType,
class MV,
class OP>
952 template <
class ScalarType,
class MV,
class OP>
958 if (initialized_ ==
false) {
964 int searchDim = blockSize_*numBlocks_;
965 if (problem_->isHermitian() ==
false) {
974 while (tester_->checkStatus(
this) !=
Passed && curDim_+blockSize_ <= searchDim) {
979 int lclDim = curDim_ + blockSize_;
982 std::vector<int> curind(blockSize_);
983 for (
int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
984 Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
988 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
989 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
995 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 996 Teuchos::TimeMonitor lcltimer( *timerOp_ );
998 OPT::Apply(*Op_,*Vprev,*Vnext);
999 count_ApplyOp_ += blockSize_;
1003 Vprev = Teuchos::null;
1007 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1008 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1012 std::vector<int> prevind(lclDim);
1013 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
1014 Vprev = MVT::CloneView(*V_,prevind);
1015 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
1018 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1019 subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
1020 ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
1021 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
1022 AsubH.append( subH );
1025 if (auxVecs_.size() > 0) {
1026 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1027 AVprev.append( auxVecs_[i] );
1028 AsubH.append( Teuchos::null );
1035 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1036 subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
1037 ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
1038 int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
1045 "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
1051 Vnext = Teuchos::null;
1052 curDim_ += blockSize_;
1054 ritzVecsCurrent_ =
false;
1055 ritzValsCurrent_ =
false;
1056 schurCurrent_ =
false;
1059 if (!(iter_%stepSize_)) {
1060 computeRitzValues();
1064 if (om_->isVerbosity(
Debug ) ) {
1068 chk.checkArn =
true;
1069 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1074 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );
1078 if (om_->isVerbosity(
Debug)) {
1079 currentStatus( om_->stream(
Debug) );
1103 template <
class ScalarType,
class MV,
class OP>
1106 std::stringstream os;
1108 os.setf(std::ios::scientific, std::ios::floatfield);
1111 os <<
" Debugging checks: iteration " << iter_ << where << std::endl;
1114 std::vector<int> lclind(curDim_);
1115 for (
int i=0; i<curDim_; i++) lclind[i] = i;
1116 std::vector<int> bsind(blockSize_);
1117 for (
int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
1119 Teuchos::RCP<const MV> lclV,lclF;
1120 Teuchos::RCP<MV> lclAV;
1122 lclV = MVT::CloneView(*V_,lclind);
1123 lclF = MVT::CloneView(*V_,bsind);
1127 tmp = orthman_->orthonormError(*lclV);
1128 os <<
" >> Error in V^H M V == I : " << tmp << std::endl;
1130 tmp = orthman_->orthonormError(*lclF);
1131 os <<
" >> Error in F^H M F == I : " << tmp << std::endl;
1133 tmp = orthman_->orthogError(*lclV,*lclF);
1134 os <<
" >> Error in V^H M F == 0 : " << tmp << std::endl;
1136 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1138 tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
1139 os <<
" >> Error in V^H M Aux[" << i <<
"] == 0 : " << tmp << std::endl;
1141 tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
1142 os <<
" >> Error in F^H M Aux[" << i <<
"] == 0 : " << tmp << std::endl;
1150 lclAV = MVT::Clone(*V_,curDim_);
1152 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1153 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1155 OPT::Apply(*Op_,*lclV,*lclAV);
1159 Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
1160 MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
1163 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
1164 blockSize_,curDim_, curDim_ );
1165 MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
1168 std::vector<MagnitudeType> arnNorms( curDim_ );
1169 orthman_->norm( *lclAV, arnNorms );
1171 for (
int i=0; i<curDim_; i++) {
1172 os <<
" >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i <<
"]|| : " << arnNorms[i] << std::endl;
1178 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1179 tmp = orthman_->orthonormError(*auxVecs_[i]);
1180 os <<
" >> Error in Aux[" << i <<
"]^H M Aux[" << i <<
"] == I : " << tmp << std::endl;
1181 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1182 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1183 os <<
" >> Error in Aux[" << i <<
"]^H M Aux[" << j <<
"] == 0 : " << tmp << std::endl;
1204 template <
class ScalarType,
class MV,
class OP>
1212 if (!ritzValsCurrent_) {
1215 computeSchurForm(
false );
1232 template <
class ScalarType,
class MV,
class OP>
1235 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1236 Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
1239 TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
1240 "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
1242 TEUCHOS_TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
1243 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
1247 if (curDim_ && initialized_) {
1250 if (!ritzVecsCurrent_) {
1253 if (!schurCurrent_) {
1256 computeSchurForm(
true );
1261 TEUCHOS_TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
1262 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1264 Teuchos::LAPACK<int,ScalarType> lapack;
1265 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1275 std::vector<int> curind( curDim_ );
1276 for (
int i=0; i<curDim_; i++) { curind[i] = i; }
1277 Teuchos::RCP<const MV> Vtemp = MVT::CloneView( *V_, curind );
1278 if (problem_->isHermitian()) {
1280 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
1283 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
1286 bool complexRitzVals =
false;
1287 for (
int i=0; i<numRitzVecs_; i++) {
1288 if (ritzIndex_[i] != 0) {
1289 complexRitzVals =
true;
1293 if (complexRitzVals)
1294 om_->stream(
Warnings) <<
" Eigenproblem is Hermitian and complex eigenvalues have converged, corresponding eigenvectors will be incorrect!!!" 1299 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1302 Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
1305 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
1313 int lwork = 3*curDim_;
1314 std::vector<ScalarType> work( lwork );
1315 std::vector<MagnitudeType> rwork( curDim_ );
1319 ScalarType vl[ ldvl ];
1320 Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
1321 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1322 copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1323 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1324 "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1327 Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
1330 curind.resize( numRitzVecs_ );
1331 Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
1332 MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
1335 std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
1336 MVT::MvNorm( *view_ritzVectors, ritzNrm );
1339 tmpritzVectors_ = Teuchos::null;
1340 view_ritzVectors = Teuchos::null;
1343 ScalarType ritzScale = ST_ONE;
1344 for (
int i=0; i<numRitzVecs_; i++) {
1347 if (ritzIndex_[i] == 1 ) {
1348 ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
1349 std::vector<int> newind(2);
1350 newind[0] = i; newind[1] = i+1;
1351 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1352 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1353 MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1360 std::vector<int> newind(1);
1362 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1363 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1364 MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1371 ritzVecsCurrent_ =
true;
1380 template <
class ScalarType,
class MV,
class OP>
1382 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1383 "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
1389 template <
class ScalarType,
class MV,
class OP>
1407 template <
class ScalarType,
class MV,
class OP>
1411 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1412 Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
1419 if (!schurCurrent_) {
1423 if (!ritzValsCurrent_) {
1424 Teuchos::LAPACK<int,ScalarType> lapack;
1425 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1426 Teuchos::BLAS<int,ScalarType> blas;
1427 Teuchos::BLAS<int,MagnitudeType> blas_mag;
1430 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1433 schurH_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
1447 int lwork = 3*curDim_;
1448 std::vector<ScalarType> work( lwork );
1449 std::vector<MagnitudeType> rwork( curDim_ );
1450 std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
1451 std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
1452 std::vector<int> bwork( curDim_ );
1453 int info = 0, sdim = 0;
1455 lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
1456 &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork,
1457 &rwork[0], &bwork[0], &info );
1459 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1460 "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1466 bool hermImagDetected =
false;
1467 if (problem_->isHermitian()) {
1468 for (
int i=0; i<curDim_; i++)
1470 if (tmp_iRitzValues[i] != MT_ZERO)
1472 hermImagDetected =
true;
1476 if (hermImagDetected) {
1478 om_->stream(
Warnings) <<
" Eigenproblem is Hermitian and complex eigenvalues have been detected!!!" << std::endl;
1480 Teuchos::SerialDenseMatrix<int,ScalarType> localH( Teuchos::View, *H_, curDim_, curDim_ );
1481 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > tLocalH;
1482 if (Teuchos::ScalarTraits<ScalarType>::isComplex)
1483 tLocalH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::CONJ_TRANS ) );
1485 tLocalH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::TRANS ) );
1486 (*tLocalH) -= localH;
1487 MagnitudeType normF = tLocalH->normFrobenius();
1488 MagnitudeType norm1 = tLocalH->normOne();
1489 om_->stream(
Warnings) <<
" Symmetry error in projected eigenproblem: || S - S' ||_F = " << normF
1490 <<
", || S - S' ||_1 = " << norm1 << std::endl;
1508 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
1509 blockSize_, curDim_, curDim_ );
1512 Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
1513 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1514 curB.values(), curB.stride(), subQ.values(), subQ.stride(),
1515 ST_ZERO, subB.values(), subB.stride() );
1519 ScalarType* b_ptr = subB.values();
1520 if (problem_->isHermitian() && !hermImagDetected) {
1524 for (
int i=0; i<curDim_ ; i++) {
1525 ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
1534 ScalarType vl[ ldvl ];
1535 Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
1536 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1537 S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1539 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1540 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1548 Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
1549 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1550 subB.values(), subB.stride(), S.values(), S.stride(),
1551 ST_ZERO, ritzRes.values(), ritzRes.stride() );
1585 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1586 Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
1589 if (problem_->isHermitian() && !hermImagDetected) {
1592 sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_);
1594 while ( i < curDim_ ) {
1596 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
1597 ritzIndex_.push_back(0);
1604 sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
1609 std::vector<MagnitudeType> ritz2( curDim_ );
1610 for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
1611 blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
1614 ritzValsCurrent_ =
true;
1630 sortSchurForm( *schurH_, *Q_, ritzOrder_ );
1633 schurCurrent_ =
true;
1644 template <
class ScalarType,
class MV,
class OP>
1646 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
1647 std::vector<int>& order )
1650 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1651 Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
1662 int i = 0, nevtemp = 0;
1664 std::vector<int> offset2( curDim_ );
1665 std::vector<int> order2( curDim_ );
1668 Teuchos::LAPACK<int,ScalarType> lapack;
1669 int lwork = 3*curDim_;
1670 std::vector<ScalarType> work( lwork );
1672 while (i < curDim_) {
1673 if ( ritzIndex_[i] != 0 ) {
1674 offset2[nevtemp] = 0;
1675 for (
int j=i; j<curDim_; j++) {
1676 if (order[j] > order[i]) { offset2[nevtemp]++; }
1678 order2[nevtemp] = order[i];
1681 offset2[nevtemp] = 0;
1682 for (
int j=i; j<curDim_; j++) {
1683 if (order[j] > order[i]) { offset2[nevtemp]++; }
1685 order2[nevtemp] = order[i];
1690 ScalarType *ptr_h = H.values();
1691 ScalarType *ptr_q = Q.values();
1692 int ldh = H.stride(), ldq = Q.stride();
1694 for (i=nevtemp-1; i>=0; i--) {
1695 lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i],
1696 1, &work[0], &info );
1697 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1698 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1704 template <
class ScalarType,
class MV,
class OP>
1709 os.setf(std::ios::scientific, std::ios::floatfield);
1711 os <<
"================================================================================" << endl;
1713 os <<
" BlockKrylovSchur Solver Status" << endl;
1715 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
1716 os <<
"The number of iterations performed is " <<iter_<<endl;
1717 os <<
"The block size is " << blockSize_<<endl;
1718 os <<
"The number of blocks is " << numBlocks_<<endl;
1719 os <<
"The current basis size is " << curDim_<<endl;
1720 os <<
"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1721 os <<
"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1723 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1727 os <<
"CURRENT RITZ VALUES "<<endl;
1728 if (ritzIndex_.size() != 0) {
1729 int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
1730 if (problem_->isHermitian()) {
1731 os << std::setw(20) <<
"Ritz Value" 1732 << std::setw(20) <<
"Ritz Residual" 1734 os <<
"--------------------------------------------------------------------------------"<<endl;
1735 for (
int i=0; i<numPrint; i++) {
1736 os << std::setw(20) << ritzValues_[i].realpart
1737 << std::setw(20) << ritzResiduals_[i]
1741 os << std::setw(24) <<
"Ritz Value" 1742 << std::setw(30) <<
"Ritz Residual" 1744 os <<
"--------------------------------------------------------------------------------"<<endl;
1745 for (
int i=0; i<numPrint; i++) {
1747 os << std::setw(15) << ritzValues_[i].realpart;
1748 if (ritzValues_[i].imagpart < MT_ZERO) {
1749 os <<
" - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
1751 os <<
" + i" << std::setw(15) << ritzValues_[i].imagpart;
1753 os << std::setw(20) << ritzResiduals_[i] << endl;
1757 os << std::setw(20) <<
"[ NONE COMPUTED ]" << endl;
1761 os <<
"================================================================================" << endl;
void computeRitzValues()
Compute the Ritz values using the current Krylov factorization.
BlockKrylovSchurOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
void setStepSize(int stepSize)
Set the step size.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the current Ritz residual 2-norms.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors.
This class defines the interface required by an eigensolver and status test class to compute solution...
bool isRitzVecsCurrent() const
Get the status of the Ritz vectors currently stored in the eigensolver.
virtual ~BlockKrylovSchur()
BlockKrylovSchur destructor.
Declaration of basic traits for the multivector type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
void iterate()
This method performs Block Krylov-Schur iterations until the status test indicates the need to stop o...
Virtual base class which defines basic traits for the operator type.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
static void scaleRitzVectors(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, Teuchos::SerialDenseMatrix< int, ScalarType > *S)
Helper function for correctly scaling the eigenvectors of the projected eigenproblem.
Structure to contain pointers to BlockKrylovSchur state variables.
An exception class parent to all Anasazi exceptions.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest 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...
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
void computeRitzVectors()
Compute the Ritz vectors using the current Krylov factorization.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
int getNumIters() const
Get the current iteration count.
void setNumRitzVectors(int numRitzVecs)
Set the number of Ritz vectors to compute.
void computeSchurForm(const bool sort=true)
Compute the Schur form of the projected eigenproblem from the current Krylov factorization.
static void sortRitzValues(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rRV, const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, std::vector< Value< ScalarType > > *RV, std::vector< int > *RO, std::vector< int > *RI)
Helper function for correctly storing the Ritz values when the eigenproblem is non-Hermitian.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
BlockKrylovSchur(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< OrthoManager< ScalarType, MV > > &ortho, Teuchos::ParameterList ¶ms)
BlockKrylovSchur constructor with eigenproblem, solver utilities, and parameter list of solver option...
Output managers remove the need for the eigensolver to know any information about the required output...
void resetNumIters()
Reset the iteration count.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
int getStepSize() const
Get the step size.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values.
std::vector< int > getRitzIndex()
Get the Ritz index vector.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
Teuchos::RCP< const MulVec > V
The current Krylov basis.
BlockKrylovSchurInitFailure is thrown when the BlockKrylovSchur solver is unable to generate an initi...
BlockKrylovSchurState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
bool isSchurCurrent() const
Get the status of the Schur form currently stored in the eigensolver.
bool isRitzValsCurrent() const
Get the status of the Ritz values currently stored in the eigensolver.
Types and exceptions used within Anasazi solvers and interfaces.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
int getNumRitzVectors() const
Get the number of Ritz vectors to compute.
static void computeRitzResiduals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, const Teuchos::SerialDenseMatrix< int, ScalarType > &S, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *RR)
Helper function for correctly computing the Ritz residuals of the projected eigenproblem.
void setBlockSize(int blockSize)
Set the blocksize.
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
int curDim
The current dimension of the reduction.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
Common interface of stopping criteria for Anasazi's solvers.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.