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" 120 template<
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_;
266 template<
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");
425 template<
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_;
663 # ifdef TEUCHOS_DEBUG 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);
679 # ifdef TEUCHOS_DEBUG 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.");
693 # ifdef TEUCHOS_DEBUG 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);
713 # ifdef TEUCHOS_DEBUG 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);
726 # ifdef TEUCHOS_DEBUG 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));
756 # ifdef TEUCHOS_DEBUG 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) {
781 # ifdef TEUCHOS_DEBUG 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;
822 # ifdef TEUCHOS_DEBUG 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();
1264 template <
class ScalarType,
class MV,
class OP>
1269 globalTest_ = global;
1272 template <
class ScalarType,
class MV,
class OP>
1273 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1279 template <
class ScalarType,
class MV,
class OP>
1287 template <
class ScalarType,
class MV,
class OP>
1288 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1294 template <
class ScalarType,
class MV,
class OP>
1299 lockingTest_ = locking;
1302 template <
class ScalarType,
class MV,
class OP>
1303 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1306 return lockingTest_;
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Pure virtual base class which describes the basic interface for a solver manager. ...
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
A special StatusTest for printing other status tests.
BlockDavidsonSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockDavidsonSolMgr.
This class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
Status test for forming logical combinations of other status tests.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
An exception class parent to all Anasazi exceptions.
Output managers remove the need for the eigensolver to know any information about the required output...
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
Implementation of the block Davidson method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi's templated, static class providing utilities for the solvers.
int numVecs
The number of computed eigenpairs.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Abstract class definition for Anasazi Output Managers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
ReturnType
Enumerated type used to pass back information from a solver manager.
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::RCP< const MV > V
The basis for the Krylov space.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A status test for testing the norm of the eigenvectors residuals.
virtual ~BlockDavidsonSolMgr()
Destructor.
Traits class which defines basic operations on multivectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
int curDim
The current dimension of the solver.
Structure to contain pointers to BlockDavidson state variables.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
int getNumIters() const
Get the iteration count for the most recent call to solve().
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Status test for forming logical combinations of other status tests.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
Abstract class definition for Anasazi output stream.
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson eigensolver...
Common interface of stopping criteria for Anasazi's solvers.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
Class which provides internal utilities for the Anasazi solvers.