00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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 TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00168 TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
00169 TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
00170 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 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 TEST_FOR_EXCEPTION(tol_ <= 0,
00179 AnasaziError,
00180 "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
00181
00182
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 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
00208 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00209
00210 Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verb_) );
00211
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
00229 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
00230 = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
00231
00232 Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
00233 = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00234
00235 Teuchos::ParameterList plist;
00236 plist.set("Block Size",blockSize_);
00237 plist.set("Full Ortho",true);
00238
00239
00240 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
00241 = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00242
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
00253 if (nev - numfound < blockSize_) {
00254 norm->setQuorum(nev-numfound);
00255 }
00256
00257
00258 try {
00259 lobpcg_solver->iterate();
00260 }
00261 catch (const std::exception &e) {
00262
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
00271 if (norm->getStatus() == Passed) {
00272
00273 int num = norm->howMany();
00274
00275 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
00280 if (num + numfound > nev) {
00281 num = nev - numfound;
00282 ind.resize(num);
00283 }
00284
00285
00286 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00287
00288 foundvecs.push_back(newvecs);
00289
00290 Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
00291 auxvecs.push_back(newvecs);
00292
00293 lobpcg_solver->setAuxVecs(auxvecs);
00294
00295
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
00312 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
00313
00314 foundvecs.push_back(newvecs);
00315
00316
00317
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;
00328 }
00329 else {
00330 TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
00331 }
00332 }
00333
00334 TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
00335
00336
00337 Eigensolution<ScalarType,MV> sol;
00338 sol.numVecs = numfound;
00339 if (numfound > 0) {
00340
00341 sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
00342 }
00343 else {
00344 sol.Evecs = Teuchos::null;
00345 }
00346 sol.Espace = sol.Evecs;
00347
00348 std::vector<MagnitudeType> vals(numfound);
00349 sol.Evals.resize(numfound);
00350
00351 sol.index.resize(numfound,0);
00352
00353 int curttl = 0;
00354 for (unsigned int i=0; i<foundvals.size(); i++) {
00355 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
00360 MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
00361
00362 std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
00363
00364 curttl += lclnum;
00365 }
00366 TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
00367
00368
00369 if (numfound > 0) {
00370 std::vector<int> order(sol.numVecs);
00371 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
00372
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
00378 SolverUtils<ScalarType,MV,OP> msutils;
00379 msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
00380 }
00381
00382
00383 lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00384
00385
00386 Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00387
00388
00389 problem_->setSolution(sol);
00390 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00391
00392
00393 numIters_ = lobpcg_solver->getNumIters();
00394
00395
00396 if (sol.numVecs < nev) return Unconverged;
00397 return Converged;
00398 }
00399
00400
00401
00402
00403 }
00404
00405 #endif