46 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
47 #define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
55 #include "AnasaziHelperTraits.hpp"
88 template <
class ScalarType,
class MulVec>
108 H(Teuchos::null),
S(Teuchos::null),
146 template <
class ScalarType,
class MV,
class OP>
286 std::vector<Value<ScalarType> > ret = ritzValues_;
287 ret.resize(ritzIndex_.size());
302 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms() {
303 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
311 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms() {
312 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
320 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms() {
321 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
322 ret.resize(ritzIndex_.size());
346 void setSize(
int blockSize,
int numBlocks);
372 if (!initialized_)
return 0;
377 int getMaxSubspaceDim()
const {
return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
443 typedef typename std::vector<ScalarType>::iterator STiter;
444 typedef typename std::vector<MagnitudeType>::iterator MTiter;
445 const MagnitudeType MT_ONE;
446 const MagnitudeType MT_ZERO;
447 const MagnitudeType NANVAL;
448 const ScalarType ST_ONE;
449 const ScalarType ST_ZERO;
457 CheckList() : checkV(false), checkArn(false), checkAux(false) {};
462 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
465 std::vector<int>& order );
482 timerCompSF_, timerSortSF_,
483 timerCompRitzVec_, timerOrtho_;
539 bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
542 std::vector<Value<ScalarType> > ritzValues_;
543 std::vector<MagnitudeType> ritzResiduals_;
546 std::vector<int> ritzIndex_;
547 std::vector<int> ritzOrder_;
556 template <
class ScalarType,
class MV,
class OP>
565 MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
566 MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
567 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
568 ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
569 ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
577 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
578 timerOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Operation Op*x")),
579 timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Sorting Ritz values")),
580 timerCompSF_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Computing Schur form")),
581 timerSortSF_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Sorting Schur form")),
582 timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
583 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchur::Orthogonalization")),
593 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
596 ritzVecsCurrent_(false),
597 ritzValsCurrent_(false),
598 schurCurrent_(false),
602 "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
604 "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
606 "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
608 "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
610 "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
612 "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
614 "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
615 TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
616 "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
617 TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
618 "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
619 TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
620 "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
623 Op_ = problem_->getOperator();
626 TEUCHOS_TEST_FOR_EXCEPTION(!params.
isParameter(
"Step Size"), std::invalid_argument,
627 "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
628 int ss = params.
get(
"Step Size",numBlocks_);
632 int bs = params.
get(
"Block Size", 1);
633 int nb = params.
get(
"Num Blocks", 3*problem_->getNEV());
638 int numRitzVecs = params.
get(
"Number of Ritz Vectors", 0);
642 numRitzPrint_ = params.
get(
"Print Number of Ritz Values", bs);
649 template <
class ScalarType,
class MV,
class OP>
652 setSize(blockSize,numBlocks_);
658 template <
class ScalarType,
class MV,
class OP>
661 TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
662 stepSize_ = stepSize;
667 template <
class ScalarType,
class MV,
class OP>
673 TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
676 if (numRitzVecs != numRitzVecs_) {
678 ritzVectors_ = Teuchos::null;
679 ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
681 ritzVectors_ = Teuchos::null;
683 numRitzVecs_ = numRitzVecs;
684 ritzVecsCurrent_ =
false;
690 template <
class ScalarType,
class MV,
class OP>
696 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
697 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument,
"Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
698 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
703 blockSize_ = blockSize;
704 numBlocks_ = numBlocks;
712 if (problem_->getInitVec() != Teuchos::null) {
713 tmp = problem_->getInitVec();
718 "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
726 if (problem_->isHermitian()) {
727 newsd = blockSize_*numBlocks_;
729 newsd = blockSize_*numBlocks_+1;
733 "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
735 ritzValues_.resize(newsd);
736 ritzResiduals_.resize(newsd,MT_ONE);
737 ritzOrder_.resize(newsd);
739 V_ = MVT::Clone(*tmp,newsd+blockSize_);
743 initialized_ =
false;
750 template <
class ScalarType,
class MV,
class OP>
757 if (om_->isVerbosity(
Debug ) ) {
761 om_->print(
Debug, accuracyCheck(chk,
": in setAuxVecs()") );
765 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
766 numAuxVecs_ += MVT::GetNumberVecs(**i);
770 if (numAuxVecs_ > 0 && initialized_) {
771 initialized_ =
false;
784 template <
class ScalarType,
class MV,
class OP>
789 std::vector<int> bsind(blockSize_);
790 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
798 std::string errstr(
"Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
802 if (newstate.
V != Teuchos::null && newstate.
H != Teuchos::null) {
807 std::invalid_argument, errstr );
808 if (newstate.
V != V_) {
810 std::invalid_argument, errstr );
812 std::invalid_argument, errstr );
815 std::invalid_argument, errstr );
817 curDim_ = newstate.
curDim;
818 int lclDim = MVT::GetNumberVecs(*newstate.
V);
823 if (curDim_ == 0 && lclDim > blockSize_) {
824 om_->stream(
Warnings) <<
"Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
825 <<
"The block size however is only " << blockSize_ << std::endl
826 <<
"The last " << lclDim - blockSize_ <<
" vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
831 if (newstate.
V != V_) {
832 std::vector<int> nevind(lclDim);
833 for (
int i=0; i<lclDim; i++) nevind[i] = i;
834 MVT::SetBlock(*newstate.
V,nevind,*V_);
838 if (newstate.
H != H_) {
839 H_->putScalar( ST_ZERO );
846 lclH = Teuchos::null;
855 "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
857 int lclDim = MVT::GetNumberVecs(*ivec);
858 bool userand =
false;
859 if (lclDim < blockSize_) {
867 std::vector<int> dimind2(lclDim);
868 for (
int i=0; i<lclDim; i++) { dimind2[i] = i; }
874 MVT::SetBlock(*ivec,dimind2,*newV1);
877 dimind2.resize(blockSize_-lclDim);
878 for (
int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
882 MVT::MvRandom(*newV2);
892 MVT::SetBlock(*ivecV,bsind,*newV1);
899 if (auxVecs_.size() > 0) {
900 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
905 int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
907 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
910 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
914 int rank = orthman_->normalize(*newV);
916 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
923 newV = Teuchos::null;
927 ritzVecsCurrent_ =
false;
928 ritzValsCurrent_ =
false;
929 schurCurrent_ =
false;
934 if (om_->isVerbosity(
Debug ) ) {
940 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
944 if (om_->isVerbosity(
Debug)) {
945 currentStatus( om_->stream(
Debug) );
955 template <
class ScalarType,
class MV,
class OP>
965 template <
class ScalarType,
class MV,
class OP>
971 if (initialized_ ==
false) {
977 int searchDim = blockSize_*numBlocks_;
978 if (problem_->isHermitian() ==
false) {
987 while (tester_->checkStatus(
this) !=
Passed && curDim_+blockSize_ <= searchDim) {
992 int lclDim = curDim_ + blockSize_;
995 std::vector<int> curind(blockSize_);
996 for (
int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
1001 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1008 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1011 OPT::Apply(*Op_,*Vprev,*Vnext);
1012 count_ApplyOp_ += blockSize_;
1016 Vprev = Teuchos::null;
1020 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1025 std::vector<int> prevind(lclDim);
1026 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
1027 Vprev = MVT::CloneView(*V_,prevind);
1038 if (auxVecs_.size() > 0) {
1039 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1040 AVprev.
append( auxVecs_[i] );
1041 AsubH.
append( Teuchos::null );
1050 (
Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
1051 int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
1058 "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
1064 Vnext = Teuchos::null;
1065 curDim_ += blockSize_;
1067 ritzVecsCurrent_ =
false;
1068 ritzValsCurrent_ =
false;
1069 schurCurrent_ =
false;
1072 if (!(iter_%stepSize_)) {
1073 computeRitzValues();
1077 if (om_->isVerbosity(
Debug ) ) {
1081 chk.checkArn =
true;
1082 om_->print(
Debug, accuracyCheck(chk,
": after local update") );
1087 om_->print(
OrthoDetails, accuracyCheck(chk,
": after local update") );
1091 if (om_->isVerbosity(
Debug)) {
1092 currentStatus( om_->stream(
Debug) );
1116 template <
class ScalarType,
class MV,
class OP>
1119 std::stringstream os;
1121 os.setf(std::ios::scientific, std::ios::floatfield);
1124 os <<
" Debugging checks: iteration " << iter_ << where << std::endl;
1127 std::vector<int> lclind(curDim_);
1128 for (
int i=0; i<curDim_; i++) lclind[i] = i;
1129 std::vector<int> bsind(blockSize_);
1130 for (
int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
1135 lclV = MVT::CloneView(*V_,lclind);
1136 lclF = MVT::CloneView(*V_,bsind);
1140 tmp = orthman_->orthonormError(*lclV);
1141 os <<
" >> Error in V^H M V == I : " << tmp << std::endl;
1143 tmp = orthman_->orthonormError(*lclF);
1144 os <<
" >> Error in F^H M F == I : " << tmp << std::endl;
1146 tmp = orthman_->orthogError(*lclV,*lclF);
1147 os <<
" >> Error in V^H M F == 0 : " << tmp << std::endl;
1149 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1151 tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
1152 os <<
" >> Error in V^H M Aux[" << i <<
"] == 0 : " << tmp << std::endl;
1154 tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
1155 os <<
" >> Error in F^H M Aux[" << i <<
"] == 0 : " << tmp << std::endl;
1163 lclAV = MVT::Clone(*V_,curDim_);
1165 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1168 OPT::Apply(*Op_,*lclV,*lclAV);
1173 MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
1177 blockSize_,curDim_, curDim_ );
1178 MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
1181 std::vector<MagnitudeType> arnNorms( curDim_ );
1182 orthman_->norm( *lclAV, arnNorms );
1184 for (
int i=0; i<curDim_; i++) {
1185 os <<
" >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i <<
"]|| : " << arnNorms[i] << std::endl;
1191 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1192 tmp = orthman_->orthonormError(*auxVecs_[i]);
1193 os <<
" >> Error in Aux[" << i <<
"]^H M Aux[" << i <<
"] == I : " << tmp << std::endl;
1194 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1195 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1196 os <<
" >> Error in Aux[" << i <<
"]^H M Aux[" << j <<
"] == 0 : " << tmp << std::endl;
1217 template <
class ScalarType,
class MV,
class OP>
1225 if (!ritzValsCurrent_) {
1228 computeSchurForm(
false );
1245 template <
class ScalarType,
class MV,
class OP>
1248 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1253 "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
1256 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
1260 if (curDim_ && initialized_) {
1263 if (!ritzVecsCurrent_) {
1266 if (!schurCurrent_) {
1269 computeSchurForm(
true );
1275 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1288 std::vector<int> curind( curDim_ );
1289 for (
int i=0; i<curDim_; i++) { curind[i] = i; }
1291 if (problem_->isHermitian()) {
1296 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
1299 bool complexRitzVals =
false;
1300 for (
int i=0; i<numRitzVecs_; i++) {
1301 if (ritzIndex_[i] != 0) {
1302 complexRitzVals =
true;
1306 if (complexRitzVals)
1307 om_->stream(
Warnings) <<
" Eigenproblem is Hermitian and complex eigenvalues have converged, corresponding eigenvectors will be incorrect!!!"
1318 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
1326 int lwork = 3*curDim_;
1327 std::vector<ScalarType> work( lwork );
1328 std::vector<MagnitudeType> rwork( curDim_ );
1332 ScalarType vl[ ldvl ];
1334 lapack.
TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1335 copyQ.
values(), copyQ.
stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1337 "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1343 curind.resize( numRitzVecs_ );
1344 Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
1345 MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
1348 std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
1349 MVT::MvNorm( *view_ritzVectors, ritzNrm );
1352 tmpritzVectors_ = Teuchos::null;
1353 view_ritzVectors = Teuchos::null;
1356 ScalarType ritzScale = ST_ONE;
1357 for (
int i=0; i<numRitzVecs_; i++) {
1360 if (ritzIndex_[i] == 1 ) {
1361 ritzScale = ST_ONE/lapack_mag.
LAPY2(ritzNrm[i],ritzNrm[i+1]);
1362 std::vector<int> newind(2);
1363 newind[0] = i; newind[1] = i+1;
1364 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1365 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1366 MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1373 std::vector<int> newind(1);
1375 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1376 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1377 MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1384 ritzVecsCurrent_ =
true;
1393 template <
class ScalarType,
class MV,
class OP>
1396 "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
1402 template <
class ScalarType,
class MV,
class OP>
1420 template <
class ScalarType,
class MV,
class OP>
1424 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1432 if (!schurCurrent_) {
1436 if (!ritzValsCurrent_) {
1460 int lwork = 3*curDim_;
1461 std::vector<ScalarType> work( lwork );
1462 std::vector<MagnitudeType> rwork( curDim_ );
1463 std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
1464 std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
1465 std::vector<int> bwork( curDim_ );
1466 int info = 0, sdim = 0;
1468 lapack.
GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
1469 &tmp_iRitzValues[0], subQ.
values(), subQ.
stride(), &work[0], lwork,
1470 &rwork[0], &bwork[0], &info );
1473 "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1479 bool hermImagDetected =
false;
1480 if (problem_->isHermitian()) {
1481 for (
int i=0; i<curDim_; i++)
1483 if (tmp_iRitzValues[i] != MT_ZERO)
1485 hermImagDetected =
true;
1489 if (hermImagDetected) {
1491 om_->stream(
Warnings) <<
" Eigenproblem is Hermitian and complex eigenvalues have been detected!!!" << std::endl;
1499 (*tLocalH) -= localH;
1500 MagnitudeType normF = tLocalH->normFrobenius();
1501 MagnitudeType norm1 = tLocalH->normOne();
1502 om_->stream(
Warnings) <<
" Symmetry error in projected eigenproblem: || S - S' ||_F = " << normF
1503 <<
", || S - S' ||_1 = " << norm1 << std::endl;
1522 blockSize_, curDim_, curDim_ );
1532 ScalarType* b_ptr = subB.
values();
1533 if (problem_->isHermitian() && !hermImagDetected) {
1537 for (
int i=0; i<curDim_ ; i++) {
1538 ritzResiduals_[i] = blas.
NRM2(blockSize_, b_ptr + i*blockSize_, 1);
1547 ScalarType vl[ ldvl ];
1549 lapack.
TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1550 S.
values(), S.
stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1553 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1598 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1602 if (problem_->isHermitian() && !hermImagDetected) {
1605 sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_);
1607 while ( i < curDim_ ) {
1609 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
1610 ritzIndex_.push_back(0);
1617 sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
1622 std::vector<MagnitudeType> ritz2( curDim_ );
1623 for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
1624 blas_mag.
COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
1627 ritzValsCurrent_ =
true;
1643 sortSchurForm( *schurH_, *Q_, ritzOrder_ );
1646 schurCurrent_ =
true;
1657 template <
class ScalarType,
class MV,
class OP>
1660 std::vector<int>& order )
1663 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1675 int i = 0, nevtemp = 0;
1677 std::vector<int> offset2( curDim_ );
1678 std::vector<int> order2( curDim_ );
1682 int lwork = 3*curDim_;
1683 std::vector<ScalarType> work( lwork );
1685 while (i < curDim_) {
1686 if ( ritzIndex_[i] != 0 ) {
1687 offset2[nevtemp] = 0;
1688 for (
int j=i; j<curDim_; j++) {
1689 if (order[j] > order[i]) { offset2[nevtemp]++; }
1691 order2[nevtemp] = order[i];
1694 offset2[nevtemp] = 0;
1695 for (
int j=i; j<curDim_; j++) {
1696 if (order[j] > order[i]) { offset2[nevtemp]++; }
1698 order2[nevtemp] = order[i];
1703 ScalarType *ptr_h = H.
values();
1704 ScalarType *ptr_q = Q.
values();
1707 for (i=nevtemp-1; i>=0; i--) {
1708 int ifst = order2[i]+1+offset2[i];
1710 lapack.
TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, &ifst,
1711 &ilst, &work[0], &info );
1713 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1719 template <
class ScalarType,
class MV,
class OP>
1724 os.setf(std::ios::scientific, std::ios::floatfield);
1726 os <<
"================================================================================" << endl;
1728 os <<
" BlockKrylovSchur Solver Status" << endl;
1730 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
1731 os <<
"The number of iterations performed is " <<iter_<<endl;
1732 os <<
"The block size is " << blockSize_<<endl;
1733 os <<
"The number of blocks is " << numBlocks_<<endl;
1734 os <<
"The current basis size is " << curDim_<<endl;
1735 os <<
"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1736 os <<
"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1738 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1742 os <<
"CURRENT RITZ VALUES "<<endl;
1743 if (ritzIndex_.size() != 0) {
1744 int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
1745 if (problem_->isHermitian()) {
1746 os << std::setw(20) <<
"Ritz Value"
1747 << std::setw(20) <<
"Ritz Residual"
1749 os <<
"--------------------------------------------------------------------------------"<<endl;
1750 for (
int i=0; i<numPrint; i++) {
1751 os << std::setw(20) << ritzValues_[i].realpart
1752 << std::setw(20) << ritzResiduals_[i]
1756 os << std::setw(24) <<
"Ritz Value"
1757 << std::setw(30) <<
"Ritz Residual"
1759 os <<
"--------------------------------------------------------------------------------"<<endl;
1760 for (
int i=0; i<numPrint; i++) {
1762 os << std::setw(15) << ritzValues_[i].realpart;
1763 if (ritzValues_[i].imagpart < MT_ZERO) {
1766 os <<
" + i" << std::setw(15) << ritzValues_[i].imagpart;
1768 os << std::setw(20) << ritzResiduals_[i] << endl;
1772 os << std::setw(20) <<
"[ NONE COMPUTED ]" << endl;
1776 os <<
"================================================================================" << endl;
ScalarType * values() const
void computeRitzValues()
Compute the Ritz values using the current Krylov factorization.
BlockKrylovSchurOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
BlockKrylovSchurState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
void setStepSize(int stepSize)
Set the step size.
bool isRitzValsCurrent() const
Get the status of the Ritz values currently stored in the eigensolver.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the current Ritz residual 2-norms.
Array< T > & append(const T &x)
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...
virtual ~BlockKrylovSchur()
BlockKrylovSchur destructor.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Declaration of basic traits for the multivector type.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
T & get(const std::string &name, T def_value)
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
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.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
An exception class parent to all Anasazi exceptions.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
void TREVC(const char &SIDE, const char &HOWMNY, OrdinalType *select, const OrdinalType &n, const ScalarType *T, const OrdinalType &ldt, ScalarType *VL, const OrdinalType &ldvl, ScalarType *VR, const OrdinalType &ldvr, const OrdinalType &mm, OrdinalType *m, ScalarType *WORK, OrdinalType *info) const
int getNumIters() const
Get the current iteration count.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
bool isRitzVecsCurrent() const
Get the status of the Ritz vectors currently stored in the eigensolver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
TEUCHOS_DEPRECATED void TREXC(const char &COMPQ, const OrdinalType &n, ScalarType *T, const OrdinalType &ldt, ScalarType *Q, const OrdinalType &ldq, const OrdinalType &ifst, const OrdinalType &ilst, ScalarType *WORK, OrdinalType *info) const
void computeRitzVectors()
Compute the Ritz vectors using the current Krylov factorization.
int getStepSize() const
Get the step size.
bool isParameter(const std::string &name) const
void setNumRitzVectors(int numRitzVecs)
Set the number of Ritz vectors to compute.
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSchurCurrent() const
Get the status of the Schur form currently stored in the eigensolver.
void GEES(const char &JOBVS, const char &SORT, OrdinalType &(*ptr2func)(ScalarType *, ScalarType *), const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *sdim, ScalarType *WR, ScalarType *WI, ScalarType *VS, const OrdinalType &ldvs, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *BWORK, OrdinalType *info) const
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 getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getNumRitzVectors() const
Get the number of Ritz vectors to compute.
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values.
std::vector< int > getRitzIndex()
Get the Ritz index vector.
static magnitudeType magnitude(T a)
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...
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...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
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.
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.
OrdinalType stride() const