Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziSimpleLOBPCGSolMgr.hpp
Go to the documentation of this file.
1
2// @HEADER
3// ***********************************************************************
4//
5// Anasazi: Block Eigensolvers Package
6// Copyright 2004 Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41// @HEADER
42
43#ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
44#define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
45
51#include "AnasaziConfigDefs.hpp"
52#include "AnasaziTypes.hpp"
53
56
57#include "AnasaziLOBPCG.hpp"
58#include "AnasaziBasicSort.hpp"
67
68#include "Teuchos_TimeMonitor.hpp"
69#include "Teuchos_FancyOStream.hpp"
70
78
102namespace Anasazi {
103
104template<class ScalarType, class MV, class OP>
105class SimpleLOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
106
107 private:
109 typedef Teuchos::ScalarTraits<ScalarType> SCT;
110 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
111 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
112
113 public:
114
116
117
131 SimpleLOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
132 Teuchos::ParameterList &pl );
133
137
139
140
142 return *problem_;
143 }
144
145 int getNumIters() const {
146 return numIters_;
147 }
148
150
152
153
164
165 private:
166 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
167 Teuchos::RCP<Teuchos::FancyOStream> osp_;
168 std::string whch_;
169 MagnitudeType tol_;
170 int osProc_;
171 int verb_;
172 int blockSize_;
173 int maxIters_;
174 int numIters_;
175};
176
177
179template<class ScalarType, class MV, class OP>
181 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
182 Teuchos::ParameterList &pl ) :
183 problem_(problem),
184 whch_("LM"),
185 tol_(1e-6),
186 osProc_(0),
187 verb_(Anasazi::Errors),
188 blockSize_(0),
189 maxIters_(100),
190 numIters_(0)
191{
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.");
196
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.");
201
202 tol_ = pl.get("Convergence Tolerance",tol_);
203 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
205 "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
206
207 // Create a formatted output stream to print to.
208 // See if user requests output processor.
209 osProc_ = pl.get("Output Processor", osProc_);
210
211 // If not passed in by user, it will be chosen based upon operator type.
212 if (pl.isParameter("Output Stream")) {
213 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
214 }
215 else {
216 osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
217 }
218
219 // verbosity level
220 if (pl.isParameter("Verbosity")) {
221 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
222 verb_ = pl.get("Verbosity", verb_);
223 } else {
224 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
225 }
226 }
227
228
229 blockSize_= pl.get("Block Size",problem_->getNEV());
230 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
232 "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
233
234 maxIters_ = pl.get("Maximum Iterations",maxIters_);
235}
236
237
238
240template<class ScalarType, class MV, class OP>
243
244 // sort manager
245 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
246 // output manager
247 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp( new OutputManager<ScalarType>(verb_,osp_) );
248 // status tests
249 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
250 if (maxIters_ > 0) {
251 max = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
252 }
253 else {
254 max = Teuchos::null;
255 }
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>(
264 ));
265 // printing StatusTest
266 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
267 = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
268 // orthomanager
269 Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
270 = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
271 // parameter list
272 Teuchos::ParameterList plist;
273 plist.set("Block Size",blockSize_);
274 plist.set("Full Ortho",true);
275
276 // create an LOBPCG solver
277 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
278 = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
279 // add the auxillary vecs from the eigenproblem to the solver
280 if (problem_->getAuxVecs() != Teuchos::null) {
281 lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
282 }
283
284 int numfound = 0;
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) {
289 // reduce the strain on norm test, if we are almost done
290 if (nev - numfound < blockSize_) {
291 norm->setQuorum(nev-numfound);
292 }
293
294 // tell the solver to iterate
295 try {
296 lobpcg_solver->iterate();
297 }
298 catch (const std::exception &e) {
299 // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
300 printer->stream(Anasazi::Errors) << "Exception: " << e.what() << std::endl;
302 sol.numVecs = 0;
303 problem_->setSolution(sol);
304 throw;
305 }
306
307 // check the status tests
308 if (norm->getStatus() == Passed) {
309
310 int num = norm->howMany();
311 // if num < blockSize_, it is because we are on the last iteration: num+numfound>=nev
312 TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
313 std::logic_error,
314 "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
315 std::vector<int> ind = norm->whichVecs();
316 // just grab the ones that we need
317 if (num + numfound > nev) {
318 num = nev - numfound;
319 ind.resize(num);
320 }
321
322 // copy the converged eigenvectors
323 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
324 // store them
325 foundvecs.push_back(newvecs);
326 // add them as auxiliary vectors
327 Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
328 auxvecs.push_back(newvecs);
329 // setAuxVecs() will reset the solver to uninitialized, without messing with numIters()
330 lobpcg_solver->setAuxVecs(auxvecs);
331
332 // copy the converged eigenvalues
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;
337 }
338 foundvals.push_back(newvals);
339
340 numfound += num;
341 }
342 else if (max != Teuchos::null && max->getStatus() == Passed) {
343
344 int num = norm->howMany();
345 std::vector<int> ind = norm->whichVecs();
346
347 if (num > 0) {
348 // copy the converged eigenvectors
349 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
350 // store them
351 foundvecs.push_back(newvecs);
352 // don't bother adding them as auxiliary vectors; we have reached maxiters and are going to quit
353
354 // copy the converged eigenvalues
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;
359 }
360 foundvals.push_back(newvals);
361
362 numfound += num;
363 }
364 break; // while(numfound < nev)
365 }
366 else {
367 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
368 }
369 } // end of while(numfound < nev)
370
371 TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
372
373 // create contiguous storage for all eigenvectors, eigenvalues
375 sol.numVecs = numfound;
376 if (numfound > 0) {
377 // allocate space for eigenvectors
378 sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
379 }
380 else {
381 sol.Evecs = Teuchos::null;
382 }
383 sol.Espace = sol.Evecs;
384 // allocate space for eigenvalues
385 std::vector<MagnitudeType> vals(numfound);
386 sol.Evals.resize(numfound);
387 // all real eigenvalues: set index vectors [0,...,numfound-1]
388 sol.index.resize(numfound,0);
389 // store eigenvectors, eigenvalues
390 int curttl = 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;
396 // put the eigenvectors
397 MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
398 // put the eigenvalues
399 std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
400
401 curttl += lclnum;
402 }
403 TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
404
405 // sort the eigenvalues and permute the eigenvectors appropriately
406 if (numfound > 0) {
407 std::vector<int> order(sol.numVecs);
408 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
409 // store the values in the Eigensolution
410 for (int i=0; i<sol.numVecs; i++) {
411 sol.Evals[i].realpart = vals[i];
412 sol.Evals[i].imagpart = MT::zero();
413 }
414 // now permute the eigenvectors according to order
416 msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
417 }
418
419 // print final summary
420 lobpcg_solver->currentStatus(printer->stream(FinalSummary));
421
422 // print timing information
423#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
424 if ( printer->isVerbosity( TimingDetails ) ) {
425 Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
426 }
427#endif
428
429 // send the solution to the eigenproblem
430 problem_->setSolution(sol);
431 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
432
433 // get the number of iterations performed for this solve.
434 numIters_ = lobpcg_solver->getNumIters();
435
436 // return from SolMgr::solve()
437 if (sol.numVecs < nev) return Unconverged;
438 return Converged;
439}
440
441
442
443
444} // end Anasazi namespace
445
446#endif /* ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP */
Basic implementation of the Anasazi::SortManager class.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method.
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
Class which provides internal utilities for the Anasazi solvers.
Status test for forming logical combinations of other status tests.
Status test for testing the number of iterations.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
This class defines the interface required by an eigensolver and status test class to compute solution...
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration,...
Traits class which defines basic operations on multivectors.
Output managers remove the need for the eigensolver to know any information about the required output...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SimpleLOBPCGSolMgr provides a simple solver manager over the LOBPCG eigensolver.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
SimpleLOBPCGSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for SimpleLOBPCGSolMgr.
int getNumIters() const
Get the iteration count for the most recent call to solve().
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Anasazi's templated, static class providing utilities for the solvers.
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
Status test for forming logical combinations of other status tests.
A status test for testing the number of iterations.
A special StatusTest for printing other status tests.
A status test for testing the norm of the eigenvectors residuals.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Output managers remove the need for the eigensolver to know any information about the required output...