42#ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
43#define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
66#include "Teuchos_BLAS.hpp"
67#include "Teuchos_LAPACK.hpp"
68#include "Teuchos_TimeMonitor.hpp"
69#include "Teuchos_FancyOStream.hpp"
120template<
class ScalarType,
class MV,
class OP>
126 typedef Teuchos::ScalarTraits<ScalarType> SCT;
127 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
128 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
161 Teuchos::ParameterList &pl );
187 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
188 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
240 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
242 std::string whch_, ortho_;
244 MagnitudeType convtol_, locktol_;
247 bool relconvtol_, rellocktol_;
248 int blockSize_, numBlocks_, numIters_;
252 int numRestartBlocks_;
253 enum ResType convNorm_, lockNorm_;
255 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
257 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
258 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
259 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
261 Teuchos::RCP<OutputManager<ScalarType> > printer_;
266template<
class ScalarType,
class MV,
class OP>
269 Teuchos::ParameterList &pl ) :
273 convtol_(MT::prec()),
283 inSituRestart_(false),
285#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
286 , _timerSolve(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr::solve()")),
287 _timerRestarting(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr restarting")),
288 _timerLocking(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr locking"))
291 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
292 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
293 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
294 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
299 whch_ = pl.get(
"Which",whch_);
300 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR",std::invalid_argument,
"Invalid sorting string.");
303 ortho_ = pl.get(
"Orthogonalization",ortho_);
304 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB") {
309 convtol_ = pl.get(
"Convergence Tolerance",convtol_);
310 relconvtol_ = pl.get(
"Relative Convergence Tolerance",relconvtol_);
311 strtmp = pl.get(
"Convergence Norm",std::string(
"2"));
313 convNorm_ = RES_2NORM;
315 else if (strtmp ==
"M") {
316 convNorm_ = RES_ORTH;
319 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
320 "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
324 useLocking_ = pl.get(
"Use Locking",useLocking_);
325 rellocktol_ = pl.get(
"Relative Locking Tolerance",rellocktol_);
327 locktol_ = convtol_/10;
328 locktol_ = pl.get(
"Locking Tolerance",locktol_);
329 strtmp = pl.get(
"Locking Norm",std::string(
"2"));
331 lockNorm_ = RES_2NORM;
333 else if (strtmp ==
"M") {
334 lockNorm_ = RES_ORTH;
337 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
338 "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
342 maxRestarts_ = pl.get(
"Maximum Restarts",maxRestarts_);
345 blockSize_ = pl.get(
"Block Size",problem_->getNEV());
346 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
347 "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
348 numBlocks_ = pl.get(
"Num Blocks",2);
349 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
350 "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
354 maxLocked_ = pl.get(
"Max Locked",problem_->getNEV());
359 if (maxLocked_ == 0) {
362 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
363 "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
364 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
365 std::invalid_argument,
366 "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
367 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_)*blockSize_ + maxLocked_ >
MVT::GetGlobalLength(*problem_->getInitVec()),
368 std::invalid_argument,
369 "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
372 lockQuorum_ = pl.get(
"Locking Quorum",lockQuorum_);
373 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
374 std::invalid_argument,
375 "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
379 numRestartBlocks_ = pl.get(
"Num Restart Blocks",numRestartBlocks_);
380 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
381 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
382 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
383 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
386 if (pl.isParameter(
"In Situ Restarting")) {
387 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
388 inSituRestart_ = pl.get(
"In Situ Restarting",inSituRestart_);
390 inSituRestart_ = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
397 int osProc = pl.get(
"Output Processor", 0);
400 Teuchos::RCP<Teuchos::FancyOStream> osp;
403 if (pl.isParameter(
"Output Stream")) {
404 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
412 if (pl.isParameter(
"Verbosity")) {
413 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
414 verbosity = pl.get(
"Verbosity", verbosity);
416 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
425template<
class ScalarType,
class MV,
class OP>
431 const int nev = problem_->getNEV();
434 Teuchos::RCP<Teuchos::FancyOStream>
435 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
436 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
437 *out <<
"Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
448 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
449 if (globalTest_ == Teuchos::null) {
453 convtest = globalTest_;
455 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
458 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
460 if (lockingTest_ == Teuchos::null) {
464 locktest = lockingTest_;
468 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
469 alltests.push_back(ordertest);
470 if (locktest != Teuchos::null) alltests.push_back(locktest);
471 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
473 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
476 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
477 if ( printer_->isVerbosity(
Debug) ) {
486 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
487 if (ortho_==
"SVQB") {
489 }
else if (ortho_==
"DGKS") {
492 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS",std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
497 Teuchos::ParameterList plist;
498 plist.set(
"Block Size",blockSize_);
499 plist.set(
"Num Blocks",numBlocks_);
503 Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
506 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
507 if (probauxvecs != Teuchos::null) {
508 bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
515 int curNumLocked = 0;
516 Teuchos::RCP<MV> lockvecs;
522 if (maxLocked_ > 0) {
523 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
525 std::vector<MagnitudeType> lockvals;
573 Teuchos::RCP<MV> workMV;
574 if (inSituRestart_ ==
false) {
576 if (useLocking_==
true || maxRestarts_ > 0) {
577 workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
581 workMV = Teuchos::null;
590 if (maxRestarts_ > 0) {
591 workMV = MVT::Clone(*problem_->getInitVec(),1);
594 workMV = Teuchos::null;
599 const ScalarType ONE = SCT::one();
600 const ScalarType ZERO = SCT::zero();
601 Teuchos::LAPACK<int,ScalarType> lapack;
602 Teuchos::BLAS<int,ScalarType> blas;
607 problem_->setSolution(sol);
613#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614 Teuchos::TimeMonitor slvtimer(*_timerSolve);
620 bd_solver->iterate();
627 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
628 throw AnasaziError(
"Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
635 else if (ordertest->getStatus() ==
Passed ) {
646 else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
648#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
649 Teuchos::TimeMonitor restimer(*_timerRestarting);
652 if ( numRestarts >= maxRestarts_ ) {
657 printer_->stream(
IterationDetails) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
660 int curdim = state.
curDim;
661 int newdim = numRestartBlocks_*blockSize_;
665 std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
666 *out <<
"Ritz values from solver:\n";
667 for (
unsigned int i=0; i<ritzvalues.size(); i++) {
668 *out << ritzvalues[i].realpart <<
" ";
676 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
677 std::vector<MagnitudeType> theta(curdim);
681 std::stringstream os;
682 os <<
"KK before HEEV...\n"
683 << printMat(*state.
KK) <<
"\n";
687 int info = msutils::directSolver(curdim,*state.
KK,Teuchos::null,S,theta,rank,10);
688 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
689 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
690 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
691 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
695 std::stringstream os;
696 *out <<
"Ritz values from heev(KK):\n";
697 for (
unsigned int i=0; i<theta.size(); i++) *out << theta[i] <<
" ";
698 os <<
"\nRitz vectors from heev(KK):\n"
699 << printMat(S) <<
"\n";
707 std::vector<int> order(curdim);
708 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
711 msutils::permuteVectors(order,S);
715 std::stringstream os;
716 *out <<
"Ritz values from heev(KK) after sorting:\n";
717 std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out,
" "));
718 os <<
"\nRitz vectors from heev(KK) after sorting:\n"
719 << printMat(S) <<
"\n";
725 Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
728 std::stringstream os;
729 os <<
"Significant primitive Ritz vectors:\n"
730 << printMat(Sr) <<
"\n";
736 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
738 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
739 KKold(Teuchos::View,*state.
KK,curdim,curdim);
742 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
743 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
744 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
746 teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
747 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
748 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
750 for (
int j=0; j<newdim-1; ++j) {
751 for (
int i=j+1; i<newdim; ++i) {
752 newKK(i,j) = SCT::conjugate(newKK(j,i));
758 std::stringstream os;
760 << printMat(newKK) <<
"\n";
768 rstate.
KK = Teuchos::rcpFromRef(newKK);
775 rstate.
KX = state.
KX;
776 rstate.
MX = state.
MX;
778 rstate.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
780 if (inSituRestart_ ==
true) {
782 *out <<
"Beginning in-situ restart...\n";
786 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.
V);
790 std::vector<ScalarType> tau(newdim), work(newdim);
792 lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
793 TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
794 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
795 if (printer_->isVerbosity(
Debug)) {
796 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
797 for (
int j=0; j<newdim; j++) {
798 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
799 for (
int i=j+1; i<newdim; i++) {
803 printer_->stream(
Debug) <<
"||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
810 std::vector<int> curind(curdim);
811 for (
int i=0; i<curdim; i++) curind[i] = i;
812 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
813 msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
819 rstate.
V = solverbasis;
823 *out <<
"Beginning ex-situ restart...\n";
826 std::vector<int> curind(curdim), newind(newdim);
827 for (
int i=0; i<curdim; i++) curind[i] = i;
828 for (
int i=0; i<newdim; i++) newind[i] = i;
829 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.
V,curind);
830 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*workMV ,newind);
832 MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
840 bd_solver->initialize(rstate);
847 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
849#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
850 Teuchos::TimeMonitor lcktimer(*_timerLocking);
856 const int curdim = state.
curDim;
860 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
861 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
862 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (
int)locktest->whichVecs().size(),std::logic_error,
863 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
865 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
866 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
869 std::vector<int> tmp_vector_int;
870 if (curNumLocked + locktest->howMany() > maxLocked_) {
872 tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
875 tmp_vector_int = locktest->whichVecs();
877 const std::vector<int> lockind(tmp_vector_int);
878 const int numNewLocked = lockind.size();
882 const int numUnlocked = curdim-numNewLocked;
883 tmp_vector_int.resize(curdim);
884 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
885 const std::vector<int> curind(tmp_vector_int);
886 tmp_vector_int.resize(numUnlocked);
887 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
888 const std::vector<int> unlockind(tmp_vector_int);
889 tmp_vector_int.clear();
893 if (printer_->isVerbosity(
Debug)) {
894 printer_->print(
Debug,
"Locking vectors: ");
895 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
896 printer_->print(
Debug,
"\n");
897 printer_->print(
Debug,
"Not locking vectors: ");
898 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
899 printer_->print(
Debug,
"\n");
910 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
911 std::vector<MagnitudeType> theta(curdim);
914 int info = msutils::directSolver(curdim,*state.
KK,Teuchos::null,S,theta,rank,10);
915 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
916 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
917 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
918 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
921 std::vector<int> order(curdim);
922 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
925 msutils::permuteVectors(order,S);
931 Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
932 for (
int i=0; i<numUnlocked; i++) {
933 blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
946 Teuchos::RCP<MV> defV, augV;
947 if (inSituRestart_ ==
true) {
950 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.
V);
954 Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
955 std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
957 lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
958 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
959 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
960 if (printer_->isVerbosity(
Debug)) {
961 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
962 for (
int j=0; j<numUnlocked; j++) {
963 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
964 for (
int i=j+1; i<numUnlocked; i++) {
968 printer_->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
975 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
976 msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
978 std::vector<int> defind(numUnlocked), augind(numNewLocked);
979 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
980 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
981 defV = MVT::CloneViewNonConst(*solverbasis,defind);
982 augV = MVT::CloneViewNonConst(*solverbasis,augind);
986 std::vector<int> defind(numUnlocked), augind(numNewLocked);
987 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
988 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
989 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.
V,curind);
990 defV = MVT::CloneViewNonConst(*workMV,defind);
991 augV = MVT::CloneViewNonConst(*workMV,augind);
993 MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
1007 Teuchos::RCP<const MV> curlocked, newLocked;
1008 Teuchos::RCP<MV> augTmp;
1011 if (curNumLocked > 0) {
1012 std::vector<int> curlockind(curNumLocked);
1013 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1014 curlocked = MVT::CloneView(*lockvecs,curlockind);
1017 curlocked = Teuchos::null;
1020 std::vector<int> augtmpind(numNewLocked);
1021 for (
int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1022 augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
1024 newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1030 MVT::MvRandom(*augV);
1035 Teuchos::Array<Teuchos::RCP<const MV> > against;
1036 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1037 if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1038 if (curlocked != Teuchos::null) against.push_back(curlocked);
1039 against.push_back(newLocked);
1040 against.push_back(defV);
1041 if (problem_->getM() != Teuchos::null) {
1042 OPT::Apply(*problem_->getM(),*augV,*augTmp);
1044 ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1055 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1057 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1058 KKold(Teuchos::View,*state.
KK,curdim,curdim),
1059 KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1062 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1063 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1064 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1066 teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1067 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1068 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1073 OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1074 Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1075 KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1076 MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
1077 MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
1081 defV = Teuchos::null;
1082 augV = Teuchos::null;
1085 for (
int j=0; j<curdim; ++j) {
1086 for (
int i=j+1; i<curdim; ++i) {
1087 newKK(i,j) = SCT::conjugate(newKK(j,i));
1093 augTmp = Teuchos::null;
1095 std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1096 for (
int i=0; i<numNewLocked; i++) {
1097 lockvals.push_back(allvals[lockind[i]].realpart);
1100 std::vector<int> indlock(numNewLocked);
1101 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1102 MVT::SetBlock(*newLocked,indlock,*lockvecs);
1103 newLocked = Teuchos::null;
1105 curNumLocked += numNewLocked;
1106 std::vector<int> curlockind(curNumLocked);
1107 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1108 curlocked = MVT::CloneView(*lockvecs,curlockind);
1114 ordertest->setAuxVals(lockvals);
1116 Teuchos::Array< Teuchos::RCP<const MV> > aux;
1117 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1118 aux.push_back(curlocked);
1119 bd_solver->setAuxVecs(aux);
1121 if (curNumLocked == maxLocked_) {
1123 combotest->removeTest(locktest);
1131 if (inSituRestart_) {
1139 rstate.
KK = Teuchos::rcpFromRef(newKK);
1142 bd_solver->initialize(rstate);
1151 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1156 <<
"Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1157 << err.what() << std::endl
1158 <<
"Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1164 workMV = Teuchos::null;
1166 sol.
numVecs = ordertest->howMany();
1168 sol.
Evecs = MVT::Clone(*problem_->getInitVec(),sol.
numVecs);
1171 std::vector<MagnitudeType> vals(sol.
numVecs);
1174 std::vector<int> which = ordertest->whichVecs();
1178 std::vector<int> inlocked(0), insolver(0);
1179 for (
unsigned int i=0; i<which.size(); i++) {
1180 if (which[i] >= 0) {
1181 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1182 insolver.push_back(which[i]);
1186 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1187 inlocked.push_back(which[i] + curNumLocked);
1191 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.
numVecs,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1194 if (insolver.size() > 0) {
1196 int lclnum = insolver.size();
1197 std::vector<int> tosol(lclnum);
1198 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1199 Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1200 MVT::SetBlock(*v,tosol,*sol.
Evecs);
1202 std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1203 for (
unsigned int i=0; i<insolver.size(); i++) {
1204 vals[i] = fromsolver[insolver[i]].realpart;
1209 if (inlocked.size() > 0) {
1210 int solnum = insolver.size();
1212 int lclnum = inlocked.size();
1213 std::vector<int> tosol(lclnum);
1214 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1215 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1216 MVT::SetBlock(*v,tosol,*sol.
Evecs);
1218 for (
unsigned int i=0; i<inlocked.size(); i++) {
1219 vals[i+solnum] = lockvals[inlocked[i]];
1225 std::vector<int> order(sol.
numVecs);
1226 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.
numVecs);
1228 for (
int i=0; i<sol.
numVecs; i++) {
1229 sol.
Evals[i].realpart = vals[i];
1230 sol.
Evals[i].imagpart = MT::zero();
1242 bd_solver->currentStatus(printer_->stream(
FinalSummary));
1245#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1247 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1251 problem_->setSolution(sol);
1252 printer_->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1255 numIters_ = bd_solver->getNumIters();
1264template <
class ScalarType,
class MV,
class OP>
1269 globalTest_ = global;
1272template <
class ScalarType,
class MV,
class OP>
1273const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1279template <
class ScalarType,
class MV,
class OP>
1287template <
class ScalarType,
class MV,
class OP>
1288const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1294template <
class ScalarType,
class MV,
class OP>
1299 lockingTest_ = locking;
1302template <
class ScalarType,
class MV,
class OP>
1303const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1306 return lockingTest_;
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
Implementation of the block Davidson method.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
Class which provides internal utilities for the Anasazi solvers.
Status test for forming logical combinations of other status tests.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson eigensolver.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
BlockDavidsonSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockDavidsonSolMgr.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
virtual ~BlockDavidsonSolMgr()
Destructor.
int getNumIters() const
Get the iteration count for the most recent call to solve().
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
This class defines the interface required by an eigensolver and status test class to compute solution...
Traits class which defines basic operations on multivectors.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Anasazi's templated, static class providing utilities for the solvers.
Status test for forming logical combinations of other status tests.
A special StatusTest for printing other status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
Structure to contain pointers to BlockDavidson state variables.
int curDim
The current dimension of the solver.
Teuchos::RCP< const MV > V
The basis for the Krylov space.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Output managers remove the need for the eigensolver to know any information about the required output...