10#ifndef ANASAZI_RANDOMIZED_SOLMGR_HPP
11#define ANASAZI_RANDOMIZED_SOLMGR_HPP
37#include "Teuchos_TimeMonitor.hpp"
38#include "Teuchos_FancyOStream.hpp"
39#include "Teuchos_LAPACK.hpp"
60 template<
class ScalarType,
class MV,
class OP>
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();
92 RandomizedSolMgr(
const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
93 Teuchos::ParameterList &pl );
102 const Eigenproblem<ScalarType,MV,OP>& getProblem()
const {
106 int getNumIters()
const {
110 int getNumFailed()
const {
131 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
132 Teuchos::RCP<Teuchos::FancyOStream> osp_;
137 Teuchos::RCP<Teuchos::Time> timerOrtho_;
138 Teuchos::RCP<Teuchos::Time> timerSolve_;
139 Teuchos::RCP<Teuchos::Time> timerOp_;
152 template<
class ScalarType,
class MV,
class OP>
154 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
155 Teuchos::ParameterList &pl ) :
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")),
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.");
178 whch_ = pl.get(
"Which",whch_);
179 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"LM" && whch_ !=
"LR",
181 "RandomizedSolMgr: \"Which\" parameter must be LM or LR. Other options not available.");
183 tol_ = pl.get(
"Convergence Tolerance",tol_);
184 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
186 "RandomizedSolMgr: \"Tolerance\" parameter must be strictly positive.");
190 osProc_ = pl.get(
"Output Processor", osProc_);
193 if (pl.isParameter(
"Output Stream")) {
194 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
197 osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
201 if (pl.isParameter(
"Verbosity")) {
202 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
203 verb_ = pl.get(
"Verbosity", verb_);
205 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
210 ortho_ = pl.get(
"Orthogonalization",
"SVQB");
212 blockSize_= pl.get(
"Block Size",problem_->getNEV());
213 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
215 "RandomizedSolMgr: \"Block Size\" parameter must be strictly positive.");
217 maxIters_ = pl.get(
"Maximum Iterations",maxIters_);
218 trackResNorms_ = pl.get(
"Track Residuals",
true);
221 if (pl.isParameter(
"Orthogonalization Frequency")) {
222 if (Teuchos::isParameterType<int>(pl,
"Orthogonalization Frequency")) {
223 orthoFreq_ = pl.get(
"Orthogonalization Frequency", orthoFreq_);
225 orthoFreq_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Orthogonalization Frequency");
230 if (pl.isParameter(
"Residual Frequency")) {
231 if (Teuchos::isParameterType<int>(pl,
"Residual Frequency")) {
232 resFreq_ = pl.get(
"Residual Frequency", resFreq_);
234 resFreq_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Residual Frequency");
240 template<
class ScalarType,
class MV,
class OP>
245 Teuchos::RCP<BasicSort<MT> > sorter = Teuchos::rcp(
new BasicSort<MT> );
246 sorter->setSortType(whch_);
247 std::vector<int> order(blockSize_);
248 SolverUtils<ScalarType,MV,OP> msutils;
251 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp(
new OutputManager<ScalarType>(verb_,osp_) );
254 Eigensolution<ScalarType,MV> sol;
256 problem_->setSolution(sol);
259 Teuchos::RCP<Anasazi::OrthoManager<ScalarType,MV> > orthoMgr;
261 if (ortho_==
"SVQB") {
263 }
else if (ortho_==
"DGKS") {
265 }
else if (ortho_==
"ICGS") {
268 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS"&&ortho_!=
"ICGS",std::logic_error,
"Anasazi::RandomSolver Invalid orthogonalization type.");
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();
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()));
287 randVecs = MVT::Clone(*(problem_->getInitVec()),blockSize_);
288 MVT::MvRandom(*randVecs);
292#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
293 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
295 rank = orthoMgr->normalize(*randVecs);
296 if( rank < blockSize_ ) printer->stream(Warnings) <<
"Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
303 Teuchos::RCP<MV> TmpVecs = MVT::Clone(*randVecs,blockSize_);
304 Teuchos::SerialDenseMatrix<OT,ScalarType> H (blockSize_, blockSize_);
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_);
317 std::vector<ScalarType> work(1);
318 std::vector<MT> rwork(2*blockSize_);
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;
330#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
331 Teuchos::TimeMonitor slvtimer(*timerSolve_);
333 sol.numVecs = blockSize_;
334 sol.Evals.resize(sol.numVecs);
337 for( i = 0; i < maxIters_; i++ ){
338 if (converged ==
true) {
345#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
346 Teuchos::TimeMonitor lcltimer( *timerOp_ );
348 OPT::Apply( *(problem_->getOperator()), *randVecs, *randVecs );
351 if ((orthoFreq_ > 0 && i % orthoFreq_ == 0) || (resFreq_ > 0 && i % resFreq_ == 0)) {
353#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
354 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
356 rank = orthoMgr->normalize(*randVecs);
357 if( rank < blockSize_ ) printer->stream(Warnings) <<
"Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
361 if (resFreq_ > 0 && i % resFreq_ == 0) {
364 OPT::Apply( *(problem_->getOperator()), *randVecs, *TmpVecs );
365 MVT::MvTransMv(ONE, *randVecs, *TmpVecs, H);
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 );
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;
378 sorter->sort(evals_real,evals_imag,Teuchos::rcpFromRef(order),blockSize_);
379 msutils.permuteVectors(order, evects);
381 for( j = 0; j < blockSize_; j++){
382 EigenVals[j].realpart = evals_real[j];
383 EigenVals[j].imagpart = evals_imag[j];
387 MVT::MvTimesMatAddMv(ONE, *randVecs, evects, 0.0, *TmpVecs);
390 sol.Evecs = MVT::CloneCopy(*TmpVecs, Teuchos::Range1D(0,sol.numVecs-1));
391 sol.numVecs = blockSize_;
392 sol.Evals = EigenVals;
400 for ( j = 0; j < numev; ++j ) T(j, j) = sol.Evals[j].realpart;
402 Avecs = MVT::Clone(*evecs, numev);
403 OPT::Apply(*(problem_->getOperator()), *evecs, *Avecs);
405 MVT::MvTimesMatAddMv(-ONE, *evecs, T, ONE, *Avecs);
406 MVT::MvNorm(*Avecs, normV);
410 for ( j = 0; j < numev; ++j )
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_) {
423 if(converged ==
false)
426#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
427 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
429 rank = orthoMgr->normalize(*randVecs);
430 if( rank < blockSize_ ) printer->stream(Warnings) <<
"Warning! Anasazi::RandomSolver Random vectors did not have full rank!" << std::endl;
433 OPT::Apply( *(problem_->getOperator()), *randVecs, *TmpVecs );
434 MVT::MvTransMv(ONE, *randVecs, *TmpVecs, H);
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;
457 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
458 work.resize( lwork );
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;
465 sorter->sort(evals_real,evals_imag,Teuchos::rcpFromRef(order),blockSize_);
466 msutils.permuteVectors(order, evects);
468 for( j = 0; j < blockSize_; j++){
469 EigenVals[j].realpart = evals_real[j];
470 EigenVals[j].imagpart = evals_imag[j];
472 sol.Evals = EigenVals;
475 MVT::MvTimesMatAddMv(ONE,*randVecs,evects, 0.0,*TmpVecs);
479 sol.numVecs = blockSize_;
480 sol.Evecs = MVT::CloneCopy(*TmpVecs, Teuchos::Range1D(0,sol.numVecs-1));
481 numIters_ = maxIters_;
487 problem_->setSolution(sol);
488 printer->stream(Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
491#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
492 if ( printer->isVerbosity( TimingDetails ) ) Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
Basic implementation of the Anasazi::OrthoManager class.
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.
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 implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
The Anasazi::RandomizedSolMgr approximates largest eigenvalues/eigenvectors by performing a simple Ra...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
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.