213 RTRBase(
const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
214 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
215 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
216 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
217 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
218 Teuchos::ParameterList ¶ms,
219 const std::string &solverLabel,
bool skinnySolver
277 void initialize(RTRState<ScalarType,MV>& newstate);
305 RTRState<ScalarType,MV>
getState()
const;
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();
386 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
389 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const;
392 const Eigenproblem<ScalarType,MV,OP>&
getProblem()
const;
430 void setAuxVecs(
const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
433 Teuchos::Array<Teuchos::RCP<const MV> >
getAuxVecs()
const;
449 typedef SolverUtils<ScalarType,MV,OP> Utils;
450 typedef MultiVecTraits<ScalarType,MV> MVT;
451 typedef OperatorTraits<ScalarType,MV,OP> OPT;
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_;
609 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
610 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
611 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
612 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
613 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
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.");
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);
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()") );