Anasazi Version of the Day
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 
00123   int getNumIters() const {
00124     return numIters_;
00125   }
00126 
00128 
00130 
00131     
00140   ReturnType solve();
00142 
00143   private:
00144   Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00145   std::string whch_; 
00146   MagnitudeType tol_;
00147   int verb_;
00148   int blockSize_;
00149   int maxIters_;
00150   int numIters_;
00151 };
00152 
00153 
00155 template<class ScalarType, class MV, class OP>
00156 SimpleLOBPCGSolMgr<ScalarType,MV,OP>::SimpleLOBPCGSolMgr( 
00157         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00158         Teuchos::ParameterList &pl ) : 
00159   problem_(problem),
00160   whch_("LM"),
00161   tol_(1e-6),
00162   verb_(Anasazi::Errors),
00163   blockSize_(0),
00164   maxIters_(100),
00165   numIters_(0)
00166 {
00167   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00168   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00169   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00170   TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00171 
00172   whch_ = pl.get("Which","SR");
00173   TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
00174                      AnasaziError,
00175                      "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
00176 
00177   tol_ = pl.get("Convergence Tolerance",tol_);
00178   TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
00179                      AnasaziError,
00180                      "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
00181 
00182   // verbosity level
00183   if (pl.isParameter("Verbosity")) {
00184     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00185       verb_ = pl.get("Verbosity", verb_);
00186     } else {
00187       verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00188     }
00189   }
00190 
00191 
00192   blockSize_= pl.get("Block Size",problem_->getNEV());
00193   TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
00194                      AnasaziError,
00195                      "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
00196 
00197   maxIters_ = pl.get("Maximum Iterations",maxIters_);
00198 }
00199 
00200 
00201 
00203 template<class ScalarType, class MV, class OP>
00204 ReturnType 
00205 SimpleLOBPCGSolMgr<ScalarType,MV,OP>::solve() {
00206 
00207   // sort manager
00208   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00209   // output manager
00210   Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verb_) );
00211   // status tests
00212   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
00213   if (maxIters_ > 0) {
00214     max = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00215   }
00216   else {
00217     max = Teuchos::null;
00218   }
00219   Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm 
00220       = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(tol_) );
00221   Teuchos::Array< Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00222   alltests.push_back(norm);
00223   if (max != Teuchos::null) alltests.push_back(max);
00224   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combo 
00225       = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>(
00226               StatusTestCombo<ScalarType,MV,OP>::OR, alltests
00227         ));
00228   // printing StatusTest
00229   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
00230     = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
00231   // orthomanager
00232   Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho 
00233     = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00234   // parameter list
00235   Teuchos::ParameterList plist;
00236   plist.set("Block Size",blockSize_);
00237   plist.set("Full Ortho",true);
00238 
00239   // create an LOBPCG solver
00240   Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver 
00241     = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00242   // add the auxillary vecs from the eigenproblem to the solver
00243   if (problem_->getAuxVecs() != Teuchos::null) {
00244     lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
00245   }
00246 
00247   int numfound = 0;
00248   int nev = problem_->getNEV();
00249   Teuchos::Array< Teuchos::RCP<MV> > foundvecs;
00250   Teuchos::Array< Teuchos::RCP< std::vector<MagnitudeType> > > foundvals;
00251   while (numfound < nev) {
00252     // reduce the strain on norm test, if we are almost done
00253     if (nev - numfound < blockSize_) {
00254       norm->setQuorum(nev-numfound);
00255     }
00256 
00257     // tell the solver to iterate
00258     try {
00259       lobpcg_solver->iterate();
00260     }
00261     catch (const std::exception &e) {
00262       // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
00263       printer->stream(Anasazi::Errors) << "Exception: " << e.what() << std::endl;
00264       Eigensolution<ScalarType,MV> sol;
00265       sol.numVecs = 0;
00266       problem_->setSolution(sol);
00267       throw;
00268     }
00269 
00270     // check the status tests
00271     if (norm->getStatus() == Passed) {
00272 
00273       int num = norm->howMany();
00274       // if num < blockSize_, it is because we are on the last iteration: num+numfound>=nev
00275       TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
00276                          std::logic_error,
00277                          "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
00278       std::vector<int> ind = norm->whichVecs();
00279       // just grab the ones that we need
00280       if (num + numfound > nev) {
00281         num = nev - numfound;
00282         ind.resize(num);
00283       }
00284 
00285       // copy the converged eigenvectors
00286       Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00287       // store them
00288       foundvecs.push_back(newvecs);
00289       // add them as auxiliary vectors
00290       Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
00291       auxvecs.push_back(newvecs);
00292       // setAuxVecs() will reset the solver to uninitialized, without messing with numIters()
00293       lobpcg_solver->setAuxVecs(auxvecs);
00294 
00295       // copy the converged eigenvalues
00296       Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
00297       std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
00298       for (int i=0; i<num; i++) {
00299         (*newvals)[i] = all[ind[i]].realpart;
00300       }
00301       foundvals.push_back(newvals);
00302 
00303       numfound += num;
00304     }
00305     else if (max != Teuchos::null && max->getStatus() == Passed) {
00306 
00307       int num = norm->howMany();
00308       std::vector<int> ind = norm->whichVecs();
00309       
00310       if (num > 0) {
00311         // copy the converged eigenvectors
00312         Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00313         // store them
00314         foundvecs.push_back(newvecs);
00315         // don't bother adding them as auxiliary vectors; we have reached maxiters and are going to quit
00316         
00317         // copy the converged eigenvalues
00318         Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
00319         std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
00320         for (int i=0; i<num; i++) {
00321           (*newvals)[i] = all[ind[i]].realpart;
00322         }
00323         foundvals.push_back(newvals);
00324   
00325         numfound += num;
00326       }
00327       break;  // while(numfound < nev)
00328     }
00329     else {
00330       TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
00331     }
00332   } // end of while(numfound < nev)
00333 
00334   TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
00335 
00336   // create contiguous storage for all eigenvectors, eigenvalues
00337   Eigensolution<ScalarType,MV> sol;
00338   sol.numVecs = numfound;
00339   if (numfound > 0) {
00340     // allocate space for eigenvectors
00341     sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
00342   }
00343   else {
00344     sol.Evecs = Teuchos::null;
00345   }
00346   sol.Espace = sol.Evecs;
00347   // allocate space for eigenvalues
00348   std::vector<MagnitudeType> vals(numfound);
00349   sol.Evals.resize(numfound);
00350   // all real eigenvalues: set index vectors [0,...,numfound-1]
00351   sol.index.resize(numfound,0);
00352   // store eigenvectors, eigenvalues
00353   int curttl = 0;
00354   for (unsigned int i=0; i<foundvals.size(); i++) {
00355     TEUCHOS_TEST_FOR_EXCEPTION((signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
00356     unsigned int lclnum = foundvals[i]->size();
00357     std::vector<int> lclind(lclnum);
00358     for (unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
00359     // put the eigenvectors
00360     MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
00361     // put the eigenvalues
00362     std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
00363 
00364     curttl += lclnum;
00365   }
00366   TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
00367 
00368   // sort the eigenvalues and permute the eigenvectors appropriately
00369   if (numfound > 0) {
00370     std::vector<int> order(sol.numVecs);
00371     sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
00372     // store the values in the Eigensolution
00373     for (int i=0; i<sol.numVecs; i++) {
00374       sol.Evals[i].realpart = vals[i];
00375       sol.Evals[i].imagpart = MT::zero();
00376     }
00377     // now permute the eigenvectors according to order
00378     SolverUtils<ScalarType,MV,OP> msutils;
00379     msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
00380   }
00381 
00382   // print final summary
00383   lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00384 
00385   // print timing information
00386 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00387   if ( printer->isVerbosity( TimingDetails ) ) {
00388     Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
00389   }
00390 #endif
00391 
00392   // send the solution to the eigenproblem
00393   problem_->setSolution(sol);
00394   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00395 
00396   // get the number of iterations performed for this solve.
00397   numIters_ = lobpcg_solver->getNumIters();
00398  
00399   // return from SolMgr::solve()
00400   if (sol.numVecs < nev) return Unconverged;
00401   return Converged;
00402 }
00403 
00404 
00405 
00406 
00407 } // end Anasazi namespace
00408 
00409 #endif /* ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends