Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziRTRSolMgr.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef ANASAZI_RTR_SOLMGR_HPP
43#define ANASAZI_RTR_SOLMGR_HPP
44
50#include "AnasaziConfigDefs.hpp"
51#include "AnasaziTypes.hpp"
52
56
57#include "AnasaziIRTR.hpp"
58#include "AnasaziSIRTR.hpp"
59#include "AnasaziBasicSort.hpp"
68
69#include "Teuchos_TimeMonitor.hpp"
70#include "Teuchos_FancyOStream.hpp"
71
81namespace Anasazi {
82
83template<class ScalarType, class MV, class OP>
84class RTRSolMgr : public SolverManager<ScalarType,MV,OP> {
85
86 private:
87 typedef MultiVecTraits<ScalarType,MV> MVT;
88 typedef OperatorTraits<ScalarType,MV,OP> OPT;
89 typedef Teuchos::ScalarTraits<ScalarType> SCT;
90 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
91 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
92
93 public:
94
96
97
116 RTRSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
117 Teuchos::ParameterList &pl );
118
120 virtual ~RTRSolMgr() {};
122
124
125
127 const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
128 return *problem_;
129 }
130
136 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
137 return Teuchos::tuple(_timerSolve);
138 }
139
141 int getNumIters() const {
142 return numIters_;
143 }
144
145
147
149
150
161
162 private:
163 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
164 std::string whch_;
165 std::string ortho_;
166
167 bool skinny_;
168 MagnitudeType convtol_;
169 int maxIters_;
170 bool relconvtol_;
171 enum ResType convNorm_;
172 int numIters_;
173 int numICGS_;
174 int blkSize_;
175
176 Teuchos::RCP<Teuchos::Time> _timerSolve;
177 Teuchos::RCP<OutputManager<ScalarType> > printer_;
178 Teuchos::ParameterList& pl_;
179};
180
181
183template<class ScalarType, class MV, class OP>
185 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
186 Teuchos::ParameterList &pl ) :
187 problem_(problem),
188 whch_("SR"),
189 skinny_(true),
190 convtol_(MT::prec()),
191 maxIters_(100),
192 relconvtol_(true),
193 numIters_(-1),
194#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
195 _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: RTRSolMgr::solve()")),
196#endif
197 pl_(pl)
198{
199 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
200 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
201 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
202 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
203
204 std::string strtmp;
205
206 whch_ = pl_.get("Which","SR");
207 TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SR" && whch_ != "LR",
208 std::invalid_argument, "Anasazi::RTRSolMgr: Invalid sorting string. RTR solvers compute only LR or SR.");
209
210 // convergence tolerance
211 convtol_ = pl_.get("Convergence Tolerance",convtol_);
212 relconvtol_ = pl_.get("Relative Convergence Tolerance",relconvtol_);
213 strtmp = pl_.get("Convergence Norm",std::string("2"));
214 if (strtmp == "2") {
215 convNorm_ = RES_2NORM;
216 }
217 else if (strtmp == "M") {
218 convNorm_ = RES_ORTH;
219 }
220 else {
221 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
222 "Anasazi::RTRSolMgr: Invalid Convergence Norm.");
223 }
224
225
226 // maximum number of (outer) iterations
227 maxIters_ = pl_.get("Maximum Iterations",maxIters_);
228
229 // skinny solver or not
230 skinny_ = pl_.get("Skinny Solver",skinny_);
231
232 // number if ICGS iterations
233 numICGS_ = pl_.get("Num ICGS",2);
234
235 // block size
236 blkSize_ = pl_.get("Block Size", problem_->getNEV());
237
238 // Create a formatted output stream to print to.
239 // See if user requests output processor.
240 int osProc = pl.get("Output Processor", 0);
241
242 // If not passed in by user, it will be chosen based upon operator type.
243 Teuchos::RCP<Teuchos::FancyOStream> osp;
244
245 if (pl.isParameter("Output Stream")) {
246 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
247 }
248 else {
249 osp = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc);
250 }
251
252 int verbosity = Anasazi::Errors;
253 if (pl_.isParameter("Verbosity")) {
254 if (Teuchos::isParameterType<int>(pl_,"Verbosity")) {
255 verbosity = pl_.get("Verbosity", verbosity);
256 } else {
257 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl_,"Verbosity");
258 }
259 }
260 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity,osp) );
261
262}
263
264
265// solve()
266template<class ScalarType, class MV, class OP>
269
270 using std::endl;
271
272 // typedef SolverUtils<ScalarType,MV,OP> msutils; // unused
273
274 const int nev = problem_->getNEV();
275
276 // clear num iters
277 numIters_ = -1;
278
279#ifdef TEUCHOS_DEBUG
280 Teuchos::RCP<Teuchos::FancyOStream>
281 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
282 out->setShowAllFrontMatter(false).setShowProcRank(true);
283 *out << "Entering Anasazi::RTRSolMgr::solve()\n";
284#endif
285
287 // Sort manager
288 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
289
291 // Status tests
292 //
293 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > maxtest;
294 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > ordertest;
295 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > combotest;
296 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
297 // maximum iters
298 if (maxIters_ > 0) {
299 maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
300 }
301 else {
302 maxtest = Teuchos::null;
303 }
304 //
305 // convergence
306 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
307 ordertest = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
308 //
309 // combo
310 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
311 alltests.push_back(ordertest);
312 if (maxtest != Teuchos::null) alltests.push_back(maxtest);
313 combotest = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
314 //
315 // printing StatusTest
316 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
317 if ( printer_->isVerbosity(Debug) ) {
318 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
319 }
320 else {
321 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
322 }
323
325 // Orthomanager
326 Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho
327 = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM(),numICGS_) );
328
330 // create an RTR solver
331 // leftmost or rightmost?
332 bool leftMost = true;
333 if (whch_ == "LR" || whch_ == "LM") {
334 leftMost = false;
335 }
336 pl_.set<bool>("Leftmost",leftMost);
337 Teuchos::RCP<RTRBase<ScalarType,MV,OP> > rtr_solver;
338 if (skinny_ == false) {
339 // "hefty" IRTR
340 rtr_solver = Teuchos::rcp( new IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
341 }
342 else {
343 // "skinny" IRTR
344 rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
345 }
346 // set any auxiliary vectors defined in the problem
347 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
348 if (probauxvecs != Teuchos::null) {
349 rtr_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
350 }
351
352 TEUCHOS_TEST_FOR_EXCEPTION(rtr_solver->getBlockSize() < problem_->getNEV(),std::logic_error,
353 "Anasazi::RTRSolMgr requires block size >= requested number of eigenvalues.");
354
355 int numfound = 0;
356 Teuchos::RCP<MV> foundvecs;
357 std::vector<MagnitudeType> foundvals;
358
359 // tell the solver to iterate
360 try {
361#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
362 Teuchos::TimeMonitor slvtimer(*_timerSolve);
363#endif
364 rtr_solver->iterate();
365 numIters_ = rtr_solver->getNumIters();
366 }
367 catch (const std::exception &e) {
368 // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
369 printer_->stream(Anasazi::Errors) << "Exception: " << e.what() << endl;
370 Eigensolution<ScalarType,MV> sol;
371 sol.numVecs = 0;
372 problem_->setSolution(sol);
373 throw;
374 }
375
376 // check the status tests
377 if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed))
378 {
379 int num = convtest->howMany();
380 if (num > 0) {
381 std::vector<int> ind = convtest->whichVecs();
382 // copy the converged eigenvectors
383 foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
384 // copy the converged eigenvalues
385 foundvals.resize(num);
386 std::vector<Value<ScalarType> > all = rtr_solver->getRitzValues();
387 for (int i=0; i<num; i++) {
388 foundvals[i] = all[ind[i]].realpart;
389 }
390 numfound = num;
391 }
392 }
393 else {
394 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::RTRSolMgr::solve(): solver returned without satisfy status test.");
395 }
396
397 // create contiguous storage for all eigenvectors, eigenvalues
398 Eigensolution<ScalarType,MV> sol;
399 sol.numVecs = numfound;
400 sol.Evecs = foundvecs;
401 sol.Espace = sol.Evecs;
402 sol.Evals.resize(sol.numVecs);
403 for (int i=0; i<sol.numVecs; i++) {
404 sol.Evals[i].realpart = foundvals[i];
405 }
406 // all real eigenvalues: set index vectors [0,...,numfound-1]
407 sol.index.resize(numfound,0);
408
409 // print final summary
410 rtr_solver->currentStatus(printer_->stream(FinalSummary));
411
412 // print timing information
413#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
414 if ( printer_->isVerbosity( TimingDetails ) ) {
415 Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
416 }
417#endif
418
419 // send the solution to the eigenproblem
420 problem_->setSolution(sol);
421 printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
422
423 // return from SolMgr::solve()
424 if (sol.numVecs < nev) return Unconverged;
425 return Converged;
426}
427
428
429
430
431} // end Anasazi namespace
432
433#endif /* ANASAZI_RTR_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...
Basic implementation of the Anasazi::OrthoManager class.
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
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.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Types and exceptions used within Anasazi solvers and interfaces.
The Anasazi::RTRSolMgr provides a simple solver manager over the RTR eigensolver. For more informatio...
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
RTRSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for RTRSolMgr.
int getNumIters() const
Get the iteration count for the most recent call to solve.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
virtual ~RTRSolMgr()
Destructor.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Status test for forming logical combinations of other status tests.
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.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
Output managers remove the need for the eigensolver to know any information about the required output...