477 typedef SolverUtils<ScalarType,MV,OP> msutils;
479 const int nev = problem_->getNEV();
485 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp(
new BasicSort<MagnitudeType>(whch_) );
489 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp(
new OutputManager<ScalarType>(verbosity_,osp_) );
495 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
497 maxtest = Teuchos::rcp(
new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
500 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
501 if (globalTest_ == Teuchos::null) {
502 convtest = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
505 convtest = globalTest_;
507 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
508 = Teuchos::rcp(
new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
510 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
512 if (lockingTest_ == Teuchos::null) {
513 locktest = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
516 locktest = lockingTest_;
520 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
521 alltests.push_back(ordertest);
522 if (locktest != Teuchos::null) alltests.push_back(locktest);
523 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
524 if (maxtest != Teuchos::null) alltests.push_back(maxtest);
525 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
528 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
529 if ( printer->isVerbosity(
Debug) ) {
530 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,
Passed+
Failed+
Undefined ) );
533 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,
Passed ) );
538 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
539 if (ortho_==
"SVQB") {
540 ortho = Teuchos::rcp(
new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
541 }
else if (ortho_==
"DGKS") {
542 ortho = Teuchos::rcp(
new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
543 }
else if (ortho_==
"ICGS") {
544 ortho = Teuchos::rcp(
new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
546 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS"&&ortho_!=
"ICGS",std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
551 Teuchos::ParameterList plist;
552 plist.set(
"Block Size",blockSize_);
553 plist.set(
"Full Ortho",fullOrtho_);
557 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
558 = Teuchos::rcp(
new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
560 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
561 if (probauxvecs != Teuchos::null) {
562 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
569 int curNumLocked = 0;
570 Teuchos::RCP<MV> lockvecs;
572 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
574 std::vector<MagnitudeType> lockvals;
583 Teuchos::RCP<MV> workMV;
584 if (fullOrtho_ ==
false && recover_ ==
true) {
585 workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
587 else if (useLocking_) {
588 if (problem_->getM() != Teuchos::null) {
589 workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
592 workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
597 Eigensolution<ScalarType,MV> sol;
599 problem_->setSolution(sol);
602 if (state_ != Teuchos::null) {
603 lobpcg_solver->initialize(*state_);
608#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
609 Teuchos::TimeMonitor slvtimer(*_timerSolve);
615 lobpcg_solver->iterate();
622 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
623 throw AnasaziError(
"Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
630 else if (ordertest->getStatus() ==
Passed || (maxtest != Teuchos::null && maxtest->getStatus() ==
Passed) ) {
641 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
643#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
644 Teuchos::TimeMonitor lcktimer(*_timerLocking);
648 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
649 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
650 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (
int)locktest->whichVecs().size(),std::logic_error,
651 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
652 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
653 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
655 int numnew = locktest->howMany();
656 std::vector<int> indnew = locktest->whichVecs();
659 if (curNumLocked + numnew > maxLocked_) {
660 numnew = maxLocked_ - curNumLocked;
661 indnew.resize(numnew);
666 bool hadP = lobpcg_solver->hasP();
670 printer->print(
Debug,
"Locking vectors: ");
671 for (
unsigned int i=0; i<indnew.size(); i++) {printer->stream(
Debug) <<
" " << indnew[i];}
672 printer->print(
Debug,
"\n");
674 std::vector<MagnitudeType> newvals(numnew);
675 Teuchos::RCP<const MV> newvecs;
679 newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
681 std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
682 for (
int i=0; i<numnew; i++) {
683 newvals[i] = allvals[indnew[i]].realpart;
688 std::vector<int> indlock(numnew);
689 for (
int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
690 MVT::SetBlock(*newvecs,indlock,*lockvecs);
691 newvecs = Teuchos::null;
694 lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
695 curNumLocked += numnew;
698 std::vector<int> indlock(curNumLocked);
699 for (
int i=0; i<curNumLocked; i++) indlock[i] = i;
700 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
701 if (probauxvecs != Teuchos::null) {
702 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
705 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
709 ordertest->setAuxVals(lockvals);
712 LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
713 Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
718 std::vector<int> bsind(blockSize_);
719 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
720 newstateX = MVT::CloneViewNonConst(*workMV,bsind);
721 MVT::SetBlock(*state.X,bsind,*newstateX);
723 if (state.MX != Teuchos::null) {
724 std::vector<int> block3(blockSize_);
725 for (
int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
726 newstateMX = MVT::CloneViewNonConst(*workMV,block3);
727 MVT::SetBlock(*state.MX,bsind,*newstateMX);
732 Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
733 MVT::MvRandom(*newX);
735 if (newstateMX != Teuchos::null) {
736 Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
737 OPT::Apply(*problem_->getM(),*newX,*newMX);
741 Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
742 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
744 ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
749 std::vector<int> block2(blockSize_);
750 for (
int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
751 newstateP = MVT::CloneViewNonConst(*workMV,block2);
752 MVT::SetBlock(*state.P,bsind,*newstateP);
754 if (state.MP != Teuchos::null) {
755 std::vector<int> block4(blockSize_);
756 for (
int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
757 newstateMP = MVT::CloneViewNonConst(*workMV,block4);
758 MVT::SetBlock(*state.MP,bsind,*newstateMP);
763 curauxvecs.push_back(newstateX);
764 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
768 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
772 LOBPCGState<ScalarType,MV> newstate;
773 newstate.X = newstateX;
774 newstate.MX = newstateMX;
775 newstate.P = newstateP;
776 newstate.MP = newstateMP;
777 lobpcg_solver->initialize(newstate);
780 if (curNumLocked == maxLocked_) {
782 combotest->removeTest(locktest);
786 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
795 if (fullOrtho_==
true || recover_==
false) {
799 printer->stream(
Warnings) <<
"Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
800 <<
"Will not try to recover." << std::endl;
803 printer->stream(
Warnings) <<
"Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
804 <<
"Full orthogonalization is off; will try to recover." << std::endl;
809 LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
810 Teuchos::RCP<MV> restart, Krestart, Mrestart;
811 int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
812 bool hasM = problem_->getM() != Teuchos::null;
814 std::vector<int> recind(localsize);
815 for (
int i=0; i<localsize; i++) recind[i] = i;
816 restart = MVT::CloneViewNonConst(*workMV,recind);
819 std::vector<int> recind(localsize);
820 for (
int i=0; i<localsize; i++) recind[i] = localsize+i;
821 Krestart = MVT::CloneViewNonConst(*workMV,recind);
834 std::vector<int> blk1(blockSize_);
835 for (
int i=0; i < blockSize_; i++) blk1[i] = i;
836 MVT::SetBlock(*curstate.X,blk1,*restart);
840 MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
846 std::vector<int> blk2(blockSize_);
847 for (
int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
848 MVT::SetBlock(*curstate.H,blk2,*restart);
852 MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
856 if (localsize == 3*blockSize_) {
857 std::vector<int> blk3(blockSize_);
858 for (
int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
859 MVT::SetBlock(*curstate.P,blk3,*restart);
863 MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
867 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
868 Teuchos::Array<Teuchos::RCP<const MV> > Q;
870 if (curNumLocked > 0) {
871 std::vector<int> indlock(curNumLocked);
872 for (
int i=0; i<curNumLocked; i++) indlock[i] = i;
873 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
874 Q.push_back(curlocked);
876 if (probauxvecs != Teuchos::null) {
877 Q.push_back(probauxvecs);
880 int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
881 if (rank < blockSize_) {
883 printer->stream(
Errors) <<
"Error! Recovered basis only rank " << rank <<
". Block size is " << blockSize_ <<
".\n"
884 <<
"Recovery failed." << std::endl;
888 if (rank < localsize) {
890 std::vector<int> redind(localsize);
891 for (
int i=0; i<localsize; i++) redind[i] = i;
893 restart = MVT::CloneViewNonConst(*restart,redind);
894 Krestart = MVT::CloneViewNonConst(*Krestart,redind);
902 Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
903 std::vector<MagnitudeType> theta(localsize);
907 MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
910 OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
913 MVT::MvTransMv(1.0,*restart,*Krestart,KK);
915 msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
916 if (rank < blockSize_) {
917 printer->stream(
Errors) <<
"Error! Recovered basis of rank " << rank <<
" produced only " << rank <<
"ritz vectors.\n"
918 <<
"Block size is " << blockSize_ <<
".\n"
919 <<
"Recovery failed." << std::endl;
926 Teuchos::BLAS<int,ScalarType> blas;
927 std::vector<int> order(rank);
929 sorter->sort( theta, Teuchos::rcpFromRef(order),rank );
931 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
932 msutils::permuteVectors(order,curS);
935 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
938 LOBPCGState<ScalarType,MV> newstate;
939 Teuchos::RCP<MV> newX;
941 std::vector<int> bsind(blockSize_);
942 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
943 newX = MVT::CloneViewNonConst(*Krestart,bsind);
945 MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
948 theta.resize(blockSize_);
949 newstate.T = Teuchos::rcpFromRef(theta);
951 lobpcg_solver->initialize(newstate);
955 <<
"Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
956 << err.what() << std::endl
957 <<
"Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
963 sol.numVecs = ordertest->howMany();
964 if (sol.numVecs > 0) {
965 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
966 sol.Espace = sol.Evecs;
967 sol.Evals.resize(sol.numVecs);
968 std::vector<MagnitudeType> vals(sol.numVecs);
971 std::vector<int> which = ordertest->whichVecs();
975 std::vector<int> inlocked(0), insolver(0);
976 for (
unsigned int i=0; i<which.size(); i++) {
978 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
979 insolver.push_back(which[i]);
983 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
984 inlocked.push_back(which[i] + curNumLocked);
988 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.numVecs,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
991 if (insolver.size() > 0) {
993 int lclnum = insolver.size();
994 std::vector<int> tosol(lclnum);
995 for (
int i=0; i<lclnum; i++) tosol[i] = i;
996 Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
997 MVT::SetBlock(*v,tosol,*sol.Evecs);
999 std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
1000 for (
unsigned int i=0; i<insolver.size(); i++) {
1001 vals[i] = fromsolver[insolver[i]].realpart;
1006 if (inlocked.size() > 0) {
1007 int solnum = insolver.size();
1009 int lclnum = inlocked.size();
1010 std::vector<int> tosol(lclnum);
1011 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1012 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1013 MVT::SetBlock(*v,tosol,*sol.Evecs);
1015 for (
unsigned int i=0; i<inlocked.size(); i++) {
1016 vals[i+solnum] = lockvals[inlocked[i]];
1022 std::vector<int> order(sol.numVecs);
1023 sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
1025 for (
int i=0; i<sol.numVecs; i++) {
1026 sol.Evals[i].realpart = vals[i];
1027 sol.Evals[i].imagpart = MT::zero();
1030 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1034 sol.index.resize(sol.numVecs,0);
1039 lobpcg_solver->currentStatus(printer->stream(
FinalSummary));
1042#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1044 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
1048 problem_->setSolution(sol);
1049 printer->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1052 numIters_ = lobpcg_solver->getNumIters();
1054 if (sol.numVecs < nev) {