33 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_HPP 34 #define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP 42 #include "AnasaziHelperTraits.hpp" 75 template <
class ScalarType,
class MulVec>
133 template <
class ScalarType,
class MV,
class OP>
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());
333 void setSize(
int blockSize,
int numBlocks);
359 if (!initialized_)
return 0;
364 int getMaxSubspaceDim()
const {
return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
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;
452 std::vector<int>& order );
469 timerCompSF_, timerSortSF_,
470 timerCompRitzVec_, timerOrtho_;
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>
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")),
583 ritzVecsCurrent_(false),
584 ritzValsCurrent_(false),
585 schurCurrent_(false),
589 "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
591 "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
593 "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
595 "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
597 "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
599 "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
601 "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
603 "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
605 "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
607 "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
610 Op_ = problem_->getOperator();
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;
699 if (problem_->getInitVec() != Teuchos::null) {
700 tmp = problem_->getInitVec();
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;
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_);
730 initialized_ =
false;
737 template <
class ScalarType,
class MV,
class OP>
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) {
794 std::invalid_argument, errstr );
795 if (newstate.
V != V_) {
797 std::invalid_argument, errstr );
799 std::invalid_argument, errstr );
802 std::invalid_argument, errstr );
804 curDim_ = newstate.
curDim;
805 int lclDim = MVT::GetNumberVecs(*newstate.
V);
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 );
833 lclH = Teuchos::null;
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; }
861 MVT::SetBlock(*ivec,dimind2,*newV1);
864 dimind2.resize(blockSize_-lclDim);
865 for (
int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
869 MVT::MvRandom(*newV2);
879 MVT::SetBlock(*ivecV,bsind,*newV1);
886 if (auxVecs_.size() > 0) {
887 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 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; }
988 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
995 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 998 OPT::Apply(*Op_,*Vprev,*Vnext);
999 count_ApplyOp_ += blockSize_;
1003 Vprev = Teuchos::null;
1007 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1012 std::vector<int> prevind(lclDim);
1013 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
1014 Vprev = MVT::CloneView(*V_,prevind);
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 );
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; }
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 1155 OPT::Apply(*Op_,*lclV,*lclAV);
1160 MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
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 1240 "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
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 );
1262 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1275 std::vector<int> curind( curDim_ );
1276 for (
int i=0; i<curDim_; i++) { curind[i] = i; }
1278 if (problem_->isHermitian()) {
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!!!" 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 ];
1321 lapack.
TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1322 copyQ.
values(), copyQ.
stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1324 "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
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>
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 1419 if (!schurCurrent_) {
1423 if (!ritzValsCurrent_) {
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 );
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;
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;
1509 blockSize_, curDim_, curDim_ );
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 ];
1536 lapack.
TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1537 S.
values(), S.
stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1540 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ <<
") returned info " << info <<
" != 0.");
1585 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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>
1647 std::vector<int>& order )
1650 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1662 int i = 0, nevtemp = 0;
1664 std::vector<int> offset2( curDim_ );
1665 std::vector<int> order2( curDim_ );
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();
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 );
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) {
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.
ScalarType LAPY2(const ScalarType x, const ScalarType y) const
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the current Ritz residual 2-norms.
Array< T > & append(const T &x)
ScalarType * values() const
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.
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.
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
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
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.
void TREXC(const char COMPQ, const OrdinalType n, ScalarType *T, const OrdinalType ldt, ScalarType *Q, const OrdinalType ldq, OrdinalType ifst, OrdinalType ilst, ScalarType *WORK, OrdinalType *info) const
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
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
static magnitudeType magnitude(T a)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void COPY(const OrdinalType n, const ScalarType *x, const OrdinalType incx, ScalarType *y, const OrdinalType incy) const
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.
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
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.
bool isParameter(const std::string &name) const
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...
OrdinalType stride() const
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.