Anasazi Version of the Day
AnasaziLOBPCGSolMgr.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
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 Teuchos::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 // Constructor
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 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00278   , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr::solve()")),
00279   _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr locking"))
00280 #endif
00281 {
00282   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00283   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00284   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00285   TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00286 
00287 
00288   std::string strtmp;
00289 
00290   // which values to solve for
00291   whch_ = pl.get("Which",whch_);
00292   TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
00293       std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
00294 
00295   // which orthogonalization to use
00296   ortho_ = pl.get("Orthogonalization",ortho_);
00297   if (ortho_ != "DGKS" && ortho_ != "SVQB") {
00298     ortho_ = "SVQB";
00299   }
00300 
00301   // convergence tolerance
00302   convtol_ = pl.get("Convergence Tolerance",convtol_);
00303   relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00304   strtmp = pl.get("Convergence Norm",std::string("2"));
00305   if (strtmp == "2") {
00306     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00307   }
00308   else if (strtmp == "M") {
00309     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00310   }
00311   else {
00312     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, 
00313         "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
00314   }
00315 
00316   
00317   // locking tolerance
00318   useLocking_ = pl.get("Use Locking",useLocking_);
00319   rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00320   // default: should be less than convtol_
00321   locktol_ = convtol_/10;
00322   locktol_ = pl.get("Locking Tolerance",locktol_);
00323   strtmp = pl.get("Locking Norm",std::string("2"));
00324   if (strtmp == "2") {
00325     lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00326   }
00327   else if (strtmp == "M") {
00328     lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00329   }
00330   else {
00331     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, 
00332         "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
00333   }
00334 
00335   // maximum number of iterations
00336   maxIters_ = pl.get("Maximum Iterations",maxIters_);
00337 
00338   // block size: default is nev()
00339   blockSize_ = pl.get("Block Size",problem_->getNEV());
00340   TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00341                      "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
00342 
00343   // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
00344   if (useLocking_) {
00345     maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00346   }
00347   else {
00348     maxLocked_ = 0;
00349   }
00350   if (maxLocked_ == 0) {
00351     useLocking_ = false;
00352   }
00353   TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00354                      "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
00355   TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(), 
00356                      std::invalid_argument,
00357                      "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
00358 
00359   if (useLocking_) {
00360     lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00361     TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00362                        std::invalid_argument,
00363                        "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
00364   }
00365 
00366   // full orthogonalization: default true
00367   fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
00368 
00369   // verbosity level
00370   if (pl.isParameter("Verbosity")) {
00371     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00372       verbosity_ = pl.get("Verbosity", verbosity_);
00373     } else {
00374       verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00375     }
00376   }
00377 
00378   // recover from LOBPCGRitzFailure
00379   recover_ = pl.get("Recover",recover_);
00380 
00381   // get (optionally) an initial state
00382   if (pl.isParameter("Init")) {
00383     state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
00384   }
00385 }
00386 
00387 
00388 // solve()
00389 template<class ScalarType, class MV, class OP>
00390 ReturnType 
00391 LOBPCGSolMgr<ScalarType,MV,OP>::solve() {
00392 
00393   typedef SolverUtils<ScalarType,MV,OP> msutils;
00394 
00395   const int nev = problem_->getNEV();
00396 
00397 
00398 
00400   // Sort manager
00401   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00402 
00404   // Output manager
00405   Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
00406 
00408   // Status tests
00409   //
00410   // maximum number of iterations: optional test
00411   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
00412   if (maxIters_ > 0) {
00413     maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00414   }
00415   // convergence
00416   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00417   if (globalTest_ == Teuchos::null) {
00418     convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00419   }
00420   else {
00421     convtest = globalTest_;
00422   }
00423   Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest 
00424     = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00425   // locking
00426   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
00427   if (useLocking_) {
00428     if (lockingTest_ == Teuchos::null) {
00429       locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
00430     }
00431     else {
00432       locktest = lockingTest_;
00433     }
00434   }
00435   // for a non-short-circuited OR test, the order doesn't matter
00436   Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00437   alltests.push_back(ordertest);
00438   if (locktest != Teuchos::null) alltests.push_back(locktest);
00439   if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00440   if (maxtest != Teuchos::null) alltests.push_back(maxtest);
00441   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00442     = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00443   // printing StatusTest
00444   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00445   if ( printer->isVerbosity(Debug) ) {
00446     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
00447   }
00448   else {
00449     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00450   }
00451 
00453   // Orthomanager
00454   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho; 
00455   if (ortho_=="SVQB") {
00456     ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00457   } else if (ortho_=="DGKS") {
00458     ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00459   } else {
00460     TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
00461   }
00462 
00464   // Parameter list
00465   Teuchos::ParameterList plist;
00466   plist.set("Block Size",blockSize_);
00467   plist.set("Full Ortho",fullOrtho_);
00468 
00470   // LOBPCG solver
00471   Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver 
00472     = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00473   // set any auxiliary vectors defined in the problem
00474   Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00475   if (probauxvecs != Teuchos::null) {
00476     lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00477   }
00478 
00480   // Storage
00481   // 
00482   // lockvecs will contain eigenvectors that have been determined "locked" by the status test
00483   int curNumLocked = 0;
00484   Teuchos::RCP<MV> lockvecs;
00485   if (useLocking_) {
00486     lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00487   }
00488   std::vector<MagnitudeType> lockvals;
00489   // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
00490   // it will be partitioned in these cases as follows:
00491   // for LOBPCGRitzFailure recovery:
00492   // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
00493   // total size: 2*3*blocksize
00494   // for locking
00495   // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
00496   // total size: 2*blocksize or 4*blocksize
00497   Teuchos::RCP<MV> workMV;
00498   if (fullOrtho_ == false && recover_ == true) {
00499     workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
00500   }
00501   else if (useLocking_) {
00502     if (problem_->getM() != Teuchos::null) {
00503       workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
00504     }
00505     else {
00506       workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
00507     }
00508   }
00509 
00510   // initialize the solution to nothing in case we throw an exception
00511   Eigensolution<ScalarType,MV> sol;
00512   sol.numVecs = 0;
00513   problem_->setSolution(sol);
00514 
00515   // initialize the solver if the user specified a state
00516   if (state_ != Teuchos::null) {
00517     lobpcg_solver->initialize(*state_);
00518   }
00519 
00520   // enter solve() iterations
00521   {
00522 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00523     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00524 #endif
00525 
00526     // tell the lobpcg_solver to iterate
00527     while (1) {
00528       try {
00529         lobpcg_solver->iterate();
00530 
00532         //
00533         // check user-specified debug test; if it passed, return an exception
00534         //
00536         if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
00537           throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
00538         }
00540         //
00541         // check convergence first
00542         //
00544         else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
00545           // we have convergence or not
00546           // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
00547           // ordertest->howMany() will tell us how many
00548           break;
00549         }
00551         //
00552         // check locking if we didn't converge
00553         //
00555         else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00556 
00557 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00558           Teuchos::TimeMonitor lcktimer(*_timerLocking);
00559 #endif
00560 
00561           // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
00562           TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
00563               "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
00564           TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
00565               "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
00566           TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
00567               "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
00568           // get the indices
00569           int numnew = locktest->howMany();
00570           std::vector<int> indnew = locktest->whichVecs();
00571 
00572           // don't lock more than maxLocked_; we didn't allocate enough space.
00573           if (curNumLocked + numnew > maxLocked_) {
00574             numnew = maxLocked_ - curNumLocked;
00575             indnew.resize(numnew);
00576           }
00577 
00578           // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
00579           // store the hasP() state for use below
00580           bool hadP = lobpcg_solver->hasP();
00581 
00582           {
00583             // debug printing
00584             printer->print(Debug,"Locking vectors: ");
00585             for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
00586             printer->print(Debug,"\n");
00587           }
00588           std::vector<MagnitudeType> newvals(numnew);
00589           Teuchos::RCP<const MV> newvecs;
00590           {
00591             // work in a local scope, to hide the variabes needed for extracting this info
00592             // get the vectors
00593             newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
00594             // get the values
00595             std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
00596             for (int i=0; i<numnew; i++) {
00597               newvals[i] = allvals[indnew[i]].realpart;
00598             }
00599           }
00600           // put newvecs into lockvecs
00601           {
00602             std::vector<int> indlock(numnew);
00603             for (int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
00604             MVT::SetBlock(*newvecs,indlock,*lockvecs);
00605             newvecs = Teuchos::null;
00606           }
00607           // put newvals into lockvals
00608           lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
00609           curNumLocked += numnew;
00610           // add locked vecs as aux vecs, along with aux vecs from problem
00611           {
00612             std::vector<int> indlock(curNumLocked);
00613             for (int i=0; i<curNumLocked; i++) indlock[i] = i;
00614             Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00615             if (probauxvecs != Teuchos::null) {
00616               lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
00617             }
00618             else {
00619               lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
00620             }
00621           }
00622           // add locked vals to ordertest
00623           ordertest->setAuxVals(lockvals);
00624           // fill out the empty state in the solver
00625           {
00626             LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
00627             Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
00628             //
00629             // workMV will be partitioned as follows: workMV = [X P MX MP], 
00630             //
00631             // make a copy of the current X,MX state
00632             std::vector<int> bsind(blockSize_); 
00633             for (int i=0; i<blockSize_; i++) bsind[i] = i;
00634             newstateX = MVT::CloneViewNonConst(*workMV,bsind);
00635             MVT::SetBlock(*state.X,bsind,*newstateX);
00636 
00637             if (state.MX != Teuchos::null) {
00638               std::vector<int> block3(blockSize_);
00639               for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
00640               newstateMX = MVT::CloneViewNonConst(*workMV,block3);
00641               MVT::SetBlock(*state.MX,bsind,*newstateMX);
00642             }
00643             //
00644             // get select part, set to random, apply M
00645             {
00646               Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
00647               MVT::MvRandom(*newX);
00648 
00649               if (newstateMX != Teuchos::null) {
00650                 Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
00651                 OPT::Apply(*problem_->getM(),*newX,*newMX);
00652               }
00653             }
00654 
00655             Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
00656             Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
00657             // ortho X against the aux vectors
00658             ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
00659 
00660             if (hadP) {
00661               //
00662               // get P and optionally MP, orthogonalize against X and auxiliary vectors
00663               std::vector<int> block2(blockSize_);
00664               for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
00665               newstateP = MVT::CloneViewNonConst(*workMV,block2);
00666               MVT::SetBlock(*state.P,bsind,*newstateP);
00667 
00668               if (state.MP != Teuchos::null) {
00669                 std::vector<int> block4(blockSize_);
00670                 for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
00671                 newstateMP = MVT::CloneViewNonConst(*workMV,block4);
00672                 MVT::SetBlock(*state.MP,bsind,*newstateMP);
00673               }
00674 
00675               if (fullOrtho_) {
00676                 // ortho P against the new aux vectors and new X
00677                 curauxvecs.push_back(newstateX);
00678                 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00679               }
00680               else {
00681                 // ortho P against the new aux vectors
00682                 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00683               }
00684             }
00685             // set the new state
00686             LOBPCGState<ScalarType,MV> newstate;
00687             newstate.X  = newstateX;
00688             newstate.MX = newstateMX;
00689             newstate.P  = newstateP;
00690             newstate.MP = newstateMP;
00691             lobpcg_solver->initialize(newstate);
00692           }
00693 
00694           if (curNumLocked == maxLocked_) {
00695             // disable locking now; remove locking test from combo test
00696             combotest->removeTest(locktest);
00697           }
00698         }
00699         else {
00700           TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
00701         }
00702       }
00704       //
00705       // check Ritz Failure
00706       //
00708       catch (const LOBPCGRitzFailure &re) {
00709         if (fullOrtho_==true || recover_==false) {
00710           // if we are already using full orthogonalization, there isn't much we can do here. 
00711           // the most recent information in the status tests is still valid, and can be used to extract/return the 
00712           // eigenpairs that have converged.
00713           printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
00714             << "Will not try to recover." << std::endl;
00715           break; // while(1)
00716         }
00717         printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
00718           << "Full orthogonalization is off; will try to recover." << std::endl;
00719         // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
00720         // if there aren't enough, break and quit with what we have
00721         //
00722         // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
00723         LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
00724         Teuchos::RCP<MV> restart, Krestart, Mrestart;
00725         int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
00726         bool hasM = problem_->getM() != Teuchos::null;
00727         {
00728           std::vector<int> recind(localsize);
00729           for (int i=0; i<localsize; i++) recind[i] = i;
00730           restart = MVT::CloneViewNonConst(*workMV,recind);
00731         }
00732         {
00733           std::vector<int> recind(localsize);
00734           for (int i=0; i<localsize; i++) recind[i] = localsize+i;
00735           Krestart = MVT::CloneViewNonConst(*workMV,recind);
00736         }
00737         if (hasM) {
00738           Mrestart = Krestart;
00739         }
00740         else {
00741           Mrestart = restart;
00742         }
00743         //
00744         // set restart = [X H P] and Mrestart = M*[X H P]
00745         //
00746         // put X into [0 , blockSize)
00747         {
00748           std::vector<int> blk1(blockSize_);
00749           for (int i=0; i < blockSize_; i++) blk1[i] = i;
00750           MVT::SetBlock(*curstate.X,blk1,*restart);
00751 
00752           // put MX into [0 , blockSize)
00753           if (hasM) {
00754             MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
00755           }
00756         }
00757         //
00758         // put H into [blockSize_ , 2*blockSize)
00759         {
00760           std::vector<int> blk2(blockSize_);
00761           for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
00762           MVT::SetBlock(*curstate.H,blk2,*restart);
00763 
00764           // put MH into [blockSize_ , 2*blockSize)
00765           if (hasM) {
00766             MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
00767           }
00768         }
00769         // optionally, put P into [2*blockSize,3*blockSize)
00770         if (localsize == 3*blockSize_) {
00771           std::vector<int> blk3(blockSize_);
00772           for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
00773           MVT::SetBlock(*curstate.P,blk3,*restart);
00774 
00775           // put MP into [2*blockSize,3*blockSize)
00776           if (hasM) {
00777             MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
00778           }
00779         }
00780         // project against auxvecs and locked vecs, and orthonormalize the basis
00781         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
00782         Teuchos::Array<Teuchos::RCP<const MV> > Q;
00783         {
00784           if (curNumLocked > 0) {
00785             std::vector<int> indlock(curNumLocked);
00786             for (int i=0; i<curNumLocked; i++) indlock[i] = i;
00787             Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00788             Q.push_back(curlocked);
00789           }
00790           if (probauxvecs != Teuchos::null) {
00791             Q.push_back(probauxvecs);
00792           }
00793         }
00794         int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
00795         if (rank < blockSize_) {
00796           // quit
00797           printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
00798             << "Recovery failed." << std::endl;
00799           break;
00800         }
00801         // reduce multivec size if necessary
00802         if (rank < localsize) {
00803           localsize = rank;
00804           std::vector<int> redind(localsize);
00805           for (int i=0; i<localsize; i++) redind[i] = i;
00806           // grab the first part of restart,Krestart
00807           restart = MVT::CloneViewNonConst(*restart,redind);
00808           Krestart = MVT::CloneViewNonConst(*Krestart,redind);
00809           if (hasM) {
00810             Mrestart = Krestart;
00811           }
00812           else {
00813             Mrestart = restart;
00814           }
00815         }
00816         Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
00817         std::vector<MagnitudeType> theta(localsize);
00818         // project the matrices
00819         //
00820         // MM = restart^H M restart
00821         MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
00822         // 
00823         // compute Krestart = K*restart
00824         OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
00825         //
00826         // KK = restart^H K restart
00827         MVT::MvTransMv(1.0,*restart,*Krestart,KK);
00828         rank = localsize;
00829         msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
00830         if (rank < blockSize_) {
00831           printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
00832             << "Block size is " << blockSize_ << ".\n"
00833             << "Recovery failed." << std::endl;
00834           break;
00835         }
00836         theta.resize(rank);
00837         //
00838         // sort the ritz values using the sort manager
00839         {
00840           Teuchos::BLAS<int,ScalarType> blas;
00841           std::vector<int> order(rank);
00842           // sort
00843           sorter->sort( theta, Teuchos::rcpFromRef(order),rank );   // don't catch exception
00844           // Sort the primitive ritz vectors
00845           Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
00846           msutils::permuteVectors(order,curS);
00847         }
00848         //
00849         Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
00850         //
00851         // compute the ritz vectors: store them in Krestart
00852         LOBPCGState<ScalarType,MV> newstate;
00853         Teuchos::RCP<MV> newX; 
00854         {
00855           std::vector<int> bsind(blockSize_);
00856           for (int i=0; i<blockSize_; i++) bsind[i] = i;
00857           newX = MVT::CloneViewNonConst(*Krestart,bsind);
00858         }
00859         MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
00860         // send X and theta into the solver
00861         newstate.X = newX;
00862         theta.resize(blockSize_);
00863         newstate.T = Teuchos::rcpFromRef(theta);
00864         // initialize
00865         lobpcg_solver->initialize(newstate);
00866       }
00867       catch (const AnasaziError &err) {
00868         printer->stream(Errors) 
00869           << "Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
00870           << err.what() << std::endl
00871           << "Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
00872         return Unconverged;
00873       }
00874       // don't catch any other exceptions
00875     }
00876 
00877     sol.numVecs = ordertest->howMany();
00878     if (sol.numVecs > 0) {
00879       sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
00880       sol.Espace = sol.Evecs;
00881       sol.Evals.resize(sol.numVecs);
00882       std::vector<MagnitudeType> vals(sol.numVecs);
00883 
00884       // copy them into the solution
00885       std::vector<int> which = ordertest->whichVecs();
00886       // indices between [0,blockSize) refer to vectors/values in the solver
00887       // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
00888       // everything has already been ordered by the solver; we just have to partition the two references
00889       std::vector<int> inlocked(0), insolver(0);
00890       for (unsigned int i=0; i<which.size(); i++) {
00891         if (which[i] >= 0) {
00892           TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
00893           insolver.push_back(which[i]);
00894         }
00895         else {
00896           // sanity check
00897           TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
00898           inlocked.push_back(which[i] + curNumLocked);
00899         }
00900       }
00901 
00902       TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
00903 
00904       // set the vecs,vals in the solution
00905       if (insolver.size() > 0) {
00906         // set vecs
00907         int lclnum = insolver.size();
00908         std::vector<int> tosol(lclnum);
00909         for (int i=0; i<lclnum; i++) tosol[i] = i;
00910         Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
00911         MVT::SetBlock(*v,tosol,*sol.Evecs);
00912         // set vals
00913         std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
00914         for (unsigned int i=0; i<insolver.size(); i++) {
00915           vals[i] = fromsolver[insolver[i]].realpart;
00916         }
00917       }
00918 
00919       // get the vecs,vals from locked storage
00920       if (inlocked.size() > 0) {
00921         int solnum = insolver.size();
00922         // set vecs
00923         int lclnum = inlocked.size();
00924         std::vector<int> tosol(lclnum);
00925         for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
00926         Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
00927         MVT::SetBlock(*v,tosol,*sol.Evecs);
00928         // set vals
00929         for (unsigned int i=0; i<inlocked.size(); i++) {
00930           vals[i+solnum] = lockvals[inlocked[i]];
00931         }
00932       }
00933 
00934       // sort the eigenvalues and permute the eigenvectors appropriately
00935       {
00936         std::vector<int> order(sol.numVecs);
00937         sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
00938         // store the values in the Eigensolution
00939         for (int i=0; i<sol.numVecs; i++) {
00940           sol.Evals[i].realpart = vals[i];
00941           sol.Evals[i].imagpart = MT::zero();
00942         }
00943         // now permute the eigenvectors according to order
00944         msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
00945       }
00946 
00947       // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
00948       sol.index.resize(sol.numVecs,0);
00949     }
00950   }
00951 
00952   // print final summary
00953   lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00954 
00955   // print timing information
00956 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00957   if ( printer->isVerbosity( TimingDetails ) ) {
00958     Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
00959   }
00960 #endif
00961 
00962   problem_->setSolution(sol);
00963   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00964 
00965   // get the number of iterations performed in this call to solve.
00966   numIters_ = lobpcg_solver->getNumIters();
00967 
00968   if (sol.numVecs < nev) {
00969     return Unconverged; // return from LOBPCGSolMgr::solve() 
00970   }
00971   return Converged; // return from LOBPCGSolMgr::solve() 
00972 }
00973 
00974 
00975 template <class ScalarType, class MV, class OP>
00976 void 
00977 LOBPCGSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
00978     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global) 
00979 {
00980   globalTest_ = global;
00981 }
00982 
00983 template <class ScalarType, class MV, class OP>
00984 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
00985 LOBPCGSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const 
00986 {
00987   return globalTest_;
00988 }
00989 
00990 template <class ScalarType, class MV, class OP>
00991 void 
00992 LOBPCGSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
00993     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
00994 {
00995   debugTest_ = debug;
00996 }
00997 
00998 template <class ScalarType, class MV, class OP>
00999 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
01000 LOBPCGSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
01001 {
01002   return debugTest_;
01003 }
01004 
01005 template <class ScalarType, class MV, class OP>
01006 void 
01007 LOBPCGSolMgr<ScalarType,MV,OP>::setLockingStatusTest(
01008     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking) 
01009 {
01010   lockingTest_ = locking;
01011 }
01012 
01013 template <class ScalarType, class MV, class OP>
01014 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
01015 LOBPCGSolMgr<ScalarType,MV,OP>::getLockingStatusTest() const 
01016 {
01017   return lockingTest_;
01018 }
01019 
01020 } // end Anasazi namespace
01021 
01022 #endif /* ANASAZI_LOBPCG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends