42 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
43 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
70 #include "Teuchos_FancyOStream.hpp"
155 template<
class ScalarType,
class MV,
class OP>
213 std::vector<Value<ScalarType> > ret( ritzValues_ );
224 return Teuchos::tuple(timerSolve_, timerRestarting_);
270 std::string whch_, ortho_;
271 MagnitudeType ortho_kappa_;
273 MagnitudeType convtol_;
275 bool relconvtol_,conjSplit_;
276 int blockSize_, numBlocks_, stepSize_, nevBlocks_, xtra_nevBlocks_;
279 bool inSituRestart_, dynXtraNev_;
282 std::vector<Value<ScalarType> > ritzValues_;
296 template<
class ScalarType,
class MV,
class OP>
314 verbosity_(Anasazi::
Errors),
315 inSituRestart_(false),
319 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
320 ,timerSolve_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr::solve()")),
321 timerRestarting_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: BlockKrylovSchurSolMgr restarting"))
326 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
328 const int nev = problem_->getNEV();
331 convtol_ = pl.
get(
"Convergence Tolerance",
MT::prec());
332 relconvtol_ = pl.
get(
"Relative Convergence Tolerance",relconvtol_);
335 maxRestarts_ = pl.
get(
"Maximum Restarts",maxRestarts_);
338 blockSize_ = pl.
get(
"Block Size",1);
340 "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
343 xtra_nevBlocks_ = pl.
get(
"Extra NEV Blocks",0);
344 if (nev%blockSize_) {
345 nevBlocks_ = nev/blockSize_ + 1;
347 nevBlocks_ = nev/blockSize_;
354 if (Teuchos::isParameterType<bool>(pl,
"Dynamic Extra NEV")) {
355 dynXtraNev_ = pl.
get(
"Dynamic Extra NEV",dynXtraNev_);
357 dynXtraNev_ = ( Teuchos::getParameter<int>(pl,
"Dynamic Extra NEV") != 0 );
361 numBlocks_ = pl.
get(
"Num Blocks",3*nevBlocks_);
363 "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
366 std::invalid_argument,
367 "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
371 stepSize_ = pl.
get(
"Step Size", (maxRestarts_+1)*(numBlocks_+1));
373 stepSize_ = pl.
get(
"Step Size", numBlocks_+1);
376 "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
380 sort_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,
"Sort Manager");
383 whch_ = pl.
get(
"Which",whch_);
384 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR" && whch_ !=
"SI" && whch_ !=
"LI",
385 std::invalid_argument,
"Invalid sorting string.");
390 ortho_ = pl.
get(
"Orthogonalization",ortho_);
391 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB") {
396 ortho_kappa_ = pl.
get(
"Orthogonalization Constant",ortho_kappa_);
400 osProc_ = pl.
get(
"Output Processor", osProc_);
404 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
412 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
413 verbosity_ = pl.
get(
"Verbosity", verbosity_);
415 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
421 if (Teuchos::isParameterType<bool>(pl,
"In Situ Restarting")) {
422 inSituRestart_ = pl.
get(
"In Situ Restarting",inSituRestart_);
424 inSituRestart_ = ( Teuchos::getParameter<int>(pl,
"In Situ Restarting") != 0 );
428 printNum_ = pl.
get<
int>(
"Print Number of Ritz Values",printNum_);
433 template<
class ScalarType,
class MV,
class OP>
437 const int nev = problem_->getNEV();
454 if (globalTest_ == Teuchos::null) {
458 convtest = globalTest_;
466 if (debugTest_ != Teuchos::null) alltests.
push_back(debugTest_);
472 if ( printer->isVerbosity(
Debug) ) {
482 if (ortho_==
"SVQB") {
484 }
else if (ortho_==
"DGKS") {
485 if (ortho_kappa_ <= 0) {
492 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS",std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
498 plist.
set(
"Block Size",blockSize_);
499 plist.
set(
"Num Blocks",numBlocks_);
500 plist.
set(
"Step Size",stepSize_);
501 if (printNum_ == -1) {
502 plist.
set(
"Print Number of Ritz Values",nevBlocks_*blockSize_);
505 plist.
set(
"Print Number of Ritz Values",printNum_);
514 if (probauxvecs != Teuchos::null) {
525 int maxXtraBlocks = 0;
526 if ( dynXtraNev_ ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / blockSize_ ) / 2;
529 if (maxRestarts_ > 0) {
530 if (inSituRestart_==
true) {
532 workMV = MVT::Clone( *problem_->getInitVec(), 1 );
535 if (problem_->isHermitian()) {
536 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + blockSize_ );
539 workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + 1 + blockSize_ );
543 workMV = Teuchos::null;
549 problem_->setSolution(sol);
552 int cur_nevBlocks = 0;
556 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
563 bks_solver->iterate();
570 if ( ordertest->getStatus() ==
Passed ) {
588 else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
589 (!problem_->isHermitian() && !conjSplit_ && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
592 if (!bks_solver->isSchurCurrent()) {
593 bks_solver->computeSchurForm(
true );
596 outputtest->checkStatus( &*bks_solver );
600 if ( numRestarts >= maxRestarts_ || ordertest->getStatus() ==
Passed) {
605 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
610 int numConv = ordertest->howMany();
611 cur_nevBlocks = nevBlocks_*blockSize_;
614 int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/blockSize_, xtra_nevBlocks_) );
616 cur_nevBlocks += moreNevBlocks * blockSize_;
617 else if ( xtra_nevBlocks_ )
618 cur_nevBlocks += xtra_nevBlocks_ * blockSize_;
626 printer->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
627 printer->stream(
Debug) <<
" - Current NEV blocks is " << cur_nevBlocks <<
", the minimum is " << nevBlocks_*blockSize_ << std::endl;
630 ritzValues_ = bks_solver->getRitzValues();
636 int curDim = oldState.
curDim;
639 std::vector<int> ritzIndex = bks_solver->getRitzIndex();
640 if (ritzIndex[cur_nevBlocks-1]==1) {
649 if (problem_->isHermitian() && conjSplit_)
652 <<
" Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
654 <<
" Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
664 std::vector<int> curind( curDim );
665 for (
int i=0; i<curDim; i++) { curind[i] = i; }
677 if (inSituRestart_) {
684 std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
686 lapack.
GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
688 "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
690 std::vector<ScalarType> d(cur_nevBlocks);
691 for (
int j=0; j<copyQnev.numCols(); j++) {
692 d[j] = copyQnev(j,j);
694 if (printer->isVerbosity(
Debug)) {
696 for (
int j=0; j<R.
numCols(); j++) {
697 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
698 for (
int i=j+1; i<R.
numRows(); i++) {
702 printer->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.
normFrobenius() << std::endl;
708 curind.resize(curDim);
709 for (
int i=0; i<curDim; i++) curind[i] = i;
712 msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
716 curind.resize(cur_nevBlocks);
717 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
720 MVT::MvScale(*newV,d);
723 curind.resize(blockSize_);
724 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
725 newF = MVT::CloneViewNonConst( *solverbasis, curind );
729 curind.resize(cur_nevBlocks);
730 for (
int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
733 MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
734 tmp_newV = Teuchos::null;
736 curind.resize(blockSize_);
737 for (
int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
738 newF = MVT::CloneViewNonConst( *workMV, curind );
742 curind.resize(blockSize_);
743 for (
int i=0; i<blockSize_; i++) { curind[i] = curDim + i; }
745 for (
int i=0; i<blockSize_; i++) { curind[i] = i; }
746 MVT::SetBlock( *oldF, curind, *newF );
747 newF = Teuchos::null;
773 if (inSituRestart_) {
774 newstate.
V = oldState.
V;
779 newstate.
curDim = cur_nevBlocks;
780 bks_solver->initialize(newstate);
790 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
795 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
796 << err.what() << std::endl
797 <<
"Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
804 workMV = Teuchos::null;
807 ritzValues_ = bks_solver->getRitzValues();
809 sol.
numVecs = ordertest->howMany();
810 printer->stream(
Debug) <<
"ordertest->howMany() : " << sol.
numVecs << std::endl;
811 std::vector<int> whichVecs = ordertest->whichVecs();
817 std::vector<int> tmpIndex = bks_solver->getRitzIndex();
818 for (
int i=0; i<(int)ritzValues_.size(); ++i) {
819 printer->stream(
Debug) << ritzValues_[i].realpart <<
" + i " << ritzValues_[i].imagpart <<
", Index = " << tmpIndex[i] << std::endl;
821 printer->stream(
Debug) <<
"Number of converged eigenpairs (before) = " << sol.
numVecs << std::endl;
822 for (
int i=0; i<sol.
numVecs; ++i) {
823 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
825 if (tmpIndex[whichVecs[sol.
numVecs-1]]==1) {
826 printer->stream(
Debug) <<
"There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
827 whichVecs.push_back(whichVecs[sol.
numVecs-1]+1);
829 for (
int i=0; i<sol.
numVecs; ++i) {
830 printer->stream(
Debug) <<
"whichVecs[" << i <<
"] = " << whichVecs[i] <<
", tmpIndex[" << whichVecs[i] <<
"] = " << tmpIndex[whichVecs[i]] << std::endl;
834 bool keepMore =
false;
836 printer->stream(
Debug) <<
"Number of converged eigenpairs (after) = " << sol.
numVecs << std::endl;
837 printer->stream(
Debug) <<
"whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.
numVecs-1] <<
" > " << sol.
numVecs-1 << std::endl;
840 numEvecs = whichVecs[sol.
numVecs-1]+1;
841 printer->stream(
Debug) <<
"keepMore = true; numEvecs = " << numEvecs << std::endl;
845 bks_solver->setNumRitzVectors(numEvecs);
846 bks_solver->computeRitzVectors();
852 sol.
index = bks_solver->getRitzIndex();
853 sol.
Evals = bks_solver->getRitzValues();
854 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
864 std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
865 for (
int vec_i=0; vec_i<sol.
numVecs; ++vec_i) {
866 sol.
index[vec_i] = tmpIndex[whichVecs[vec_i]];
867 sol.
Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
869 sol.
Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
878 bks_solver->currentStatus(printer->stream(
FinalSummary));
881 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
887 problem_->setSolution(sol);
888 printer->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
891 numIters_ = bks_solver->getNumIters();
900 template <
class ScalarType,
class MV,
class OP>
905 globalTest_ = global;
908 template <
class ScalarType,
class MV,
class OP>
915 template <
class ScalarType,
class MV,
class OP>
923 template <
class ScalarType,
class MV,
class OP>
ScalarType * values() const
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. ...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
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 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.
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
An exception class parent to all Anasazi exceptions.
Output managers remove the need for the eigensolver to know any information about the required output...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
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.
bool isParameter(const std::string &name) const
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
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.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
Output managers remove the need for the eigensolver to know any information about the required output...
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...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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 Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
OrdinalType numCols() const
void push_back(const value_type &x)
void GEQRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
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.
Abstract class definition for Anasazi output stream.
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...
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.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
OrdinalType stride() const
int getNumIters() const
Get the iteration count for the most recent call to solve().
OrdinalType numRows() const
Class which provides internal utilities for the Anasazi solvers.