46#ifndef ANASAZI_RTRBASE_HPP
47#define ANASAZI_RTRBASE_HPP
54#include "Teuchos_ScalarTraits.hpp"
59#include "Teuchos_LAPACK.hpp"
60#include "Teuchos_BLAS.hpp"
61#include "Teuchos_SerialDenseMatrix.hpp"
62#include "Teuchos_ParameterList.hpp"
63#include "Teuchos_TimeMonitor.hpp"
123 template <
class ScalarType,
class MV>
126 Teuchos::RCP<const MV>
X;
128 Teuchos::RCP<const MV>
AX;
130 Teuchos::RCP<const MV>
BX;
132 Teuchos::RCP<const MV>
R;
134 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >
T;
138 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
rho;
139 RTRState() :
X(Teuchos::null),
AX(Teuchos::null),
BX(Teuchos::null),
140 R(Teuchos::null),
T(Teuchos::null),
rho(0) {};
201 template <
class ScalarType,
class MV,
class OP>
214 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
218 Teuchos::ParameterList ¶ms,
219 const std::string &solverLabel,
bool skinnySolver
348 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
356 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
363 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
389 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
430 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
433 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const;
452 typedef Teuchos::ScalarTraits<ScalarType> SCT;
453 typedef typename SCT::magnitudeType MagnitudeType;
454 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
455 const MagnitudeType ONE;
456 const MagnitudeType ZERO;
457 const MagnitudeType NANVAL;
458 typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
459 typedef typename std::vector<ScalarType>::iterator vecSTiter;
464 bool checkX, checkAX, checkBX;
465 bool checkEta, checkAEta, checkBEta;
466 bool checkR, checkQ, checkBR;
467 bool checkZ, checkPBX;
468 CheckList() : checkX(false),checkAX(false),checkBX(false),
469 checkEta(false),checkAEta(false),checkBEta(false),
470 checkR(false),checkQ(false),checkBR(false),
471 checkZ(false), checkPBX(false) {};
476 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
478 virtual void solveTRSubproblem() = 0;
480 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(
const MV &xi)
const;
481 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(
const MV &xi,
const MV &zeta)
const;
482 void ginnersep(
const MV &xi, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d)
const;
483 void ginnersep(
const MV &xi,
const MV &zeta, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d)
const;
487 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
488 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
489 const Teuchos::RCP<OutputManager<ScalarType> > om_;
490 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
491 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_;
495 Teuchos::RCP<const OP> AOp_;
496 Teuchos::RCP<const OP> BOp_;
497 Teuchos::RCP<const OP> Prec_;
498 bool hasBOp_, hasPrec_, olsenPrec_;
502 Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
504 timerLocalProj_, timerDS_,
505 timerLocalUpdate_, timerCompRes_,
506 timerOrtho_, timerInit_;
511 int counterAOp_, counterBOp_, counterPrec_;
574 Teuchos::RCP<MV> V_, BV_, PBV_,
577 delta_, Adelta_, Bdelta_, Hdelta_,
579 Teuchos::RCP<const MV> X_, BX_;
582 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
589 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
592 bool Rnorms_current_, R2norms_current_;
595 MagnitudeType conv_kappa_, conv_theta_;
607 template <
class ScalarType,
class MV,
class OP>
610 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
614 Teuchos::ParameterList ¶ms,
615 const std::string &solverLabel,
bool skinnySolver
617 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
618 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
619 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
626#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
627 timerAOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation A*x")),
628 timerBOp_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation B*x")),
629 timerPrec_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Operation Prec*x")),
630 timerSort_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Sorting eigenvalues")),
631 timerLocalProj_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Local projection")),
632 timerDS_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Direct solve")),
633 timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Local update")),
634 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Computing residuals")),
635 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Orthogonalization")),
636 timerInit_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: "+solverLabel+
"::Initialization")),
645 isSkinny_(skinnySolver),
647 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
650 Rnorms_current_(false),
651 R2norms_current_(false),
657 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
658 "Anasazi::RTRBase::constructor: user passed null problem pointer.");
659 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
660 "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
661 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
662 "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
663 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
664 "Anasazi::RTRBase::constructor: user passed null status test pointer.");
665 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
666 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
667 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() ==
false, std::invalid_argument,
668 "Anasazi::RTRBase::constructor: problem is not set.");
669 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
670 "Anasazi::RTRBase::constructor: problem is not Hermitian.");
673 AOp_ = problem_->getOperator();
674 TEUCHOS_TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
675 "Anasazi::RTRBase::constructor: problem provides no A matrix.");
676 BOp_ = problem_->getM();
677 Prec_ = problem_->getPrec();
678 hasBOp_ = (BOp_ != Teuchos::null);
679 hasPrec_ = (Prec_ != Teuchos::null);
680 olsenPrec_ = params.get<
bool>(
"Olsen Prec",
true);
682 TEUCHOS_TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
683 "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
686 int bs = params.get(
"Block Size", problem_->getNEV());
690 leftMost_ = params.get(
"Leftmost",leftMost_);
692 conv_kappa_ = params.get(
"Kappa Convergence",conv_kappa_);
693 TEUCHOS_TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
694 "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
695 conv_theta_ = params.get(
"Theta Convergence",conv_theta_);
696 TEUCHOS_TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
697 "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
703 template <
class ScalarType,
class MV,
class OP>
707#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
708 Teuchos::TimeMonitor lcltimer( *timerInit_ );
715 Teuchos::RCP<const MV> tmp;
721 if (blockSize_ > 0) {
725 tmp = problem_->getInitVec();
726 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
727 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
730 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetGlobalLength(*tmp), std::invalid_argument,
731 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
734 if (blockSize == blockSize_) {
753 if (blockSize > blockSize_)
757 Teuchos::RCP<const MV> Q;
758 std::vector<int> indQ(numAuxVecs_);
759 for (
int i=0; i<numAuxVecs_; i++) indQ[i] = i;
761 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
762 "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
764 if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
765 V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
766 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
769 if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
770 BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
771 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
777 if (hasPrec_ && olsenPrec_) {
778 if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
779 PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
780 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
786 std::vector<int> indV(numAuxVecs_+blockSize);
787 for (
int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
789 V_ = MVT::CloneCopy(*V_,indV);
792 BV_ = MVT::CloneCopy(*BV_,indV);
798 if (hasPrec_ && olsenPrec_) {
799 PBV_ = MVT::CloneCopy(*PBV_,indV);
803 if (blockSize < blockSize_) {
805 blockSize_ = blockSize;
807 theta_.resize(blockSize_);
808 ritz2norms_.resize(blockSize_);
809 Rnorms_.resize(blockSize_);
810 R2norms_.resize(blockSize_);
814 std::vector<int> ind(blockSize_);
815 for (
int i=0; i<blockSize_; i++) ind[i] = i;
823 R_ = MVT::CloneCopy(*R_ ,ind);
824 eta_ = MVT::CloneCopy(*eta_ ,ind);
825 delta_ = MVT::CloneCopy(*delta_ ,ind);
826 Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
828 AX_ = MVT::CloneCopy(*AX_ ,ind);
829 Aeta_ = MVT::CloneCopy(*Aeta_ ,ind);
830 Adelta_ = MVT::CloneCopy(*Adelta_,ind);
834 Aeta_ = Teuchos::null;
835 Adelta_ = Teuchos::null;
840 Beta_ = MVT::CloneCopy(*Beta_,ind);
841 Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
844 Beta_ = Teuchos::null;
845 Bdelta_ = Teuchos::null;
855 if (hasPrec_ || isSkinny_) {
856 Z_ = MVT::Clone(*V_,blockSize_);
868 eta_ = Teuchos::null;
869 Aeta_ = Teuchos::null;
870 delta_ = Teuchos::null;
871 Adelta_ = Teuchos::null;
872 Hdelta_ = Teuchos::null;
873 Beta_ = Teuchos::null;
874 Bdelta_ = Teuchos::null;
877 R_ = MVT::Clone(*tmp,blockSize_);
878 eta_ = MVT::Clone(*tmp,blockSize_);
879 delta_ = MVT::Clone(*tmp,blockSize_);
880 Hdelta_ = MVT::Clone(*tmp,blockSize_);
882 AX_ = MVT::Clone(*tmp,blockSize_);
883 Aeta_ = MVT::Clone(*tmp,blockSize_);
884 Adelta_ = MVT::Clone(*tmp,blockSize_);
889 Beta_ = MVT::Clone(*tmp,blockSize_);
890 Bdelta_ = MVT::Clone(*tmp,blockSize_);
900 if (hasPrec_ || isSkinny_) {
901 Z_ = MVT::Clone(*tmp,blockSize_);
910 initialized_ =
false;
913 theta_.resize(blockSize,NANVAL);
914 ritz2norms_.resize(blockSize,NANVAL);
915 Rnorms_.resize(blockSize,NANVAL);
916 R2norms_.resize(blockSize,NANVAL);
921 eta_ = Teuchos::null;
922 Aeta_ = Teuchos::null;
923 delta_ = Teuchos::null;
924 Adelta_ = Teuchos::null;
925 Hdelta_ = Teuchos::null;
926 Beta_ = Teuchos::null;
927 Bdelta_ = Teuchos::null;
931 R_ = MVT::Clone(*tmp,blockSize);
932 eta_ = MVT::Clone(*tmp,blockSize);
933 delta_ = MVT::Clone(*tmp,blockSize);
934 Hdelta_ = MVT::Clone(*tmp,blockSize);
936 AX_ = MVT::Clone(*tmp,blockSize);
937 Aeta_ = MVT::Clone(*tmp,blockSize);
938 Adelta_ = MVT::Clone(*tmp,blockSize);
943 Beta_ = MVT::Clone(*tmp,blockSize);
944 Bdelta_ = MVT::Clone(*tmp,blockSize);
951 if (hasPrec_ || isSkinny_) {
952 Z_ = MVT::Clone(*tmp,blockSize);
957 blockSize_ = blockSize;
963 std::vector<int> indX(blockSize_);
964 for (
int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
965 X_ = MVT::CloneView(*V_,indX);
967 BX_ = MVT::CloneView(*BV_,indX);
978 template <
class ScalarType,
class MV,
class OP>
980 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
981 "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
988 template <
class ScalarType,
class MV,
class OP>
996 template <
class ScalarType,
class MV,
class OP>
998 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
1002 auxVecs_.reserve(auxvecs.size());
1005 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1006 numAuxVecs_ += MVT::GetNumberVecs(**v);
1010 if (numAuxVecs_ > 0 && initialized_) {
1011 initialized_ =
false;
1016 BX_ = Teuchos::null;
1019 BV_ = Teuchos::null;
1020 PBV_ = Teuchos::null;
1023 if (numAuxVecs_ > 0) {
1024 V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1026 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1027 std::vector<int> ind(MVT::GetNumberVecs(**v));
1028 for (
size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
1029 MVT::SetBlock(**v,ind,*V_);
1030 auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
1032 TEUCHOS_TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
1033 "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
1036 BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1037 OPT::Apply(*BOp_,*V_,*BV_);
1042 if (hasPrec_ && olsenPrec_) {
1043 PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1044 OPT::Apply(*Prec_,*BV_,*V_);
1049 if (om_->isVerbosity(
Debug ) ) {
1053 om_->print(
Debug, accuracyCheck(chk,
"in setAuxVecs()") );
1070 template <
class ScalarType,
class MV,
class OP>
1079 BX_ = Teuchos::null;
1080 Teuchos::RCP<MV> X, BX;
1082 std::vector<int> ind(blockSize_);
1083 for (
int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1084 X = MVT::CloneViewNonConst(*V_,ind);
1086 BX = MVT::CloneViewNonConst(*BV_,ind);
1093#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1094 Teuchos::TimeMonitor inittimer( *timerInit_ );
1097 std::vector<int> bsind(blockSize_);
1098 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
1117 if (newstate.
X != Teuchos::null) {
1118 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X),
1119 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
1121 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
X) < blockSize_,
1122 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
1125 MVT::SetBlock(*newstate.
X,bsind,*X);
1134 if (newstate.
AX != Teuchos::null) {
1135 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
AX) != MVT::GetGlobalLength(*X),
1136 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
1138 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
AX) < blockSize_,
1139 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
1140 MVT::SetBlock(*newstate.
AX,bsind,*AX_);
1144#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1145 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1147 OPT::Apply(*AOp_,*X,*AX_);
1148 counterAOp_ += blockSize_;
1151 newstate.
R = Teuchos::null;
1157 if (newstate.
BX != Teuchos::null) {
1158 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
BX) != MVT::GetGlobalLength(*X),
1159 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
1161 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
BX) < blockSize_,
1162 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
1163 MVT::SetBlock(*newstate.
BX,bsind,*BX);
1167#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1168 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1170 OPT::Apply(*BOp_,*X,*BX);
1171 counterBOp_ += blockSize_;
1174 newstate.
R = Teuchos::null;
1179 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error,
"Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1187 newstate.
R = Teuchos::null;
1188 newstate.
T = Teuchos::null;
1191 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1192 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1193 "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1195 int initSize = MVT::GetNumberVecs(*ivec);
1196 if (initSize > blockSize_) {
1198 initSize = blockSize_;
1199 std::vector<int> ind(blockSize_);
1200 for (
int i=0; i<blockSize_; i++) ind[i] = i;
1201 ivec = MVT::CloneView(*ivec,ind);
1206 std::vector<int> ind(initSize);
1207 for (
int i=0; i<initSize; i++) ind[i] = i;
1208 MVT::SetBlock(*ivec,ind,*X);
1211 if (blockSize_ > initSize) {
1212 std::vector<int> ind(blockSize_ - initSize);
1213 for (
int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1214 Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind);
1221#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1224 OPT::Apply(*BOp_,*X,*BX);
1225 counterBOp_ += blockSize_;
1229 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error,
"Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1233 if (numAuxVecs_ > 0) {
1234#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1235 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1237 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1238 int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
1240 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1243#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1244 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1246 int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
1248 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1256#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1257 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1259 OPT::Apply(*AOp_,*X,*AX_);
1260 counterAOp_ += blockSize_;
1267 if (newstate.
T != Teuchos::null) {
1268 TEUCHOS_TEST_FOR_EXCEPTION( (
signed int)(newstate.
T->size()) < blockSize_,
1269 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1270 for (
int i=0; i<blockSize_; i++) {
1271 theta_[i] = (*newstate.
T)[i];
1276 Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
1277 BB(blockSize_,blockSize_),
1278 S(blockSize_,blockSize_);
1281#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1282 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1284 MVT::MvTransMv(ONE,*X,*AX_,AA);
1286 MVT::MvTransMv(ONE,*X,*BX,BB);
1289 nevLocal_ = blockSize_;
1295#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1296 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1298 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
1301#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1302 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1304 ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
1307 "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
1308 TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,
RTRInitFailure,
"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
1313#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1314 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1317 std::vector<int> order(blockSize_);
1320 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
1323 Utils::permuteVectors(order,S);
1328 Teuchos::BLAS<int,ScalarType> blas;
1329 Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
1333 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
1334 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1339 for (
int i=0; i<blockSize_; i++) {
1340 blas.SCAL(blockSize_,theta_[i],RR[i],1);
1342 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
1343 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1344 for (
int i=0; i<blockSize_; i++) {
1345 ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
1353#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1354 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1357 MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );
1358 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
1360 MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );
1361 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
1364 MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );
1365 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
1374 std::vector<int> ind(blockSize_);
1375 for (
int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1376 this->X_ = MVT::CloneView(
static_cast<const MV&
>(*V_),ind);
1378 this->BX_ = MVT::CloneView(
static_cast<const MV&
>(*BV_),ind);
1381 this->BX_ = this->X_;
1387 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
1390 if (newstate.
R != Teuchos::null) {
1391 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) < blockSize_,
1392 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
1393 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
1394 std::invalid_argument,
"Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
1395 MVT::SetBlock(*newstate.
R,bsind,*R_);
1398#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1399 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1402 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
1403 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1405 for (
int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1406 MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
1410 Rnorms_current_ =
false;
1411 R2norms_current_ =
false;
1416 AX_ = Teuchos::null;
1420 initialized_ =
true;
1422 if (om_->isVerbosity(
Debug ) ) {
1430 om_->print(
Debug, accuracyCheck(chk,
"after initialize()") );
1434 template <
class ScalarType,
class MV,
class OP>
1446 template <
class ScalarType,
class MV,
class OP>
1447 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1449 if (Rnorms_current_ ==
false) {
1451 orthman_->norm(*R_,Rnorms_);
1452 Rnorms_current_ =
true;
1460 template <
class ScalarType,
class MV,
class OP>
1461 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1463 if (R2norms_current_ ==
false) {
1465 MVT::MvNorm(*R_,R2norms_);
1466 R2norms_current_ =
true;
1498 template <
class ScalarType,
class MV,
class OP>
1501 using std::setprecision;
1502 using std::scientific;
1505 std::stringstream os;
1508 os <<
" Debugging checks: " << where << endl;
1511 if (chk.checkX && initialized_) {
1512 tmp = orthman_->orthonormError(*X_);
1513 os <<
" >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl;
1514 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1515 tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
1516 os <<
" >> Error in X^H B Q[" << i <<
"] == 0 : " << scientific << setprecision(10) << tmp << endl;
1519 if (chk.checkAX && !isSkinny_ && initialized_) {
1520 tmp = Utils::errorEquality(*X_, *AX_, AOp_);
1521 os <<
" >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl;
1523 if (chk.checkBX && hasBOp_ && initialized_) {
1524 Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
1525 OPT::Apply(*BOp_,*this->X_,*BX);
1526 tmp = Utils::errorEquality(*BX, *BX_);
1527 os <<
" >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl;
1531 if (chk.checkEta && initialized_) {
1532 tmp = orthman_->orthogError(*X_,*eta_);
1533 os <<
" >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl;
1534 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1535 tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
1536 os <<
" >> Error in Eta^H B Q[" << i <<
"]==0 : " << scientific << setprecision(10) << tmp << endl;
1539 if (chk.checkAEta && !isSkinny_ && initialized_) {
1540 tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
1541 os <<
" >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
1543 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
1544 tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
1545 os <<
" >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl;
1549 if (chk.checkR && initialized_) {
1550 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1551 MVT::MvTransMv(ONE,*X_,*R_,xTx);
1552 tmp = xTx.normFrobenius();
1553 os <<
" >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl;
1558 if (chk.checkBR && hasBOp_ && initialized_) {
1559 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1560 MVT::MvTransMv(ONE,*BX_,*R_,xTx);
1561 tmp = xTx.normFrobenius();
1562 os <<
" >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl;
1566 if (chk.checkZ && initialized_) {
1567 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1568 MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
1569 tmp = xTx.normFrobenius();
1570 os <<
" >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl;
1575 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1576 tmp = orthman_->orthonormError(*auxVecs_[i]);
1577 os <<
" >> Error in Q[" << i <<
"]^H B Q[" << i <<
"]==I: " << scientific << setprecision(10) << tmp << endl;
1578 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1579 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1580 os <<
" >> Error in Q[" << i <<
"]^H B Q[" << j <<
"]==0: " << scientific << setprecision(10) << tmp << endl;
1591 template <
class ScalarType,
class MV,
class OP>
1595 using std::setprecision;
1596 using std::scientific;
1601 os <<
"================================================================================" << endl;
1603 os <<
" RTRBase Solver Status" << endl;
1605 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
1606 os <<
"The number of iterations performed is " << iter_ << endl;
1607 os <<
"The current block size is " << blockSize_ << endl;
1608 os <<
"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1609 os <<
"The number of operations A*x is " << counterAOp_ << endl;
1610 os <<
"The number of operations B*x is " << counterBOp_ << endl;
1611 os <<
"The number of operations Prec*x is " << counterPrec_ << endl;
1612 os <<
"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
1613 os <<
"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
1617 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
1618 os << setw(20) <<
"Eigenvalue"
1619 << setw(20) <<
"Residual(B)"
1620 << setw(20) <<
"Residual(2)"
1622 os <<
"--------------------------------------------------------------------------------"<<endl;
1623 for (
int i=0; i<blockSize_; i++) {
1624 os << scientific << setprecision(10) << setw(20) << theta_[i];
1625 if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
1626 else os << scientific << setprecision(10) << setw(20) <<
"not current";
1627 if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
1628 else os << scientific << setprecision(10) << setw(20) <<
"not current";
1632 os <<
"================================================================================" << endl;
1639 template <
class ScalarType,
class MV,
class OP>
1640 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1643 std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
1645 MagnitudeType ret = 0;
1646 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1655 template <
class ScalarType,
class MV,
class OP>
1656 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1657 RTRBase<ScalarType,MV,OP>::ginner(
const MV &xi,
const MV &zeta)
const
1659 std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
1660 MVT::MvDot(xi,zeta,d);
1661 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
1667 template <
class ScalarType,
class MV,
class OP>
1668 void RTRBase<ScalarType,MV,OP>::ginnersep(
1670 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d)
const
1673 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1681 template <
class ScalarType,
class MV,
class OP>
1682 void RTRBase<ScalarType,MV,OP>::ginnersep(
1683 const MV &xi,
const MV &zeta,
1684 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d)
const
1686 std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
1687 MVT::MvDot(xi,zeta,dC);
1688 vecMTiter iR=d.begin();
1689 vecSTiter iS=dC.begin();
1690 for (; iR != d.end(); iR++, iS++) {
1691 (*iR) = SCT::real(*iS);
1695 template <
class ScalarType,
class MV,
class OP>
1700 template <
class ScalarType,
class MV,
class OP>
1705 template <
class ScalarType,
class MV,
class OP>
1710 template <
class ScalarType,
class MV,
class OP>
1715 template <
class ScalarType,
class MV,
class OP>
1718 if (!initialized_)
return 0;
1722 template <
class ScalarType,
class MV,
class OP>
1723 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1726 std::vector<MagnitudeType> ret = ritz2norms_;
1727 ret.resize(nevLocal_);
1731 template <
class ScalarType,
class MV,
class OP>
1732 std::vector<Value<ScalarType> >
1735 std::vector<Value<ScalarType> > ret(nevLocal_);
1736 for (
int i=0; i<nevLocal_; i++) {
1737 ret[i].realpart = theta_[i];
1738 ret[i].imagpart = ZERO;
1743 template <
class ScalarType,
class MV,
class OP>
1744 Teuchos::RCP<const MV>
1750 template <
class ScalarType,
class MV,
class OP>
1756 template <
class ScalarType,
class MV,
class OP>
1762 template <
class ScalarType,
class MV,
class OP>
1772 state.
BX = Teuchos::null;
1776 state.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta_));
1780 template <
class ScalarType,
class MV,
class OP>
1783 return initialized_;
1786 template <
class ScalarType,
class MV,
class OP>
1789 std::vector<int> ret(nevLocal_,0);
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Class which provides internal utilities for the Anasazi solvers.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
This class defines the interface required by an eigensolver and status test class to compute solution...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Output managers remove the need for the eigensolver to know any information about the required output...
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers....
virtual void iterate()=0
This method performs RTR iterations until the status test indicates the need to stop or an error occu...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For RTR, this always returns getBlockSiz...
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void resetNumIters()
Reset the iteration count.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
virtual void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
int getNumIters() const
Get the current iteration count.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
virtual ~RTRBase()
RTRBase destructor
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the Ritz residuals.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
RTRBase(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< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms, const std::string &solverLabel, bool skinnySolver)
RTRBase constructor with eigenproblem, solver utilities, and parameter list of solver options.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
RTRState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
RTRInitFailure is thrown when the RTR solver is unable to generate an initial iterate in the RTRBase:...
RTROrthoFailure is thrown when an orthogonalization attempt fails.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
Anasazi's templated, static class providing utilities for the solvers.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Structure to contain pointers to RTR state variables.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MV > BX
The image of the current eigenvectors under B, or Teuchos::null if B was not specified.
Teuchos::ScalarTraits< ScalarType >::magnitudeType rho
The current rho value. This is only valid if the debugging level of verbosity is enabled.
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > AX
The image of the current eigenvectors under A, or Teuchos::null is we implement a skinny solver.