181 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
182 Teuchos::ParameterList &pl ) :
192 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
193 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
194 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
195 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
197 whch_ = pl.get(
"Which",
"SR");
198 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR",
200 "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
202 tol_ = pl.get(
"Convergence Tolerance",tol_);
203 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
205 "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
209 osProc_ = pl.get(
"Output Processor", osProc_);
212 if (pl.isParameter(
"Output Stream")) {
213 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
220 if (pl.isParameter(
"Verbosity")) {
221 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
222 verb_ = pl.get(
"Verbosity", verb_);
224 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
229 blockSize_= pl.get(
"Block Size",problem_->getNEV());
230 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
232 "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
234 maxIters_ = pl.get(
"Maximum Iterations",maxIters_);
245 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp(
new BasicSort<MagnitudeType>(whch_) );
247 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp(
new OutputManager<ScalarType>(verb_,osp_) );
249 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
251 max = Teuchos::rcp(
new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
256 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm
257 = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(tol_) );
258 Teuchos::Array< Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
259 alltests.push_back(norm);
260 if (max != Teuchos::null) alltests.push_back(max);
261 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combo
262 = Teuchos::rcp(
new StatusTestCombo<ScalarType,MV,OP>(
266 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
267 = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,
Passed ) );
269 Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
270 = Teuchos::rcp(
new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
272 Teuchos::ParameterList plist;
273 plist.set(
"Block Size",blockSize_);
274 plist.set(
"Full Ortho",
true);
277 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
278 = Teuchos::rcp(
new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
280 if (problem_->getAuxVecs() != Teuchos::null) {
281 lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
285 int nev = problem_->getNEV();
286 Teuchos::Array< Teuchos::RCP<MV> > foundvecs;
287 Teuchos::Array< Teuchos::RCP< std::vector<MagnitudeType> > > foundvals;
288 while (numfound < nev) {
290 if (nev - numfound < blockSize_) {
291 norm->setQuorum(nev-numfound);
296 lobpcg_solver->iterate();
298 catch (
const std::exception &e) {
300 printer->stream(
Anasazi::Errors) <<
"Exception: " << e.what() << std::endl;
301 Eigensolution<ScalarType,MV> sol;
303 problem_->setSolution(sol);
308 if (norm->getStatus() ==
Passed) {
310 int num = norm->howMany();
312 TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
314 "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
315 std::vector<int> ind = norm->whichVecs();
317 if (num + numfound > nev) {
318 num = nev - numfound;
323 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
325 foundvecs.push_back(newvecs);
327 Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
328 auxvecs.push_back(newvecs);
330 lobpcg_solver->setAuxVecs(auxvecs);
333 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
334 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
335 for (
int i=0; i<num; i++) {
336 (*newvals)[i] = all[ind[i]].realpart;
338 foundvals.push_back(newvals);
342 else if (max != Teuchos::null && max->getStatus() ==
Passed) {
344 int num = norm->howMany();
345 std::vector<int> ind = norm->whichVecs();
349 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
351 foundvecs.push_back(newvecs);
355 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
356 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
357 for (
int i=0; i<num; i++) {
358 (*newvals)[i] = all[ind[i]].realpart;
360 foundvals.push_back(newvals);
367 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
371 TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
374 Eigensolution<ScalarType,MV> sol;
375 sol.numVecs = numfound;
378 sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
381 sol.Evecs = Teuchos::null;
383 sol.Espace = sol.Evecs;
385 std::vector<MagnitudeType> vals(numfound);
386 sol.Evals.resize(numfound);
388 sol.index.resize(numfound,0);
391 for (
unsigned int i=0; i<foundvals.size(); i++) {
392 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
393 unsigned int lclnum = foundvals[i]->size();
394 std::vector<int> lclind(lclnum);
395 for (
unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
397 MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
399 std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
403 TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
407 std::vector<int> order(sol.numVecs);
408 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
410 for (
int i=0; i<sol.numVecs; i++) {
411 sol.Evals[i].realpart = vals[i];
412 sol.Evals[i].imagpart = MT::zero();
415 SolverUtils<ScalarType,MV,OP> msutils;
416 msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
420 lobpcg_solver->currentStatus(printer->stream(
FinalSummary));
423#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
425 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
430 problem_->setSolution(sol);
431 printer->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
434 numIters_ = lobpcg_solver->getNumIters();