33 #ifndef ANASAZI_TRACEMIN_DAVIDSON_HPP 34 #define ANASAZI_TRACEMIN_DAVIDSON_HPP 43 #include "Teuchos_ScalarTraits.hpp" 44 #include "Teuchos_SerialDenseMatrix.hpp" 45 #include "Teuchos_ParameterList.hpp" 46 #include "Teuchos_TimeMonitor.hpp" 66 template <
class ScalarType,
class MV,
class OP>
79 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
83 Teuchos::ParameterList ¶ms
92 typedef Teuchos::ScalarTraits<ScalarType> SCT;
93 typedef typename SCT::magnitudeType MagnitudeType;
96 void addToBasis(
const Teuchos::RCP<const MV> Delta);
98 void harmonicAddToBasis(
const Teuchos::RCP<const MV> Delta);
113 template <
class ScalarType,
class MV,
class OP>
116 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
120 Teuchos::ParameterList ¶ms
122 TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params)
132 template <
class ScalarType,
class MV,
class OP>
136 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
137 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
141 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
143 Teuchos::RCP<MV> lclV, lclKV, lclMV;
145 Teuchos::Array< Teuchos::RCP<const MV> > projVecs = this->auxVecs_;
148 for(
int i=0; i<this->curDim_; i++)
150 lclV = MVT::CloneViewNonConst(*this->V_,curind);
151 projVecs.push_back(lclV);
154 for (
int i=0; i<this->blockSize_; ++i)
155 newind[i] = this->curDim_ + i;
156 lclV = MVT::CloneViewNonConst(*this->V_,newind);
159 MVT::SetBlock(*Delta,newind,*this->V_);
160 this->curDim_ += this->blockSize_;
166 Teuchos::Array< Teuchos::RCP<const MV> > MprojVecs = this->MauxVecs_;
167 lclMV = MVT::CloneViewNonConst(*this->MV_,curind);
168 MprojVecs.push_back(lclMV);
171 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
173 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 174 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
176 this->count_ApplyM_+= this->blockSize_;
177 OPT::Apply(*this->MOp_,*lclV,*lclMV);
181 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 182 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
186 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs,
187 Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
188 Teuchos::null,lclMV,MprojVecs);
191 MprojVecs.pop_back();
195 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 196 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
200 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs);
205 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
206 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
209 if(this->Op_ != Teuchos::null)
211 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 212 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
214 this->count_ApplyOp_+= this->blockSize_;
216 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
217 OPT::Apply(*this->Op_,*lclV,*lclKV);
228 template <
class ScalarType,
class MV,
class OP>
229 void TraceMinDavidson<ScalarType,MV,OP>::harmonicAddToBasis(Teuchos::RCP<const MV> Delta)
232 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
233 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
237 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
239 Teuchos::RCP<MV> lclV, lclKV, lclMV, projVecs, KprojVecs;
243 for(
int i=0; i<this->curDim_; i++)
245 projVecs = MVT::CloneViewNonConst(*this->V_,curind);
247 if(this->Op_ != Teuchos::null)
249 lclKV = MVT::CloneViewNonConst(*this->KV_,curind);
254 for (
int i=0; i<this->blockSize_; ++i)
255 newind[i] = this->curDim_ + i;
256 lclV = MVT::CloneViewNonConst(*this->V_,newind);
259 MVT::SetBlock(*Delta,newind,*this->V_);
260 this->curDim_ += this->blockSize_;
263 if(this->auxVecs_.size() > 0)
264 this->orthman_->projectMat(*lclV,this->auxVecs_);
267 if(this->Op_ != Teuchos::null)
269 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 270 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
272 this->count_ApplyOp_+= this->blockSize_;
274 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
275 OPT::Apply(*this->Op_,*lclV,*lclKV);
281 int nauxVecs = MVT::GetNumberVecs(*projVecs);
282 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(nauxVecs,this->blockSize_));
284 this->orthman_->innerProdMat(*KprojVecs,*lclKV,*gamma);
287 MVT::MvTimesMatAddMv(-this->ONE,*KprojVecs,*gamma,this->ONE,*lclKV);
290 MVT::MvTimesMatAddMv(-this->ONE,*projVecs,*gamma,this->ONE,*lclV);
293 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma2 = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
294 rank = this->orthman_->normalizeMat(*lclKV,gamma2);
297 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
298 SDsolver.setMatrix(gamma2);
300 RCP<MV> tempMV = MVT::CloneCopy(*lclV);
301 MVT::MvTimesMatAddMv(this->ONE,*tempMV,*gamma2,this->ZERO,*lclV);
303 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
304 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
309 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 310 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
312 this->count_ApplyM_+= this->blockSize_;
314 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
315 OPT::Apply(*this->MOp_,*lclV,*lclMV);
Abstract base class for trace minimization eigensolvers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
This class defines the interface required by an eigensolver and status test class to compute solution...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Output managers remove the need for the eigensolver to know any information about the required output...
TraceMinDavidson(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
TraceMinBase constructor with eigenproblem, solver utilities, and parameter list of solver options...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
This class implements a TraceMin-Davidson iteration for solving symmetric generalized eigenvalue prob...
This is an abstract base class for the trace minimization eigensolvers.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.