Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziRandomizedSolMgr.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ANASAZI_RANDOMIZED_SOLMGR_HPP
11#define ANASAZI_RANDOMIZED_SOLMGR_HPP
12
18#include "AnasaziConfigDefs.hpp"
19#include "AnasaziTypes.hpp"
20
21#include "AnasaziEigenproblem.hpp"
22#include "AnasaziSolverManager.hpp"
23
24#include "AnasaziBasicSort.hpp"
25#include "AnasaziSVQBOrthoManager.hpp"
26#include "AnasaziBasicOrthoManager.hpp"
27#include "AnasaziICGSOrthoManager.hpp"
28
29#include "AnasaziStatusTestMaxIters.hpp"
30#include "AnasaziStatusTestResNorm.hpp"
31#include "AnasaziStatusTestCombo.hpp"
32#include "AnasaziStatusTestOutput.hpp"
33#include "AnasaziOutputManager.hpp"
34#include "AnasaziOutputStreamTraits.hpp"
35#include "AnasaziSolverUtils.hpp"
36
37#include "Teuchos_TimeMonitor.hpp"
38#include "Teuchos_FancyOStream.hpp"
39#include "Teuchos_LAPACK.hpp"
40
56namespace Anasazi {
57
58 namespace Experimental {
59
60 template<class ScalarType, class MV, class OP>
61 class RandomizedSolMgr : public SolverManager<ScalarType,MV,OP> {
62
63 private:
64 typedef int OT;
65 typedef MultiVecTraits<ScalarType,MV> MVT;
66 typedef OperatorTraits<ScalarType,MV,OP> OPT;
67 typedef Teuchos::ScalarTraits<ScalarType> SCT;
68 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;
69 const ScalarType ONE = SCT::one();
70
71 public:
72
74
75
92 RandomizedSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
93 Teuchos::ParameterList &pl );
94
96 virtual ~RandomizedSolMgr() {};
98
100
101
102 const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
103 return *problem_;
104 }
105
106 int getNumIters() const {
107 return numIters_;
108 }
109
110 int getNumFailed() const {
111 return numFailed_;
112 }
113
115
117
118
127 ReturnType solve();
129
130 private:
131 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
132 Teuchos::RCP<Teuchos::FancyOStream> osp_;
133 std::string whch_;
134 MT tol_;
135 int osProc_;
136 int verb_;
137 Teuchos::RCP<Teuchos::Time> timerOrtho_;
138 Teuchos::RCP<Teuchos::Time> timerSolve_;
139 Teuchos::RCP<Teuchos::Time> timerOp_;
140 std::string ortho_;
141 int orthoFreq_;
142 int resFreq_;
143 int blockSize_;
144 int maxIters_;
145 int numIters_;
146 bool trackResNorms_;
147 int numFailed_;
148 };
149
150
152 template<class ScalarType, class MV, class OP>
153 RandomizedSolMgr<ScalarType,MV,OP>::RandomizedSolMgr(
154 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
155 Teuchos::ParameterList &pl ) :
156 problem_(problem),
157 whch_("LM"),
158 tol_(1e-6),
159 osProc_(0),
160 verb_(Anasazi::Errors),
161#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
162 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: Randomized::Orthogonalization")),
163 timerSolve_(Teuchos::TimeMonitor::getNewTimer("Anasazi: Randomized::solve()")),
164 timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: Randomized::Operation Op*x")),
165#endif
166 ortho_("SVQB"),
167 orthoFreq_(0),
168 resFreq_(0),
169 blockSize_(0),
170 maxIters_(5),
171 numIters_(0),
172 trackResNorms_(true)
173 {
174 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
175 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
176 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
177
178 whch_ = pl.get("Which",whch_);
179 TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "LM" && whch_ != "LR",
180 AnasaziError,
181 "RandomizedSolMgr: \"Which\" parameter must be LM or LR. Other options not available.");
182
183 tol_ = pl.get("Convergence Tolerance",tol_);
184 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
185 AnasaziError,
186 "RandomizedSolMgr: \"Tolerance\" parameter must be strictly positive.");
187
188 // Create a formatted output stream to print to.
189 // See if user requests output processor.
190 osProc_ = pl.get("Output Processor", osProc_);
191
192 // If not passed in by user, it will be chosen based upon operator type.
193 if (pl.isParameter("Output Stream")) {
194 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
195 }
196 else {
197 osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
198 }
199
200 // verbosity level
201 if (pl.isParameter("Verbosity")) {
202 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
203 verb_ = pl.get("Verbosity", verb_);
204 } else {
205 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
206 }
207 }
208
209 // Orthogonalization type
210 ortho_ = pl.get("Orthogonalization","SVQB");
211
212 blockSize_= pl.get("Block Size",problem_->getNEV());
213 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
214 AnasaziError,
215 "RandomizedSolMgr: \"Block Size\" parameter must be strictly positive.");
216
217 maxIters_ = pl.get("Maximum Iterations",maxIters_);
218 trackResNorms_ = pl.get("Track Residuals",true);
219
220 // How many iterations between orthogonalizations
221 if (pl.isParameter("Orthogonalization Frequency")) {
222 if (Teuchos::isParameterType<int>(pl,"Orthogonalization Frequency")) {
223 orthoFreq_ = pl.get("Orthogonalization Frequency", orthoFreq_);
224 } else {
225 orthoFreq_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Orthogonalization Frequency");
226 }
227 }
228
229 // How many iterations between checking the residuals
230 if (pl.isParameter("Residual Frequency")) {
231 if (Teuchos::isParameterType<int>(pl,"Residual Frequency")) {
232 resFreq_ = pl.get("Residual Frequency", resFreq_);
233 } else {
234 resFreq_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Residual Frequency");
235 }
236 }
237 }
238
240 template<class ScalarType, class MV, class OP>
241 ReturnType
242 RandomizedSolMgr<ScalarType,MV,OP>::solve() {
243
244 // Sort manager
245 Teuchos::RCP<BasicSort<MT> > sorter = Teuchos::rcp( new BasicSort<MT> );
246 sorter->setSortType(whch_);
247 std::vector<int> order(blockSize_); /* Permutation array for sorting the eigenvectors */
248 SolverUtils<ScalarType,MV,OP> msutils;
249
250 // Output manager
251 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp( new OutputManager<ScalarType>(verb_,osp_) );
252
253 // Eigensolution manager
254 Eigensolution<ScalarType,MV> sol;
255 sol.numVecs = 0;
256 problem_->setSolution(sol); /* In case there is an exception thrown */
257
258 // ortho manager
259 Teuchos::RCP<Anasazi::OrthoManager<ScalarType,MV> > orthoMgr;
260 int rank;
261 if (ortho_=="SVQB") {
262 orthoMgr = Teuchos::rcp( new Anasazi::SVQBOrthoManager<ScalarType,MV,OP>());
263 } else if (ortho_=="DGKS") {
264 orthoMgr = Teuchos::rcp( new Anasazi::BasicOrthoManager<ScalarType,MV,OP>());
265 } else if (ortho_=="ICGS") {
266 orthoMgr = Teuchos::rcp( new Anasazi::ICGSOrthoManager<ScalarType,MV,OP>());
267 } else {
268 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS"&&ortho_!="ICGS",std::logic_error,"Anasazi::RandomSolver Invalid orthogonalization type.");
269 }
270
271 if(blockSize_ < problem_->getNEV()){
272 printer->stream(Warnings) << "Warning! Block size smaller than number evals. Increasing Block Size to num evals." << std::endl;
273 blockSize_ = problem_->getNEV();
274 }
275
276 /* Grab some Multivector to Clone
277 * in practice, getInitVec() should always provide this, but it is possible to use a
278 * Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
279 * in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(),
280 * because we would like to clear the storage associated with V_ so we have room for the new V_ */
281 Teuchos::RCP<MV> randVecs;
282 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
283 "Anasazi::Randomized: eigenproblem did not specify initial vectors to clone from.");
284 if(MVT::GetNumberVecs(*(problem_->getInitVec()))==blockSize_){
285 randVecs = MVT::CloneCopy(*(problem_->getInitVec()));
286 } else {
287 randVecs = MVT::Clone(*(problem_->getInitVec()),blockSize_);
288 MVT::MvRandom(*randVecs);
289 }
290
291 { /* Ortho Timer */
292#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
293 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
294#endif
295 rank = orthoMgr->normalize(*randVecs);
296 if( rank < blockSize_ ) printer->stream(Warnings) << "Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
297 } /* End Ortho Timer */
298
299 /* Set up variables for residual computation ------------------------------ */
300 int i, j; // Loop variables
301
302 /* For computing H = Q^TAQ. (RR projection) */
303 Teuchos::RCP<MV> TmpVecs = MVT::Clone(*randVecs,blockSize_);
304 Teuchos::SerialDenseMatrix<OT,ScalarType> H (blockSize_, blockSize_);
305
306 /* For solving the projected eigenvalue problem. */
307 Teuchos::LAPACK<OT,ScalarType> lapack;
308 Teuchos::SerialDenseMatrix<OT,ScalarType> evects (blockSize_, blockSize_);
309 std::vector<MT> evals_real(blockSize_);
310 std::vector<MT> evals_imag(blockSize_);
311
312 /* Size of workspace and workspace for DGEEV */
313 int info = -1000;
314 ScalarType* vlr = 0;
315 const int ldv = 1;
316 int lwork = -1;
317 std::vector<ScalarType> work(1);
318 std::vector<MT> rwork(2*blockSize_);
319 numIters_ = 0;
320
321 /* For computing the residuals of the eigenproblem */
322 int numev;
323 std::vector<Value<ScalarType>> EigenVals(blockSize_);
324 Teuchos::RCP<MV> Avecs, evecs;
325 Teuchos::SerialDenseMatrix<OT,ScalarType> T (blockSize_, blockSize_ );
326 std::vector<MT> normV( blockSize_ );
327 bool converged = false;
328
329 { // Solve Timer
330#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
331 Teuchos::TimeMonitor slvtimer(*timerSolve_);
332#endif
333 sol.numVecs = blockSize_;
334 sol.Evals.resize(sol.numVecs);
335
336 // Perform multiplies by A and Rayleigh-Ritz
337 for( i = 0; i < maxIters_; i++ ){
338 if (converged == true) {
339 numFailed_ = 0;
340 numIters_ = i-1;
341 break;
342 }
343
344 { /* Begin Operator Timer */
345#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
346 Teuchos::TimeMonitor lcltimer( *timerOp_ );
347#endif
348 OPT::Apply( *(problem_->getOperator()), *randVecs, *randVecs );
349 } /* End Operator Timer */
350
351 if ((orthoFreq_ > 0 && i % orthoFreq_ == 0) || (resFreq_ > 0 && i % resFreq_ == 0)) {
352 { /* Start Ortho Timer */
353#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
354 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
355#endif
356 rank = orthoMgr->normalize(*randVecs);
357 if( rank < blockSize_ ) printer->stream(Warnings) << "Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
358 } /* End Ortho Timer */
359 } /* End Ortho */
360
361 if (resFreq_ > 0 && i % resFreq_ == 0) {
362
363 // Build the H matrix to run Rayleigh-Ritz on
364 OPT::Apply( *(problem_->getOperator()), *randVecs, *TmpVecs );
365 MVT::MvTransMv(ONE, *randVecs, *TmpVecs, H);
366
367 // Run GEEV once to find the correct size for rwork
368 lapack.GEEV('N','V',blockSize_,H.values(),H.stride(),evals_real.data(),evals_imag.data(),vlr, ldv, evects.values(), evects.stride(), &work[0], lwork, &rwork[0], &info);
369 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
370 work.resize( lwork );
371
372 // Run GEEV a second time to solve for Harmonic Ritz Values:
373 lapack.GEEV('N','V',blockSize_,H.values(),H.stride(),evals_real.data(),evals_imag.data(),vlr, ldv, evects.values(), evects.stride(), &work[0], lwork, &rwork[0], &info);
374 if(info != 0) printer->stream(Warnings) << "Warning! Anasazi::RandomSolver GEEV solve possible failure: info = " << info << std::endl;
375 lwork = -1;
376
377 // sort the eigenvalues and permute the eigenvectors appropriately
378 sorter->sort(evals_real,evals_imag,Teuchos::rcpFromRef(order),blockSize_);
379 msutils.permuteVectors(order, evects);
380
381 for( j = 0; j < blockSize_; j++){
382 EigenVals[j].realpart = evals_real[j];
383 EigenVals[j].imagpart = evals_imag[j];
384 }
385
386 // Project Evects back up to large problem.
387 MVT::MvTimesMatAddMv(ONE, *randVecs, evects, 0.0, *TmpVecs);
388
389 // Copy only the eigenvectors we asked for to soln.
390 sol.Evecs = MVT::CloneCopy(*TmpVecs, Teuchos::Range1D(0,sol.numVecs-1));
391 sol.numVecs = blockSize_;
392 sol.Evals = EigenVals;
393
394 // Extract evects/evals from solution
395 evecs = sol.Evecs;
396 numev = sol.numVecs;
397
398 // Check residuals for convergence
399 if (numev > 0 ) {
400 for ( j = 0; j < numev; ++j ) T(j, j) = sol.Evals[j].realpart;
401
402 Avecs = MVT::Clone(*evecs, numev);
403 OPT::Apply(*(problem_->getOperator()), *evecs, *Avecs);
404
405 MVT::MvTimesMatAddMv(-ONE, *evecs, T, ONE, *Avecs); /* Residuals = A*evecs - evecs*lambda */
406 MVT::MvNorm(*Avecs, normV);
407
408 numFailed_ = 0;
409 converged = true;
410 for ( j = 0; j < numev; ++j )
411 {
412 if ( SCT::magnitude(sol.Evals[j].realpart) != SCT::zero() ) normV[j] = SCT::magnitude(normV[j]/sol.Evals[j].realpart);
413 if (normV[j] > tol_) {
414 numFailed_++;
415 converged = false;
416 break;
417 }
418 }
419 } // End residual computation
420 } // End Rayleigh-Ritz solve
421 } // End subspace iterations
422
423 if(converged == false)
424 {
425 { /* Begin Ortho Timer */
426#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
427 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
428#endif
429 rank = orthoMgr->normalize(*randVecs);
430 if( rank < blockSize_ ) printer->stream(Warnings) << "Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
431 } /* End Ortho Timer */
432
433 OPT::Apply( *(problem_->getOperator()), *randVecs, *TmpVecs );
434 MVT::MvTransMv(ONE, *randVecs, *TmpVecs, H);
435
436 /* --------------------------------------------------------------------------
437 * Parameters for DGEEV:
438 * 'N' = Don't compute left eigenvectors.
439 * 'V' = Compute right eigenvectors.
440 * blockSize = Dimension of H (numEvals)
441 * H.values = H matrix (Q'AQ)
442 * H.stride = Leading dimension of H (numEvals)
443 * evals_real.data() = Array to store evals, real parts
444 * evals_imag.data() = Array to store evals, imag parts
445 * vlr = Stores left evects, so don't need this
446 * ldv = Leading dimension of vlr
447 * evects = Array to store right evects
448 * evects.stride = Leading dimension of evects
449 * work = Work array
450 * lwork = -1 means to query for array size
451 * rwork = Not referenced because ST is not complex
452 * -------------------------------------------------------------------------- */
453 //Find workspace size for DGEEV:
454 lapack.GEEV('N','V',blockSize_,H.values(),H.stride(),evals_real.data(),evals_imag.data(),vlr, ldv, evects.values(), evects.stride(), &work[0], lwork, &rwork[0], &info);
455 if(info != 0) printer->stream(IterationDetails) << "Warning!! Anasazi::RandomSolver GEEV solve possible failure: info = " << info << std::endl;
456
457 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
458 work.resize( lwork );
459
460 // Solve for Harmonic Ritz Values:
461 lapack.GEEV('N','V',blockSize_,H.values(),H.stride(),evals_real.data(),evals_imag.data(),vlr, ldv, evects.values(), evects.stride(), &work[0], lwork, &rwork[0], &info);
462 if(info != 0) printer->stream(IterationDetails) << "Warning!! Anasazi::RandomSolver GEEV solve possible failure: info = " << info << std::endl;
463
464 // Sort the eigenvalues and permute the eigenvectors appropriately
465 sorter->sort(evals_real,evals_imag,Teuchos::rcpFromRef(order),blockSize_);
466 msutils.permuteVectors(order, evects);
467
468 for( j = 0; j < blockSize_; j++){
469 EigenVals[j].realpart = evals_real[j];
470 EigenVals[j].imagpart = evals_imag[j];
471 }
472 sol.Evals = EigenVals;
473
474 // Project Evects back up to large problem and permute
475 MVT::MvTimesMatAddMv(ONE,*randVecs,evects, 0.0,*TmpVecs);
476
477 //------Post-Solve Processing----------------------------
478 //Copy only the eigenvectors we asked for to soln.
479 sol.numVecs = blockSize_;
480 sol.Evecs = MVT::CloneCopy(*TmpVecs, Teuchos::Range1D(0,sol.numVecs-1));
481 numIters_ = maxIters_;
482
483 } // End converged == false
484 } // End solve timer
485
486 // Send the solution to the eigenproblem
487 problem_->setSolution(sol);
488 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
489
490 // print timing information
491#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
492 if ( printer->isVerbosity( TimingDetails ) ) Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
493#endif
494
495 if (converged) return Converged; // Return from RandomizedSolMgr::solve()
496 return Unconverged; // Return from RandomizedSolMgr::solve()
497
498 } // End Solve function
499 } // end Experimental namespace
500} // end Anasazi namespace
501
502#endif /* ANASAZI_RANDOMIZED_SOLMGR_HPP */
AnasaziBasicOrthoManager.hpp
Basic implementation of the Anasazi::OrthoManager class.
AnasaziBasicSort.hpp
Basic implementation of the Anasazi::SortManager class.
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
AnasaziEigenproblem.hpp
Abstract base class which defines the interface required by an eigensolver and status test class to c...
AnasaziICGSOrthoManager.hpp
Basic implementation of the Anasazi::OrthoManager class.
AnasaziOutputManager.hpp
Abstract class definition for Anasazi Output Managers.
AnasaziOutputStreamTraits.hpp
Abstract class definition for Anasazi output stream.
AnasaziSVQBOrthoManager.hpp
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
AnasaziSolverManager.hpp
Pure virtual base class which describes the basic interface for a solver manager.
AnasaziSolverUtils.hpp
Class which provides internal utilities for the Anasazi solvers.
AnasaziStatusTestCombo.hpp
Status test for forming logical combinations of other status tests.
AnasaziStatusTestMaxIters.hpp
Status test for testing the number of iterations.
AnasaziStatusTestOutput.hpp
Special StatusTest for printing status tests.
AnasaziStatusTestResNorm.hpp
A status test for testing the norm of the eigenvectors residuals.
AnasaziTypes.hpp
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi::BasicOrthoManager
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Definition AnasaziBasicOrthoManager.hpp:39
Anasazi::ICGSOrthoManager
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Definition AnasaziICGSOrthoManager.hpp:39
RandomizedSolMgr
The Anasazi::RandomizedSolMgr approximates largest eigenvalues/eigenvectors by performing a simple Ra...
Anasazi::SVQBOrthoManager
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
Definition AnasaziSVQBOrthoManager.hpp:37
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Anasazi::ReturnType
ReturnType
Enumerated type used to pass back information from a solver manager.
Definition AnasaziTypes.hpp:89
Generated by doxygen 1.9.8