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 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   _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   // which values to solve for
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   // which orthogonalization to use
00294   ortho_ = pl.get("Orthogonalization",ortho_);
00295   if (ortho_ != "DGKS" && ortho_ != "SVQB") {
00296     ortho_ = "SVQB";
00297   }
00298 
00299   // convergence tolerance
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   // locking tolerance
00316   useLocking_ = pl.get("Use Locking",useLocking_);
00317   rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00318   // default: should be less than convtol_
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   // maximum number of iterations
00334   maxIters_ = pl.get("Maximum Iterations",maxIters_);
00335 
00336   // block size: default is nev()
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   // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
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   // full orthogonalization: default true
00365   fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
00366 
00367   // verbosity level
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   // recover from LOBPCGRitzFailure
00377   recover_ = pl.get("Recover",recover_);
00378 
00379   // get (optionally) an initial state
00380   if (pl.isParameter("Init")) {
00381     state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
00382   }
00383 }
00384 
00385 
00386 // solve()
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   // Sort manager
00399   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00400 
00402   // Output manager
00403   Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
00404 
00406   // Status tests
00407   //
00408   // maximum number of iterations: optional test
00409   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
00410   if (maxIters_ > 0) {
00411     maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00412   }
00413   // convergence
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   // locking
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   // for a non-short-circuited OR test, the order doesn't matter
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   // printing StatusTest
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   // Orthomanager
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   // Parameter list
00463   Teuchos::ParameterList plist;
00464   plist.set("Block Size",blockSize_);
00465   plist.set("Full Ortho",fullOrtho_);
00466 
00468   // LOBPCG solver
00469   Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver 
00470     = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00471   // set any auxiliary vectors defined in the problem
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   // Storage
00479   // 
00480   // lockvecs will contain eigenvectors that have been determined "locked" by the status test
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   // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
00488   // it will be partitioned in these cases as follows:
00489   // for LOBPCGRitzFailure recovery:
00490   // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
00491   // total size: 2*3*blocksize
00492   // for locking
00493   // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
00494   // total size: 2*blocksize or 4*blocksize
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   // initialize the solution to nothing in case we throw an exception
00509   Eigensolution<ScalarType,MV> sol;
00510   sol.numVecs = 0;
00511   problem_->setSolution(sol);
00512 
00513   // initialize the solver if the user specified a state
00514   if (state_ != Teuchos::null) {
00515     lobpcg_solver->initialize(*state_);
00516   }
00517 
00518   // enter solve() iterations
00519   {
00520     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00521 
00522     // tell the lobpcg_solver to iterate
00523     while (1) {
00524       try {
00525         lobpcg_solver->iterate();
00526 
00528         //
00529         // check user-specified debug test; if it passed, return an exception
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         // check convergence first
00538         //
00540         else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
00541           // we have convergence or not
00542           // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
00543           // ordertest->howMany() will tell us how many
00544           break;
00545         }
00547         //
00548         // check locking if we didn't converge
00549         //
00551         else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00552 
00553           Teuchos::TimeMonitor lcktimer(*_timerLocking);
00554 
00555           // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
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           // get the indices
00563           int numnew = locktest->howMany();
00564           std::vector<int> indnew = locktest->whichVecs();
00565 
00566           // don't lock more than maxLocked_; we didn't allocate enough space.
00567           if (curNumLocked + numnew > maxLocked_) {
00568             numnew = maxLocked_ - curNumLocked;
00569             indnew.resize(numnew);
00570           }
00571 
00572           // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
00573           // store the hasP() state for use below
00574           bool hadP = lobpcg_solver->hasP();
00575 
00576           {
00577             // debug printing
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             // work in a local scope, to hide the variabes needed for extracting this info
00586             // get the vectors
00587             newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
00588             // get the values
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           // put newvecs into lockvecs
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           // put newvals into lockvals
00602           lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
00603           curNumLocked += numnew;
00604           // add locked vecs as aux vecs, along with aux vecs from problem
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           // add locked vals to ordertest
00617           ordertest->setAuxVals(lockvals);
00618           // fill out the empty state in the solver
00619           {
00620             LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
00621             Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
00622             //
00623             // workMV will be partitioned as follows: workMV = [X P MX MP], 
00624             //
00625             // make a copy of the current X,MX state
00626             std::vector<int> bsind(blockSize_); 
00627             for (int i=0; i<blockSize_; i++) bsind[i] = i;
00628             newstateX = MVT::CloneViewNonConst(*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::CloneViewNonConst(*workMV,block3);
00635               MVT::SetBlock(*state.MX,bsind,*newstateMX);
00636             }
00637             //
00638             // get select part, set to random, apply M
00639             {
00640               Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
00641               MVT::MvRandom(*newX);
00642 
00643               if (newstateMX != Teuchos::null) {
00644                 Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*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             // ortho X against the aux vectors
00652             ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
00653 
00654             if (hadP) {
00655               //
00656               // get P and optionally MP, orthogonalize against X and auxiliary vectors
00657               std::vector<int> block2(blockSize_);
00658               for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
00659               newstateP = MVT::CloneViewNonConst(*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::CloneViewNonConst(*workMV,block4);
00666                 MVT::SetBlock(*state.MP,bsind,*newstateMP);
00667               }
00668 
00669               if (fullOrtho_) {
00670                 // ortho P against the new aux vectors and new X
00671                 curauxvecs.push_back(newstateX);
00672                 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00673               }
00674               else {
00675                 // ortho P against the new aux vectors
00676                 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
00677               }
00678             }
00679             // set the new state
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             // disable locking now; remove locking test from combo test
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       // check Ritz Failure
00700       //
00702       catch (const LOBPCGRitzFailure &re) {
00703         if (fullOrtho_==true || recover_==false) {
00704           // if we are already using full orthogonalization, there isn't much we can do here. 
00705           // the most recent information in the status tests is still valid, and can be used to extract/return the 
00706           // eigenpairs that have converged.
00707           printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
00708             << "Will not try to recover." << std::endl;
00709           break; // while(1)
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         // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
00714         // if there aren't enough, break and quit with what we have
00715         //
00716         // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
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::CloneViewNonConst(*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::CloneViewNonConst(*workMV,recind);
00730         }
00731         if (hasM) {
00732           Mrestart = Krestart;
00733         }
00734         else {
00735           Mrestart = restart;
00736         }
00737         //
00738         // set restart = [X H P] and Mrestart = M*[X H P]
00739         //
00740         // put X into [0 , blockSize)
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           // put MX into [0 , blockSize)
00747           if (hasM) {
00748             MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
00749           }
00750         }
00751         //
00752         // put H into [blockSize_ , 2*blockSize)
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           // put MH into [blockSize_ , 2*blockSize)
00759           if (hasM) {
00760             MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
00761           }
00762         }
00763         // optionally, put P into [2*blockSize,3*blockSize)
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           // put MP into [2*blockSize,3*blockSize)
00770           if (hasM) {
00771             MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
00772           }
00773         }
00774         // project against auxvecs and locked vecs, and orthonormalize the basis
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           // quit
00791           printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
00792             << "Recovery failed." << std::endl;
00793           break;
00794         }
00795         // reduce multivec size if necessary
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           // grab the first part of restart,Krestart
00801           restart = MVT::CloneViewNonConst(*restart,redind);
00802           Krestart = MVT::CloneViewNonConst(*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         // project the matrices
00813         //
00814         // MM = restart^H M restart
00815         MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
00816         // 
00817         // compute Krestart = K*restart
00818         OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
00819         //
00820         // KK = restart^H K restart
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         // sort the ritz values using the sort manager
00833         {
00834           Teuchos::BLAS<int,ScalarType> blas;
00835           std::vector<int> order(rank);
00836           // sort
00837           sorter->sort( theta, Teuchos::rcpFromRef(order),rank );   // don't catch exception
00838           // Sort the primitive ritz vectors
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         // compute the ritz vectors: store them in Krestart
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::CloneViewNonConst(*Krestart,bsind);
00852         }
00853         MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
00854         // send X and theta into the solver
00855         newstate.X = newX;
00856         theta.resize(blockSize_);
00857         newstate.T = Teuchos::rcpFromRef(theta);
00858         // initialize
00859         lobpcg_solver->initialize(newstate);
00860       }
00861       catch (const 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       // don't catch any other exceptions
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       // copy them into the solution
00879       std::vector<int> which = ordertest->whichVecs();
00880       // indices between [0,blockSize) refer to vectors/values in the solver
00881       // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
00882       // everything has already been ordered by the solver; we just have to partition the two references
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           // sanity check
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       // set the vecs,vals in the solution
00899       if (insolver.size() > 0) {
00900         // set vecs
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         // set vals
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       // get the vecs,vals from locked storage
00914       if (inlocked.size() > 0) {
00915         int solnum = insolver.size();
00916         // set vecs
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         // set vals
00923         for (unsigned int i=0; i<inlocked.size(); i++) {
00924           vals[i+solnum] = lockvals[inlocked[i]];
00925         }
00926       }
00927 
00928       // sort the eigenvalues and permute the eigenvectors appropriately
00929       {
00930         std::vector<int> order(sol.numVecs);
00931         sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
00932         // store the values in the Eigensolution
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         // now permute the eigenvectors according to order
00938         msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
00939       }
00940 
00941       // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
00942       sol.index.resize(sol.numVecs,0);
00943     }
00944   }
00945 
00946   // print final summary
00947   lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00948 
00949   // print timing information
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   // get the number of iterations performed in this call to solve.
00956   numIters_ = lobpcg_solver->getNumIters();
00957 
00958   if (sol.numVecs < nev) {
00959     return Unconverged; // return from LOBPCGSolMgr::solve() 
00960   }
00961   return Converged; // return from LOBPCGSolMgr::solve() 
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 } // end Anasazi namespace
01011 
01012 #endif /* ANASAZI_LOBPCG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends