29 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP 30 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP 139 template<
class ScalarType,
class MV,
class OP>
194 std::vector<Value<ScalarType> > ret( _ritzValues );
205 return Teuchos::tuple(_timerSolve, _timerRestarting);
251 std::string _whch, _ortho;
252 MagnitudeType _ortho_kappa;
254 MagnitudeType _convtol;
256 bool _relconvtol,_conjSplit;
257 int _blockSize, _numBlocks, _stepSize, _nevBlocks, _xtra_nevBlocks;
260 bool _inSituRestart, _dynXtraNev;
262 std::vector<Value<ScalarType> > _ritzValues;
274 template<
class ScalarType,
class MV,
class OP>
293 _inSituRestart(false),
296 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
297 ,_timerSolve(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr::solve()")),
298 _timerRestarting(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr restarting"))
303 TEUCHOS_TEST_FOR_EXCEPTION(_problem->getInitVec() == Teuchos::null, std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
305 const int nev = _problem->getNEV();
308 _convtol = pl.
get(
"Convergence Tolerance",
MT::prec());
309 _relconvtol = pl.
get(
"Relative Convergence Tolerance",_relconvtol);
312 _maxRestarts = pl.
get(
"Maximum Restarts",_maxRestarts);
315 _blockSize = pl.
get(
"Block Size",1);
317 "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
320 _xtra_nevBlocks = pl.
get(
"Extra NEV Blocks",0);
321 if (nev%_blockSize) {
322 _nevBlocks = nev/_blockSize + 1;
324 _nevBlocks = nev/_blockSize;
331 if (Teuchos::isParameterType<bool>(pl,
"Dynamic Extra NEV")) {
332 _dynXtraNev = pl.
get(
"Dynamic Extra NEV",_dynXtraNev);
334 _dynXtraNev = ( Teuchos::getParameter<int>(pl,
"Dynamic Extra NEV") != 0 );
338 _numBlocks = pl.
get(
"Num Blocks",3*_nevBlocks);
340 "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
343 std::invalid_argument,
344 "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
348 _stepSize = pl.
get(
"Step Size", (_maxRestarts+1)*(_numBlocks+1));
350 _stepSize = pl.
get(
"Step Size", _numBlocks+1);
353 "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
357 _sort = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,
"Sort Manager");
360 _whch = pl.
get(
"Which",_whch);
361 TEUCHOS_TEST_FOR_EXCEPTION(_whch !=
"SM" && _whch !=
"LM" && _whch !=
"SR" && _whch !=
"LR" && _whch !=
"SI" && _whch !=
"LI",
362 std::invalid_argument,
"Invalid sorting string.");
367 _ortho = pl.
get(
"Orthogonalization",_ortho);
368 if (_ortho !=
"DGKS" && _ortho !=
"SVQB") {
373 _ortho_kappa = pl.
get(
"Orthogonalization Constant",_ortho_kappa);
377 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
378 _verbosity = pl.
get(
"Verbosity", _verbosity);
380 _verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
386 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
387 _inSituRestart = pl.
get(
"In Situ Restarting",_inSituRestart);
389 _inSituRestart = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
393 _printNum = pl.
get<
int>(
"Print Number of Ritz Values",-1);
398 template<
class ScalarType,
class MV,
class OP>
402 const int nev = _problem->getNEV();
419 if (globalTest_ == Teuchos::null) {
423 convtest = globalTest_;
431 if (debugTest_ != Teuchos::null) alltests.
push_back(debugTest_);
437 if ( printer->isVerbosity(
Debug) ) {
447 if (_ortho==
"SVQB") {
449 }
else if (_ortho==
"DGKS") {
450 if (_ortho_kappa <= 0) {
457 TEUCHOS_TEST_FOR_EXCEPTION(_ortho!=
"SVQB"&&_ortho!=
"DGKS",std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
463 plist.
set(
"Block Size",_blockSize);
464 plist.
set(
"Num Blocks",_numBlocks);
465 plist.
set(
"Step Size",_stepSize);
466 if (_printNum == -1) {
467 plist.
set(
"Print Number of Ritz Values",_nevBlocks*_blockSize);
470 plist.
set(
"Print Number of Ritz Values",_printNum);
479 if (probauxvecs != Teuchos::null) {
490 int maxXtraBlocks = 0;
491 if ( _dynXtraNev ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / _blockSize ) / 2;
494 if (_maxRestarts > 0) {
495 if (_inSituRestart==
true) {
497 workMV = MVT::Clone( *_problem->getInitVec(), 1 );
500 if (_problem->isHermitian()) {
501 workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + _blockSize );
504 workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + 1 + _blockSize );
508 workMV = Teuchos::null;
514 _problem->setSolution(sol);
517 int cur_nevBlocks = 0;
521 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 528 bks_solver->iterate();
535 if ( ordertest->getStatus() ==
Passed ) {
553 else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
554 (!_problem->isHermitian() && !_conjSplit && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
557 if (!bks_solver->isSchurCurrent()) {
558 bks_solver->computeSchurForm(
true );
561 outputtest->checkStatus( &*bks_solver );
565 if ( numRestarts >= _maxRestarts || ordertest->getStatus() ==
Passed) {
570 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 575 int numConv = ordertest->howMany();
576 cur_nevBlocks = _nevBlocks*_blockSize;
579 int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/_blockSize, _xtra_nevBlocks) );
581 cur_nevBlocks += moreNevBlocks * _blockSize;
582 else if ( _xtra_nevBlocks )
583 cur_nevBlocks += _xtra_nevBlocks * _blockSize;
591 printer->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << _maxRestarts << std::endl << std::endl;
592 printer->stream(
Debug) <<
" - Current NEV blocks is " << cur_nevBlocks <<
", the minimum is " << _nevBlocks*_blockSize << std::endl;
595 _ritzValues = bks_solver->getRitzValues();
601 int curDim = oldState.
curDim;
604 std::vector<int> ritzIndex = bks_solver->getRitzIndex();
605 if (ritzIndex[cur_nevBlocks-1]==1) {
614 if (_problem->isHermitian() && _conjSplit)
617 <<
" Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!" 619 <<
" Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!" 629 std::vector<int> curind( curDim );
630 for (
int i=0; i<curDim; i++) { curind[i] = i; }
642 if (_inSituRestart) {
649 std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
651 lapack.
GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
653 "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
655 std::vector<ScalarType> d(cur_nevBlocks);
656 for (
int j=0; j<copyQnev.numCols(); j++) {
657 d[j] = copyQnev(j,j);
659 if (printer->isVerbosity(
Debug)) {
661 for (
int j=0; j<R.
numCols(); j++) {
662 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
663 for (
int i=j+1; i<R.
numRows(); i++) {
667 printer->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.
normFrobenius() << std::endl;
673 curind.resize(curDim);
674 for (
int i=0; i<curDim; i++) curind[i] = i;
677 msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
681 curind.resize(cur_nevBlocks);
682 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
685 MVT::MvScale(*newV,d);
688 curind.resize(_blockSize);
689 for (
int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
690 newF = MVT::CloneViewNonConst( *solverbasis, curind );
694 curind.resize(cur_nevBlocks);
695 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
698 MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
699 tmp_newV = Teuchos::null;
701 curind.resize(_blockSize);
702 for (
int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
703 newF = MVT::CloneViewNonConst( *workMV, curind );
707 curind.resize(_blockSize);
708 for (
int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
710 for (
int i=0; i<_blockSize; i++) { curind[i] = i; }
711 MVT::SetBlock( *oldF, curind, *newF );
712 newF = Teuchos::null;
738 if (_inSituRestart) {
739 newstate.
V = oldState.
V;
744 newstate.
curDim = cur_nevBlocks;
745 bks_solver->initialize(newstate);
755 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
760 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
761 << err.what() << std::endl
762 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
769 workMV = Teuchos::null;
772 _ritzValues = bks_solver->getRitzValues();
774 sol.
numVecs = ordertest->howMany();
775 printer->stream(
Debug) <<
"ordertest->howMany() : " << sol.
numVecs << std::endl;
776 std::vector<int> whichVecs = ordertest->whichVecs();
782 std::vector<int> tmpIndex = bks_solver->getRitzIndex();
783 for (
int i=0; i<(int)_ritzValues.size(); ++i) {
784 printer->stream(
Debug) << _ritzValues[i].realpart <<
" + i " << _ritzValues[i].imagpart <<
", Index = " << tmpIndex[i] << std::endl;
786 printer->stream(
Debug) <<
"Number of converged eigenpairs (before) = " << sol.
numVecs << std::endl;
787 for (
int i=0; i<sol.
numVecs; ++i) {
788 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
790 if (tmpIndex[whichVecs[sol.
numVecs-1]]==1) {
791 printer->stream(
Debug) <<
"There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
792 whichVecs.push_back(whichVecs[sol.
numVecs-1]+1);
794 for (
int i=0; i<sol.
numVecs; ++i) {
795 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
799 bool keepMore =
false;
801 printer->stream(
Debug) <<
"Number of converged eigenpairs (after) = " << sol.
numVecs << std::endl;
802 printer->stream(
Debug) <<
"whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.
numVecs-1] <<
" > " << sol.
numVecs-1 << std::endl;
805 numEvecs = whichVecs[sol.
numVecs-1]+1;
806 printer->stream(
Debug) <<
"keepMore = true; numEvecs = " << numEvecs << std::endl;
810 bks_solver->setNumRitzVectors(numEvecs);
811 bks_solver->computeRitzVectors();
817 sol.
index = bks_solver->getRitzIndex();
818 sol.
Evals = bks_solver->getRitzValues();
819 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
829 std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
830 for (
int vec_i=0; vec_i<sol.
numVecs; ++vec_i) {
831 sol.
index[vec_i] = tmpIndex[whichVecs[vec_i]];
832 sol.
Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
834 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
843 bks_solver->currentStatus(printer->stream(
FinalSummary));
846 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 852 _problem->setSolution(sol);
853 printer->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
856 _numIters = bks_solver->getNumIters();
865 template <
class ScalarType,
class MV,
class OP>
870 globalTest_ = global;
873 template <
class ScalarType,
class MV,
class OP>
880 template <
class ScalarType,
class MV,
class OP>
888 template <
class ScalarType,
class MV,
class OP>
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Pure virtual base class which describes the basic interface for a solver manager. ...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
ScalarType * values() const
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Status test for forming logical combinations of other status tests.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
Structure to contain pointers to BlockKrylovSchur state variables.
virtual ~BlockKrylovSchurSolMgr()
Destructor.
An exception class parent to all Anasazi exceptions.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi's templated, static class providing utilities for the solvers.
void GEQRF(const OrdinalType m, const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType lwork, OrdinalType *info) const
int numVecs
The number of computed eigenpairs.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
BlockKrylovSchurSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockKrylovSchurSolMgr.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
OrdinalType numRows() const
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
Abstract base class which defines the interface required by an eigensolver and status test class to c...
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
static magnitudeType prec()
ReturnType
Enumerated type used to pass back information from a solver manager.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
void push_back(const value_type &x)
Teuchos::RCP< const MulVec > V
The current Krylov basis.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Status test for forming logical combinations of other status tests.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
int getNumIters() const
Get the iteration count for the most recent call to solve().
bool isParameter(const std::string &name) const
OrdinalType stride() const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
Implementation of a block Krylov-Schur eigensolver.
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems...
OrdinalType numCols() const
int curDim
The current dimension of the reduction.
Common interface of stopping criteria for Anasazi's solvers.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Class which provides internal utilities for the Anasazi solvers.