00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP
00030 #define ANASAZI_LOBPCG_SOLMGR_HPP
00031
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038
00039 #include "AnasaziEigenproblem.hpp"
00040 #include "AnasaziSolverManager.hpp"
00041 #include "AnasaziSolverUtils.hpp"
00042
00043 #include "AnasaziLOBPCG.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziBasicOrthoManager.hpp"
00047 #include "AnasaziStatusTestMaxIters.hpp"
00048 #include "AnasaziStatusTestResNorm.hpp"
00049 #include "AnasaziStatusTestWithOrdering.hpp"
00050 #include "AnasaziStatusTestCombo.hpp"
00051 #include "AnasaziStatusTestOutput.hpp"
00052 #include "AnasaziBasicOutputManager.hpp"
00053 #include "Teuchos_BLAS.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055
00056
00070 namespace Anasazi {
00071
00072
00108 template<class ScalarType, class MV, class OP>
00109 class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00110
00111 private:
00112 typedef MultiVecTraits<ScalarType,MV> MVT;
00113 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00114 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00115 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00116 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00117
00118 public:
00119
00121
00122
00147 LOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00148 Teuchos::ParameterList &pl );
00149
00151 virtual ~LOBPCGSolMgr() {};
00153
00155
00156
00158 const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00159 return *problem_;
00160 }
00161
00163 int getNumIters() const {
00164 return numIters_;
00165 }
00166
00173 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00174 return tuple(_timerSolve, _timerLocking);
00175 }
00176
00177
00179
00181
00182
00209 ReturnType solve();
00210
00212 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
00213
00215 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
00216
00218 void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
00219
00221 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
00222
00224 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
00225
00227 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
00228
00230
00231 private:
00232 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00233
00234 std::string whch_, ortho_;
00235
00236 MagnitudeType convtol_, locktol_;
00237 int maxIters_, numIters_;
00238 bool useLocking_;
00239 bool relconvtol_, rellocktol_;
00240 int blockSize_;
00241 bool fullOrtho_;
00242 int maxLocked_;
00243 int verbosity_;
00244 int lockQuorum_;
00245 bool recover_;
00246 Teuchos::RCP<LOBPCGState<ScalarType,MV> > state_;
00247 enum StatusTestResNorm<ScalarType,MV,OP>::ResType convNorm_, lockNorm_;
00248
00249 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerLocking;
00250
00251 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
00252 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
00253 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
00254 };
00255
00256
00257
00258 template<class ScalarType, class MV, class OP>
00259 LOBPCGSolMgr<ScalarType,MV,OP>::LOBPCGSolMgr(
00260 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00261 Teuchos::ParameterList &pl ) :
00262 problem_(problem),
00263 whch_("SR"),
00264 ortho_("SVQB"),
00265 convtol_(MT::prec()),
00266 maxIters_(100),
00267 numIters_(0),
00268 useLocking_(false),
00269 relconvtol_(true),
00270 rellocktol_(true),
00271 blockSize_(0),
00272 fullOrtho_(true),
00273 maxLocked_(0),
00274 verbosity_(Anasazi::Errors),
00275 lockQuorum_(1),
00276 recover_(true),
00277 _timerSolve(Teuchos::TimeMonitor::getNewTimer("LOBPCGSolMgr::solve()")),
00278 _timerLocking(Teuchos::TimeMonitor::getNewTimer("LOBPCGSolMgr locking"))
00279 {
00280 TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00281 TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
00282 TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
00283 TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00284
00285
00286 std::string strtmp;
00287
00288
00289 whch_ = pl.get("Which",whch_);
00290 TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
00291 std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
00292
00293
00294 ortho_ = pl.get("Orthogonalization",ortho_);
00295 if (ortho_ != "DGKS" && ortho_ != "SVQB") {
00296 ortho_ = "SVQB";
00297 }
00298
00299
00300 convtol_ = pl.get("Convergence Tolerance",convtol_);
00301 relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00302 strtmp = pl.get("Convergence Norm",std::string("2"));
00303 if (strtmp == "2") {
00304 convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00305 }
00306 else if (strtmp == "M") {
00307 convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00308 }
00309 else {
00310 TEST_FOR_EXCEPTION(true, std::invalid_argument,
00311 "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
00312 }
00313
00314
00315
00316 useLocking_ = pl.get("Use Locking",useLocking_);
00317 rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00318
00319 locktol_ = convtol_/10;
00320 locktol_ = pl.get("Locking Tolerance",locktol_);
00321 strtmp = pl.get("Locking Norm",std::string("2"));
00322 if (strtmp == "2") {
00323 lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00324 }
00325 else if (strtmp == "M") {
00326 lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00327 }
00328 else {
00329 TEST_FOR_EXCEPTION(true, std::invalid_argument,
00330 "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
00331 }
00332
00333
00334 maxIters_ = pl.get("Maximum Iterations",maxIters_);
00335
00336
00337 blockSize_ = pl.get("Block Size",problem_->getNEV());
00338 TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00339 "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
00340
00341
00342 if (useLocking_) {
00343 maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00344 }
00345 else {
00346 maxLocked_ = 0;
00347 }
00348 if (maxLocked_ == 0) {
00349 useLocking_ = false;
00350 }
00351 TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00352 "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
00353 TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
00354 std::invalid_argument,
00355 "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
00356
00357 if (useLocking_) {
00358 lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00359 TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00360 std::invalid_argument,
00361 "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
00362 }
00363
00364
00365 fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
00366
00367
00368 if (pl.isParameter("Verbosity")) {
00369 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00370 verbosity_ = pl.get("Verbosity", verbosity_);
00371 } else {
00372 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00373 }
00374 }
00375
00376
00377 recover_ = pl.get("Recover",recover_);
00378
00379
00380 if (pl.isParameter("Init")) {
00381 state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
00382 }
00383 }
00384
00385
00386
00387 template<class ScalarType, class MV, class OP>
00388 ReturnType
00389 LOBPCGSolMgr<ScalarType,MV,OP>::solve() {
00390
00391 typedef SolverUtils<ScalarType,MV,OP> msutils;
00392
00393 const int nev = problem_->getNEV();
00394
00395
00396
00398
00399 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00400
00402
00403 Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
00404
00406
00407
00408
00409 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
00410 if (maxIters_ > 0) {
00411 maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00412 }
00413
00414 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00415 if (globalTest_ == Teuchos::null) {
00416 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00417 }
00418 else {
00419 convtest = globalTest_;
00420 }
00421 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
00422 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00423
00424 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
00425 if (useLocking_) {
00426 if (lockingTest_ == Teuchos::null) {
00427 locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
00428 }
00429 else {
00430 locktest = lockingTest_;
00431 }
00432 }
00433
00434 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00435 alltests.push_back(ordertest);
00436 if (locktest != Teuchos::null) alltests.push_back(locktest);
00437 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00438 if (maxtest != Teuchos::null) alltests.push_back(maxtest);
00439 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00440 = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00441
00442 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00443 if ( printer->isVerbosity(Debug) ) {
00444 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
00445 }
00446 else {
00447 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00448 }
00449
00451
00452 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
00453 if (ortho_=="SVQB") {
00454 ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00455 } else if (ortho_=="DGKS") {
00456 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00457 } else {
00458 TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
00459 }
00460
00462
00463 Teuchos::ParameterList plist;
00464 plist.set("Block Size",blockSize_);
00465 plist.set("Full Ortho",fullOrtho_);
00466
00468
00469 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
00470 = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00471
00472 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00473 if (probauxvecs != Teuchos::null) {
00474 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00475 }
00476
00478
00479
00480
00481 int curNumLocked = 0;
00482 Teuchos::RCP<MV> lockvecs;
00483 if (useLocking_) {
00484 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00485 }
00486 std::vector<MagnitudeType> lockvals;
00487
00488
00489
00490
00491
00492
00493
00494
00495 Teuchos::RCP<MV> workMV;
00496 if (fullOrtho_ == false && recover_ == true) {
00497 workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
00498 }
00499 else if (useLocking_) {
00500 if (problem_->getM() != Teuchos::null) {
00501 workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
00502 }
00503 else {
00504 workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
00505 }
00506 }
00507
00508
00509 Eigensolution<ScalarType,MV> sol;
00510 sol.numVecs = 0;
00511 problem_->setSolution(sol);
00512
00513
00514 if (state_ != Teuchos::null) {
00515 lobpcg_solver->initialize(*state_);
00516 }
00517
00518
00519 {
00520 Teuchos::TimeMonitor slvtimer(*_timerSolve);
00521
00522
00523 while (1) {
00524 try {
00525 lobpcg_solver->iterate();
00526
00528
00529
00530
00532 if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
00533 throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
00534 }
00536
00537
00538
00540 else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
00541
00542
00543
00544 break;
00545 }
00547
00548
00549
00551 else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00552
00553 Teuchos::TimeMonitor lcktimer(*_timerLocking);
00554
00555
00556 TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
00557 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
00558 TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
00559 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
00560 TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
00561 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
00562
00563 int numnew = locktest->howMany();
00564 std::vector<int> indnew = locktest->whichVecs();
00565
00566
00567 if (curNumLocked + numnew > maxLocked_) {
00568 numnew = maxLocked_ - curNumLocked;
00569 indnew.resize(numnew);
00570 }
00571
00572
00573
00574 bool hadP = lobpcg_solver->hasP();
00575
00576 {
00577
00578 printer->print(Debug,"Locking vectors: ");
00579 for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
00580 printer->print(Debug,"\n");
00581 }
00582 std::vector<MagnitudeType> newvals(numnew);
00583 Teuchos::RCP<const MV> newvecs;
00584 {
00585
00586
00587 newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
00588
00589 std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
00590 for (int i=0; i<numnew; i++) {
00591 newvals[i] = allvals[indnew[i]].realpart;
00592 }
00593 }
00594
00595 {
00596 std::vector<int> indlock(numnew);
00597 for (int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
00598 MVT::SetBlock(*newvecs,indlock,*lockvecs);
00599 newvecs = Teuchos::null;
00600 }
00601
00602 lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
00603 curNumLocked += numnew;
00604
00605 {
00606 std::vector<int> indlock(curNumLocked);
00607 for (int i=0; i<curNumLocked; i++) indlock[i] = i;
00608 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00609 if (probauxvecs != Teuchos::null) {
00610 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
00611 }
00612 else {
00613 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
00614 }
00615 }
00616
00617 ordertest->setAuxVals(lockvals);
00618
00619 {
00620 LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
00621 Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
00622
00623
00624
00625
00626 std::vector<int> bsind(blockSize_);
00627 for (int i=0; i<blockSize_; i++) bsind[i] = i;
00628 newstateX = MVT::CloneView(*workMV,bsind);
00629 MVT::SetBlock(*state.X,bsind,*newstateX);
00630
00631 if (state.MX != Teuchos::null) {
00632 std::vector<int> block3(blockSize_);
00633 for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
00634 newstateMX = MVT::CloneView(*workMV,block3);
00635 MVT::SetBlock(*state.MX,bsind,*newstateMX);
00636 }
00637
00638
00639 {
00640 Teuchos::RCP<MV> newX = MVT::CloneView(*newstateX,indnew);
00641 MVT::MvRandom(*newX);
00642
00643 if (newstateMX != Teuchos::null) {
00644 Teuchos::RCP<MV> newMX = MVT::CloneView(*newstateMX,indnew);
00645 OPT::Apply(*problem_->getM(),*newX,*newMX);
00646 }
00647 }
00648
00649 Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
00650 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
00651
00652 ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
00653
00654 if (hadP) {
00655
00656
00657 std::vector<int> block2(blockSize_);
00658 for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
00659 newstateP = MVT::CloneView(*workMV,block2);
00660 MVT::SetBlock(*state.P,bsind,*newstateP);
00661
00662 if (state.MP != Teuchos::null) {
00663 std::vector<int> block4(blockSize_);
00664 for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
00665 newstateMP = MVT::CloneView(*workMV,block4);
00666 MVT::SetBlock(*state.MP,bsind,*newstateMP);
00667 }
00668
00669 if (fullOrtho_) {
00670
00671 curauxvecs.push_back(newstateX);
00672 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00673 }
00674 else {
00675
00676 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00677 }
00678 }
00679
00680 LOBPCGState<ScalarType,MV> newstate;
00681 newstate.X = newstateX;
00682 newstate.MX = newstateMX;
00683 newstate.P = newstateP;
00684 newstate.MP = newstateMP;
00685 lobpcg_solver->initialize(newstate);
00686 }
00687
00688 if (curNumLocked == maxLocked_) {
00689
00690 combotest->removeTest(locktest);
00691 }
00692 }
00693 else {
00694 TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
00695 }
00696 }
00698
00699
00700
00702 catch (LOBPCGRitzFailure re) {
00703 if (fullOrtho_==true || recover_==false) {
00704
00705
00706
00707 printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
00708 << "Will not try to recover." << std::endl;
00709 break;
00710 }
00711 printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
00712 << "Full orthogonalization is off; will try to recover." << std::endl;
00713
00714
00715
00716
00717 LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
00718 Teuchos::RCP<MV> restart, Krestart, Mrestart;
00719 int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
00720 bool hasM = problem_->getM() != Teuchos::null;
00721 {
00722 std::vector<int> recind(localsize);
00723 for (int i=0; i<localsize; i++) recind[i] = i;
00724 restart = MVT::CloneView(*workMV,recind);
00725 }
00726 {
00727 std::vector<int> recind(localsize);
00728 for (int i=0; i<localsize; i++) recind[i] = localsize+i;
00729 Krestart = MVT::CloneView(*workMV,recind);
00730 }
00731 if (hasM) {
00732 Mrestart = Krestart;
00733 }
00734 else {
00735 Mrestart = restart;
00736 }
00737
00738
00739
00740
00741 {
00742 std::vector<int> blk1(blockSize_);
00743 for (int i=0; i < blockSize_; i++) blk1[i] = i;
00744 MVT::SetBlock(*curstate.X,blk1,*restart);
00745
00746
00747 if (hasM) {
00748 MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
00749 }
00750 }
00751
00752
00753 {
00754 std::vector<int> blk2(blockSize_);
00755 for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
00756 MVT::SetBlock(*curstate.H,blk2,*restart);
00757
00758
00759 if (hasM) {
00760 MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
00761 }
00762 }
00763
00764 if (localsize == 3*blockSize_) {
00765 std::vector<int> blk3(blockSize_);
00766 for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
00767 MVT::SetBlock(*curstate.P,blk3,*restart);
00768
00769
00770 if (hasM) {
00771 MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
00772 }
00773 }
00774
00775 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
00776 Teuchos::Array<Teuchos::RCP<const MV> > Q;
00777 {
00778 if (curNumLocked > 0) {
00779 std::vector<int> indlock(curNumLocked);
00780 for (int i=0; i<curNumLocked; i++) indlock[i] = i;
00781 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00782 Q.push_back(curlocked);
00783 }
00784 if (probauxvecs != Teuchos::null) {
00785 Q.push_back(probauxvecs);
00786 }
00787 }
00788 int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
00789 if (rank < blockSize_) {
00790
00791 printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
00792 << "Recovery failed." << std::endl;
00793 break;
00794 }
00795
00796 if (rank < localsize) {
00797 localsize = rank;
00798 std::vector<int> redind(localsize);
00799 for (int i=0; i<localsize; i++) redind[i] = i;
00800
00801 restart = MVT::CloneView(*restart,redind);
00802 Krestart = MVT::CloneView(*Krestart,redind);
00803 if (hasM) {
00804 Mrestart = Krestart;
00805 }
00806 else {
00807 Mrestart = restart;
00808 }
00809 }
00810 Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
00811 std::vector<MagnitudeType> theta(localsize);
00812
00813
00814
00815 MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
00816
00817
00818 OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
00819
00820
00821 MVT::MvTransMv(1.0,*restart,*Krestart,KK);
00822 rank = localsize;
00823 msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
00824 if (rank < blockSize_) {
00825 printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
00826 << "Block size is " << blockSize_ << ".\n"
00827 << "Recovery failed." << std::endl;
00828 break;
00829 }
00830 theta.resize(rank);
00831
00832
00833 {
00834 Teuchos::BLAS<int,ScalarType> blas;
00835 std::vector<int> order(rank);
00836
00837 sorter->sort( theta, Teuchos::rcpFromRef(order),rank );
00838
00839 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
00840 msutils::permuteVectors(order,curS);
00841 }
00842
00843 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
00844
00845
00846 LOBPCGState<ScalarType,MV> newstate;
00847 Teuchos::RCP<MV> newX;
00848 {
00849 std::vector<int> bsind(blockSize_);
00850 for (int i=0; i<blockSize_; i++) bsind[i] = i;
00851 newX = MVT::CloneView(*Krestart,bsind);
00852 }
00853 MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
00854
00855 newstate.X = newX;
00856 theta.resize(blockSize_);
00857 newstate.T = Teuchos::rcpFromRef(theta);
00858
00859 lobpcg_solver->initialize(newstate);
00860 }
00861 catch (AnasaziError err) {
00862 printer->stream(Errors)
00863 << "Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
00864 << err.what() << std::endl
00865 << "Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
00866 return Unconverged;
00867 }
00868
00869 }
00870
00871 sol.numVecs = ordertest->howMany();
00872 if (sol.numVecs > 0) {
00873 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
00874 sol.Espace = sol.Evecs;
00875 sol.Evals.resize(sol.numVecs);
00876 std::vector<MagnitudeType> vals(sol.numVecs);
00877
00878
00879 std::vector<int> which = ordertest->whichVecs();
00880
00881
00882
00883 std::vector<int> inlocked(0), insolver(0);
00884 for (unsigned int i=0; i<which.size(); i++) {
00885 if (which[i] >= 0) {
00886 TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
00887 insolver.push_back(which[i]);
00888 }
00889 else {
00890
00891 TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
00892 inlocked.push_back(which[i] + curNumLocked);
00893 }
00894 }
00895
00896 TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
00897
00898
00899 if (insolver.size() > 0) {
00900
00901 int lclnum = insolver.size();
00902 std::vector<int> tosol(lclnum);
00903 for (int i=0; i<lclnum; i++) tosol[i] = i;
00904 Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
00905 MVT::SetBlock(*v,tosol,*sol.Evecs);
00906
00907 std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
00908 for (unsigned int i=0; i<insolver.size(); i++) {
00909 vals[i] = fromsolver[insolver[i]].realpart;
00910 }
00911 }
00912
00913
00914 if (inlocked.size() > 0) {
00915 int solnum = insolver.size();
00916
00917 int lclnum = inlocked.size();
00918 std::vector<int> tosol(lclnum);
00919 for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
00920 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
00921 MVT::SetBlock(*v,tosol,*sol.Evecs);
00922
00923 for (unsigned int i=0; i<inlocked.size(); i++) {
00924 vals[i+solnum] = lockvals[inlocked[i]];
00925 }
00926 }
00927
00928
00929 {
00930 std::vector<int> order(sol.numVecs);
00931 sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
00932
00933 for (int i=0; i<sol.numVecs; i++) {
00934 sol.Evals[i].realpart = vals[i];
00935 sol.Evals[i].imagpart = MT::zero();
00936 }
00937
00938 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
00939 }
00940
00941
00942 sol.index.resize(sol.numVecs,0);
00943 }
00944 }
00945
00946
00947 lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00948
00949
00950 Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00951
00952 problem_->setSolution(sol);
00953 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00954
00955
00956 numIters_ = lobpcg_solver->getNumIters();
00957
00958 if (sol.numVecs < nev) {
00959 return Unconverged;
00960 }
00961 return Converged;
00962 }
00963
00964
00965 template <class ScalarType, class MV, class OP>
00966 void
00967 LOBPCGSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
00968 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
00969 {
00970 globalTest_ = global;
00971 }
00972
00973 template <class ScalarType, class MV, class OP>
00974 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00975 LOBPCGSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const
00976 {
00977 return globalTest_;
00978 }
00979
00980 template <class ScalarType, class MV, class OP>
00981 void
00982 LOBPCGSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
00983 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
00984 {
00985 debugTest_ = debug;
00986 }
00987
00988 template <class ScalarType, class MV, class OP>
00989 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00990 LOBPCGSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
00991 {
00992 return debugTest_;
00993 }
00994
00995 template <class ScalarType, class MV, class OP>
00996 void
00997 LOBPCGSolMgr<ScalarType,MV,OP>::setLockingStatusTest(
00998 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
00999 {
01000 lockingTest_ = locking;
01001 }
01002
01003 template <class ScalarType, class MV, class OP>
01004 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
01005 LOBPCGSolMgr<ScalarType,MV,OP>::getLockingStatusTest() const
01006 {
01007 return lockingTest_;
01008 }
01009
01010 }
01011
01012 #endif