29 #ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP 30 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP 52 #include "Teuchos_BLAS.hpp" 53 #include "Teuchos_LAPACK.hpp" 54 #include "Teuchos_TimeMonitor.hpp" 56 # include <Teuchos_FancyOStream.hpp> 112 template<
class ScalarType,
class MV,
class OP>
118 typedef Teuchos::ScalarTraits<ScalarType> SCT;
119 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
120 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
150 Teuchos::ParameterList &pl );
176 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
177 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
229 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
231 std::string whch_, ortho_;
233 MagnitudeType convtol_, locktol_;
236 bool relconvtol_, rellocktol_;
237 int blockSize_, numBlocks_, numIters_;
241 int numRestartBlocks_;
242 enum ResType convNorm_, lockNorm_;
244 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
246 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
247 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
248 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
250 Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
255 template<
class ScalarType,
class MV,
class OP>
258 Teuchos::ParameterList &pl ) :
262 convtol_(MT::prec()),
272 inSituRestart_(false),
274 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
275 , _timerSolve(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr::solve()")),
276 _timerRestarting(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr restarting")),
277 _timerLocking(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockDavidsonSolMgr locking"))
280 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
281 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
282 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
283 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
288 whch_ = pl.get(
"Which",whch_);
289 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR",std::invalid_argument,
"Invalid sorting string.");
292 ortho_ = pl.get(
"Orthogonalization",ortho_);
293 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB") {
298 convtol_ = pl.get(
"Convergence Tolerance",convtol_);
299 relconvtol_ = pl.get(
"Relative Convergence Tolerance",relconvtol_);
300 strtmp = pl.get(
"Convergence Norm",std::string(
"2"));
302 convNorm_ = RES_2NORM;
304 else if (strtmp ==
"M") {
305 convNorm_ = RES_ORTH;
308 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
309 "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
313 useLocking_ = pl.get(
"Use Locking",useLocking_);
314 rellocktol_ = pl.get(
"Relative Locking Tolerance",rellocktol_);
316 locktol_ = convtol_/10;
317 locktol_ = pl.get(
"Locking Tolerance",locktol_);
318 strtmp = pl.get(
"Locking Norm",std::string(
"2"));
320 lockNorm_ = RES_2NORM;
322 else if (strtmp ==
"M") {
323 lockNorm_ = RES_ORTH;
326 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
327 "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
331 maxRestarts_ = pl.get(
"Maximum Restarts",maxRestarts_);
334 blockSize_ = pl.get(
"Block Size",problem_->getNEV());
335 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
336 "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
337 numBlocks_ = pl.get(
"Num Blocks",2);
338 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
339 "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
343 maxLocked_ = pl.get(
"Max Locked",problem_->getNEV());
348 if (maxLocked_ == 0) {
351 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
352 "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
353 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
354 std::invalid_argument,
355 "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
356 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ + maxLocked_ >
MVT::GetGlobalLength(*problem_->getInitVec()),
357 std::invalid_argument,
358 "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
361 lockQuorum_ = pl.get(
"Locking Quorum",lockQuorum_);
362 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
363 std::invalid_argument,
364 "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
368 numRestartBlocks_ = pl.get(
"Num Restart Blocks",numRestartBlocks_);
369 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
370 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
371 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
372 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
375 if (pl.isParameter(
"In Situ Restarting")) {
376 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
377 inSituRestart_ = pl.get(
"In Situ Restarting",inSituRestart_);
379 inSituRestart_ = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
384 std::string fntemplate =
"";
385 bool allProcs =
false;
386 if (pl.isParameter(
"Output on all processors")) {
387 if (Teuchos::isParameterType<bool>(pl,
"Output on all processors")) {
388 allProcs = pl.get(
"Output on all processors",allProcs);
390 allProcs = ( Teuchos::getParameter<int>(pl,
"Output on all processors") != 0 );
393 fntemplate = pl.get(
"Output filename template",fntemplate);
398 MPI_Initialized(&mpiStarted);
399 if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
404 if (fntemplate !=
"") {
405 std::ostringstream MyPIDstr;
409 while ( (pos = fntemplate.find(
"%d",start)) != -1 ) {
410 fntemplate.replace(pos,2,MyPIDstr.str());
414 Teuchos::RCP<ostream> osp;
415 if (fntemplate !=
"") {
416 osp = Teuchos::rcp(
new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
418 osp = Teuchos::rcpFromRef(std::cout);
419 std::cout <<
"Anasazi::BlockDavidsonSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
423 osp = Teuchos::rcpFromRef(std::cout);
427 if (pl.isParameter(
"Verbosity")) {
428 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
429 verbosity = pl.get(
"Verbosity", verbosity);
431 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
446 template<
class ScalarType,
class MV,
class OP>
452 const int nev = problem_->getNEV();
455 Teuchos::RCP<Teuchos::FancyOStream>
456 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
457 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
458 *out <<
"Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
469 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
470 if (globalTest_ == Teuchos::null) {
474 convtest = globalTest_;
476 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
479 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
481 if (lockingTest_ == Teuchos::null) {
485 locktest = lockingTest_;
489 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
490 alltests.push_back(ordertest);
491 if (locktest != Teuchos::null) alltests.push_back(locktest);
492 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
494 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
497 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
498 if ( printer_->isVerbosity(
Debug) ) {
507 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
508 if (ortho_==
"SVQB") {
510 }
else if (ortho_==
"DGKS") {
513 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS",std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
518 Teuchos::ParameterList plist;
519 plist.set(
"Block Size",blockSize_);
520 plist.set(
"Num Blocks",numBlocks_);
524 Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
527 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
528 if (probauxvecs != Teuchos::null) {
529 bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
536 int curNumLocked = 0;
537 Teuchos::RCP<MV> lockvecs;
543 if (maxLocked_ > 0) {
544 lockvecs =
MVT::Clone(*problem_->getInitVec(),maxLocked_);
546 std::vector<MagnitudeType> lockvals;
594 Teuchos::RCP<MV> workMV;
595 if (inSituRestart_ ==
false) {
597 if (useLocking_==
true || maxRestarts_ > 0) {
598 workMV =
MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
602 workMV = Teuchos::null;
611 if (maxRestarts_ > 0) {
612 workMV =
MVT::Clone(*problem_->getInitVec(),1);
615 workMV = Teuchos::null;
620 const ScalarType ONE = SCT::one();
621 const ScalarType ZERO = SCT::zero();
622 Teuchos::LAPACK<int,ScalarType> lapack;
623 Teuchos::BLAS<int,ScalarType> blas;
628 problem_->setSolution(sol);
634 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 635 Teuchos::TimeMonitor slvtimer(*_timerSolve);
641 bd_solver->iterate();
648 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
649 throw AnasaziError(
"Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
656 else if (ordertest->getStatus() ==
Passed ) {
667 else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
669 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 670 Teuchos::TimeMonitor restimer(*_timerRestarting);
673 if ( numRestarts >= maxRestarts_ ) {
678 printer_->stream(
IterationDetails) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
681 int curdim = state.
curDim;
682 int newdim = numRestartBlocks_*blockSize_;
684 # ifdef TEUCHOS_DEBUG 686 std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
687 *out <<
"Ritz values from solver:\n";
688 for (
unsigned int i=0; i<ritzvalues.size(); i++) {
689 *out << ritzvalues[i].realpart <<
" ";
697 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
698 std::vector<MagnitudeType> theta(curdim);
700 # ifdef TEUCHOS_DEBUG 702 std::stringstream os;
703 os <<
"KK before HEEV...\n" 704 << *state.
KK <<
"\n";
708 int info = msutils::directSolver(curdim,*state.
KK,Teuchos::null,S,theta,rank,10);
709 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
710 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
711 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
712 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
714 # ifdef TEUCHOS_DEBUG 716 std::stringstream os;
717 *out <<
"Ritz values from heev(KK):\n";
718 for (
unsigned int i=0; i<theta.size(); i++) *out << theta[i] <<
" ";
719 os <<
"\nRitz vectors from heev(KK):\n" 728 std::vector<int> order(curdim);
729 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
732 msutils::permuteVectors(order,S);
734 # ifdef TEUCHOS_DEBUG 736 std::stringstream os;
737 *out <<
"Ritz values from heev(KK) after sorting:\n";
738 std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out,
" "));
739 os <<
"\nRitz vectors from heev(KK) after sorting:\n" 746 Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
747 # ifdef TEUCHOS_DEBUG 749 std::stringstream os;
750 os <<
"Significant primitive Ritz vectors:\n" 757 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
759 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
760 KKold(Teuchos::View,*state.
KK,curdim,curdim);
763 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
764 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
765 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
767 teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
768 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
769 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
771 for (
int j=0; j<newdim-1; ++j) {
772 for (
int i=j+1; i<newdim; ++i) {
773 newKK(i,j) = SCT::conjugate(newKK(j,i));
777 # ifdef TEUCHOS_DEBUG 779 std::stringstream os;
789 rstate.
KK = Teuchos::rcpFromRef(newKK);
796 rstate.
KX = state.
KX;
797 rstate.
MX = state.
MX;
799 rstate.
T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
801 if (inSituRestart_ ==
true) {
802 # ifdef TEUCHOS_DEBUG 803 *out <<
"Beginning in-situ restart...\n";
807 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.
V);
811 std::vector<ScalarType> tau(newdim), work(newdim);
813 lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
814 TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
815 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
816 if (printer_->isVerbosity(
Debug)) {
817 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
818 for (
int j=0; j<newdim; j++) {
819 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
820 for (
int i=j+1; i<newdim; i++) {
824 printer_->stream(
Debug) <<
"||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
831 std::vector<int> curind(curdim);
832 for (
int i=0; i<curdim; i++) curind[i] = i;
834 msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
840 rstate.
V = solverbasis;
843 # ifdef TEUCHOS_DEBUG 844 *out <<
"Beginning ex-situ restart...\n";
847 std::vector<int> curind(curdim), newind(newdim);
848 for (
int i=0; i<curdim; i++) curind[i] = i;
849 for (
int i=0; i<newdim; i++) newind[i] = i;
861 bd_solver->initialize(rstate);
868 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
870 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 871 Teuchos::TimeMonitor lcktimer(*_timerLocking);
877 const int curdim = state.
curDim;
881 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
882 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
883 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
884 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
886 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
887 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
890 std::vector<int> tmp_vector_int;
891 if (curNumLocked + locktest->howMany() > maxLocked_) {
893 tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
896 tmp_vector_int = locktest->whichVecs();
898 const std::vector<int> lockind(tmp_vector_int);
899 const int numNewLocked = lockind.size();
903 const int numUnlocked = curdim-numNewLocked;
904 tmp_vector_int.resize(curdim);
905 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
906 const std::vector<int> curind(tmp_vector_int);
907 tmp_vector_int.resize(numUnlocked);
908 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
909 const std::vector<int> unlockind(tmp_vector_int);
910 tmp_vector_int.clear();
914 if (printer_->isVerbosity(
Debug)) {
915 printer_->print(
Debug,
"Locking vectors: ");
916 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
917 printer_->print(
Debug,
"\n");
918 printer_->print(
Debug,
"Not locking vectors: ");
919 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
920 printer_->print(
Debug,
"\n");
931 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
932 std::vector<MagnitudeType> theta(curdim);
935 int info = msutils::directSolver(curdim,*state.
KK,Teuchos::null,S,theta,rank,10);
936 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
937 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
938 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
939 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
942 std::vector<int> order(curdim);
943 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
946 msutils::permuteVectors(order,S);
952 Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
953 for (
int i=0; i<numUnlocked; i++) {
954 blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
967 Teuchos::RCP<MV> defV, augV;
968 if (inSituRestart_ ==
true) {
971 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.
V);
975 Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
976 std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
978 lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
979 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
980 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
981 if (printer_->isVerbosity(
Debug)) {
982 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
983 for (
int j=0; j<numUnlocked; j++) {
984 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
985 for (
int i=j+1; i<numUnlocked; i++) {
989 printer_->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
997 msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
999 std::vector<int> defind(numUnlocked), augind(numNewLocked);
1000 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
1001 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
1007 std::vector<int> defind(numUnlocked), augind(numNewLocked);
1008 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
1009 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
1028 Teuchos::RCP<const MV> curlocked, newLocked;
1029 Teuchos::RCP<MV> augTmp;
1032 if (curNumLocked > 0) {
1033 std::vector<int> curlockind(curNumLocked);
1034 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1038 curlocked = Teuchos::null;
1041 std::vector<int> augtmpind(numNewLocked);
1042 for (
int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1045 newLocked =
MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1056 Teuchos::Array<Teuchos::RCP<const MV> > against;
1057 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1058 if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1059 if (curlocked != Teuchos::null) against.push_back(curlocked);
1060 against.push_back(newLocked);
1061 against.push_back(defV);
1062 if (problem_->getM() != Teuchos::null) {
1065 ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1076 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1078 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1079 KKold(Teuchos::View,*state.
KK,curdim,curdim),
1080 KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1083 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1084 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1085 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1087 teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1088 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1089 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1094 OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1095 Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1096 KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1102 defV = Teuchos::null;
1103 augV = Teuchos::null;
1106 for (
int j=0; j<curdim; ++j) {
1107 for (
int i=j+1; i<curdim; ++i) {
1108 newKK(i,j) = SCT::conjugate(newKK(j,i));
1114 augTmp = Teuchos::null;
1116 std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1117 for (
int i=0; i<numNewLocked; i++) {
1118 lockvals.push_back(allvals[lockind[i]].realpart);
1121 std::vector<int> indlock(numNewLocked);
1122 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1124 newLocked = Teuchos::null;
1126 curNumLocked += numNewLocked;
1127 std::vector<int> curlockind(curNumLocked);
1128 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1135 ordertest->setAuxVals(lockvals);
1137 Teuchos::Array< Teuchos::RCP<const MV> > aux;
1138 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1139 aux.push_back(curlocked);
1140 bd_solver->setAuxVecs(aux);
1142 if (curNumLocked == maxLocked_) {
1144 combotest->removeTest(locktest);
1152 if (inSituRestart_) {
1160 rstate.
KK = Teuchos::rcpFromRef(newKK);
1163 bd_solver->initialize(rstate);
1172 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1177 <<
"Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1178 << err.what() << std::endl
1179 <<
"Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1185 workMV = Teuchos::null;
1187 sol.
numVecs = ordertest->howMany();
1192 std::vector<MagnitudeType> vals(sol.
numVecs);
1195 std::vector<int> which = ordertest->whichVecs();
1199 std::vector<int> inlocked(0), insolver(0);
1200 for (
unsigned int i=0; i<which.size(); i++) {
1201 if (which[i] >= 0) {
1202 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1203 insolver.push_back(which[i]);
1207 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1208 inlocked.push_back(which[i] + curNumLocked);
1212 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.
numVecs,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1215 if (insolver.size() > 0) {
1217 int lclnum = insolver.size();
1218 std::vector<int> tosol(lclnum);
1219 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1220 Teuchos::RCP<const MV> v =
MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1223 std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1224 for (
unsigned int i=0; i<insolver.size(); i++) {
1225 vals[i] = fromsolver[insolver[i]].realpart;
1230 if (inlocked.size() > 0) {
1231 int solnum = insolver.size();
1233 int lclnum = inlocked.size();
1234 std::vector<int> tosol(lclnum);
1235 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1239 for (
unsigned int i=0; i<inlocked.size(); i++) {
1240 vals[i+solnum] = lockvals[inlocked[i]];
1246 std::vector<int> order(sol.
numVecs);
1247 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.
numVecs);
1249 for (
int i=0; i<sol.
numVecs; i++) {
1250 sol.
Evals[i].realpart = vals[i];
1251 sol.
Evals[i].imagpart = MT::zero();
1263 bd_solver->currentStatus(printer_->stream(
FinalSummary));
1266 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1268 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1272 problem_->setSolution(sol);
1273 printer_->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1276 numIters_ = bd_solver->getNumIters();
1285 template <
class ScalarType,
class MV,
class OP>
1290 globalTest_ = global;
1293 template <
class ScalarType,
class MV,
class OP>
1294 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1300 template <
class ScalarType,
class MV,
class OP>
1308 template <
class ScalarType,
class MV,
class OP>
1309 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1315 template <
class ScalarType,
class MV,
class OP>
1320 lockingTest_ = locking;
1323 template <
class ScalarType,
class MV,
class OP>
1324 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1327 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. ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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...
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...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
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.
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
ReturnType
Enumerated type used to pass back information from a solver manager.
Teuchos::RCP< const MV > V
The basis for the Krylov space.
A status test for testing the norm of the eigenvectors residuals.
virtual ~BlockDavidsonSolMgr()
Destructor.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
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 > > & getLockingStatusTest() const
Get the status test defining locking.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
Struct for storing an eigenproblem solution.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
int curDim
The current dimension of the solver.
Structure to contain pointers to BlockDavidson state variables.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
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.
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 ...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.