AnasaziSimpleLOBPCGSolMgr.hpp

Go to the documentation of this file.
00001 
00002 // @HEADER
00003 // ***********************************************************************
00004 //
00005 //                 Anasazi: Block Eigensolvers Package
00006 //                 Copyright (2004) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
00031 #define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
00032 
00038 #include "AnasaziConfigDefs.hpp"
00039 #include "AnasaziTypes.hpp"
00040 
00041 #include "AnasaziEigenproblem.hpp"
00042 #include "AnasaziSolverManager.hpp"
00043 
00044 #include "AnasaziLOBPCG.hpp"
00045 #include "AnasaziBasicSort.hpp"
00046 #include "AnasaziSVQBOrthoManager.hpp"
00047 #include "AnasaziStatusTestMaxIters.hpp"
00048 #include "AnasaziStatusTestResNorm.hpp"
00049 #include "AnasaziStatusTestCombo.hpp"
00050 #include "AnasaziStatusTestOutput.hpp"
00051 #include "AnasaziBasicOutputManager.hpp"
00052 #include "AnasaziSolverUtils.hpp"
00053 
00054 #include "Teuchos_TimeMonitor.hpp"
00055 
00083 namespace Anasazi {
00084 
00085 template<class ScalarType, class MV, class OP>
00086 class SimpleLOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00087 
00088   private:
00089     typedef MultiVecTraits<ScalarType,MV> MVT;
00090     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00091     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00092     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00093     
00094   public:
00095 
00097 
00098 
00109   SimpleLOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00110                              Teuchos::ParameterList &pl );
00111 
00113   virtual ~SimpleLOBPCGSolMgr() {};
00115   
00117 
00118 
00119   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00120     return *problem_;
00121   }
00122 
00124 
00126 
00127     
00136   ReturnType solve();
00138 
00139   private:
00140   Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00141   std::string whch_; 
00142   MagnitudeType tol_;
00143   int verb_;
00144   int blockSize_;
00145   int maxIters_;
00146 };
00147 
00148 
00150 template<class ScalarType, class MV, class OP>
00151 SimpleLOBPCGSolMgr<ScalarType,MV,OP>::SimpleLOBPCGSolMgr( 
00152         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00153         Teuchos::ParameterList &pl ) : 
00154   problem_(problem),
00155   whch_("LM"),
00156   tol_(1e-6),
00157   verb_(Anasazi::Errors),
00158   blockSize_(0),
00159   maxIters_(100)
00160 {
00161   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00162   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00163   TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00164   TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00165 
00166   whch_ = pl.get("Which","SR");
00167   TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
00168                      AnasaziError,
00169                      "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
00170 
00171   tol_ = pl.get("Convergence Tolerance",tol_);
00172   TEST_FOR_EXCEPTION(tol_ <= 0,
00173                      AnasaziError,
00174                      "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
00175 
00176   // verbosity level
00177   if (pl.isParameter("Verbosity")) {
00178     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00179       verb_ = pl.get("Verbosity", verb_);
00180     } else {
00181       verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00182     }
00183   }
00184 
00185 
00186   blockSize_= pl.get("Block Size",problem_->getNEV());
00187   TEST_FOR_EXCEPTION(blockSize_ <= 0,
00188                      AnasaziError,
00189                      "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
00190 
00191   maxIters_ = pl.get("Maximum Iterations",maxIters_);
00192 }
00193 
00194 
00195 
00197 template<class ScalarType, class MV, class OP>
00198 ReturnType 
00199 SimpleLOBPCGSolMgr<ScalarType,MV,OP>::solve() {
00200 
00201   // sort manager
00202   Teuchos::RCP<BasicSort<ScalarType,MV,OP> > sorter = Teuchos::rcp( new BasicSort<ScalarType,MV,OP>(whch_) );
00203   // output manager
00204   Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verb_) );
00205   // status tests
00206   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
00207   if (maxIters_ > 0) {
00208     max = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00209   }
00210   else {
00211     max = Teuchos::null;
00212   }
00213   Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm 
00214       = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(tol_) );
00215   Teuchos::Array< Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00216   alltests.push_back(norm);
00217   if (max != Teuchos::null) alltests.push_back(max);
00218   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combo 
00219       = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>(
00220               StatusTestCombo<ScalarType,MV,OP>::OR, alltests
00221         ));
00222   // printing StatusTest
00223   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
00224     = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
00225   // orthomanager
00226   Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho 
00227     = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00228   // parameter list
00229   Teuchos::ParameterList plist;
00230   plist.set("Block Size",blockSize_);
00231   plist.set("Full Ortho",true);
00232 
00233   // create an LOBPCG solver
00234   Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver 
00235     = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00236   // add the auxillary vecs from the eigenproblem to the solver
00237   if (problem_->getAuxVecs() != Teuchos::null) {
00238     lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
00239   }
00240 
00241   int numfound = 0;
00242   int nev = problem_->getNEV();
00243   Teuchos::Array< Teuchos::RCP<MV> > foundvecs;
00244   Teuchos::Array< Teuchos::RCP< std::vector<MagnitudeType> > > foundvals;
00245   while (numfound < nev) {
00246     // reduce the strain on norm test, if we are almost done
00247     if (nev - numfound < blockSize_) {
00248       norm->setQuorum(nev-numfound);
00249     }
00250 
00251     // tell the solver to iterate
00252     try {
00253       lobpcg_solver->iterate();
00254     }
00255     catch (std::exception e) {
00256       // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
00257       printer->stream(Anasazi::Errors) << "Exception: " << e.what() << std::endl;
00258       Eigensolution<ScalarType,MV> sol;
00259       sol.numVecs = 0;
00260       problem_->setSolution(sol);
00261       throw;
00262     }
00263 
00264     // check the status tests
00265     if (norm->getStatus() == Passed) {
00266 
00267       int num = norm->howMany();
00268       // if num < blockSize_, it is because we are on the last iteration: num+numfound>=nev
00269       TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
00270                          std::logic_error,
00271                          "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
00272       std::vector<int> ind = norm->whichVecs();
00273       // just grab the ones that we need
00274       if (num + numfound > nev) {
00275         num = nev - numfound;
00276         ind.resize(num);
00277       }
00278 
00279       // copy the converged eigenvectors
00280       Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00281       // store them
00282       foundvecs.push_back(newvecs);
00283       // add them as auxiliary vectors
00284       Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
00285       auxvecs.push_back(newvecs);
00286       // setAuxVecs() will reset the solver to uninitialized, without messing with numIters()
00287       lobpcg_solver->setAuxVecs(auxvecs);
00288 
00289       // copy the converged eigenvalues
00290       Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
00291       std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
00292       for (int i=0; i<num; i++) {
00293         (*newvals)[i] = all[ind[i]].realpart;
00294       }
00295       foundvals.push_back(newvals);
00296 
00297       numfound += num;
00298     }
00299     else if (max != Teuchos::null && max->getStatus() == Passed) {
00300 
00301       int num = norm->howMany();
00302       std::vector<int> ind = norm->whichVecs();
00303       
00304       if (num > 0) {
00305         // copy the converged eigenvectors
00306         Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00307         // orthornormalize to be safe
00308         ortho->normalizeMat(*newvecs,Teuchos::null,Teuchos::null);
00309         // store them
00310         foundvecs.push_back(newvecs);
00311         // don't bother adding them as auxiliary vectors; we have reached maxiters and are going to quit
00312         
00313         // copy the converged eigenvalues
00314         Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
00315         std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
00316         for (int i=0; i<num; i++) {
00317           (*newvals)[i] = all[ind[i]].realpart;
00318         }
00319         foundvals.push_back(newvals);
00320   
00321         numfound += num;
00322       }
00323       break;  // while(numfound < nev)
00324     }
00325     else {
00326       TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
00327     }
00328   } // end of while(numfound < nev)
00329 
00330   TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
00331 
00332   // create contiguous storage for all eigenvectors, eigenvalues
00333   Eigensolution<ScalarType,MV> sol;
00334   sol.numVecs = numfound;
00335   if (numfound > 0) {
00336     // allocate space for eigenvectors
00337     sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
00338   }
00339   else {
00340     sol.Evecs = Teuchos::null;
00341   }
00342   sol.Espace = sol.Evecs;
00343   // allocate space for eigenvalues
00344   std::vector<MagnitudeType> vals(numfound);
00345   sol.Evals.resize(numfound);
00346   // all real eigenvalues: set index vectors [0,...,numfound-1]
00347   sol.index.resize(numfound,0);
00348   // store eigenvectors, eigenvalues
00349   int curttl = 0;
00350   for (unsigned int i=0; i<foundvals.size(); i++) {
00351     TEST_FOR_EXCEPTION((signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
00352     unsigned int lclnum = foundvals[i]->size();
00353     std::vector<int> lclind(lclnum);
00354     for (unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
00355     // put the eigenvectors
00356     MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
00357     // put the eigenvalues
00358     copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
00359 
00360     curttl += lclnum;
00361   }
00362   TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
00363 
00364   // sort the eigenvalues and permute the eigenvectors appropriately
00365   if (numfound > 0) {
00366     std::vector<int> order(sol.numVecs);
00367     sorter->sort( lobpcg_solver.get(), sol.numVecs, vals, &order );
00368     // store the values in the Eigensolution
00369     for (int i=0; i<sol.numVecs; i++) {
00370       sol.Evals[i].realpart = vals[i];
00371       sol.Evals[i].imagpart = MT::zero();
00372     }
00373     // now permute the eigenvectors according to order
00374     SolverUtils<ScalarType,MV,OP> msutils;
00375     msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
00376   }
00377 
00378   // print final summary
00379   lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00380 
00381   // print timing information
00382   Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00383 
00384   // send the solution to the eigenproblem
00385   problem_->setSolution(sol);
00386   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00387 
00388   // return from SolMgr::solve()
00389   if (sol.numVecs < nev) return Unconverged;
00390   return Converged;
00391 }
00392 
00393 
00394 
00395 
00396 } // end Anasazi namespace
00397 
00398 #endif /* ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP */

Generated on Tue Oct 20 12:45:18 2009 for Anasazi by doxygen 1.4.7