46 #ifndef ANASAZI_TRACEMIN_DAVIDSON_HPP 47 #define ANASAZI_TRACEMIN_DAVIDSON_HPP 56 #include "Teuchos_ScalarTraits.hpp" 57 #include "Teuchos_SerialDenseMatrix.hpp" 58 #include "Teuchos_ParameterList.hpp" 59 #include "Teuchos_TimeMonitor.hpp" 79 template <
class ScalarType,
class MV,
class OP>
92 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
96 Teuchos::ParameterList ¶ms
105 typedef Teuchos::ScalarTraits<ScalarType> SCT;
106 typedef typename SCT::magnitudeType MagnitudeType;
109 void addToBasis(
const Teuchos::RCP<const MV> Delta);
111 void harmonicAddToBasis(
const Teuchos::RCP<const MV> Delta);
126 template <
class ScalarType,
class MV,
class OP>
129 const Teuchos::RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
133 Teuchos::ParameterList ¶ms
135 TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params)
145 template <
class ScalarType,
class MV,
class OP>
149 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
150 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
154 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
156 Teuchos::RCP<MV> lclV, lclKV, lclMV;
158 Teuchos::Array< Teuchos::RCP<const MV> > projVecs = this->auxVecs_;
161 for(
int i=0; i<this->curDim_; i++)
163 lclV = MVT::CloneViewNonConst(*this->V_,curind);
164 projVecs.push_back(lclV);
167 for (
int i=0; i<this->blockSize_; ++i)
168 newind[i] = this->curDim_ + i;
169 lclV = MVT::CloneViewNonConst(*this->V_,newind);
172 MVT::SetBlock(*Delta,newind,*this->V_);
173 this->curDim_ += this->blockSize_;
179 Teuchos::Array< Teuchos::RCP<const MV> > MprojVecs = this->MauxVecs_;
180 lclMV = MVT::CloneViewNonConst(*this->MV_,curind);
181 MprojVecs.push_back(lclMV);
184 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
186 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 187 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
189 this->count_ApplyM_+= this->blockSize_;
190 OPT::Apply(*this->MOp_,*lclV,*lclMV);
194 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 195 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
199 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs,
200 Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
201 Teuchos::null,lclMV,MprojVecs);
204 MprojVecs.pop_back();
208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 209 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
213 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs);
218 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
219 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
222 if(this->Op_ != Teuchos::null)
224 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 225 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
227 this->count_ApplyOp_+= this->blockSize_;
229 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
230 OPT::Apply(*this->Op_,*lclV,*lclKV);
241 template <
class ScalarType,
class MV,
class OP>
242 void TraceMinDavidson<ScalarType,MV,OP>::harmonicAddToBasis(Teuchos::RCP<const MV> Delta)
245 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
246 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
250 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
252 Teuchos::RCP<MV> lclV, lclKV, lclMV, projVecs, KprojVecs;
256 for(
int i=0; i<this->curDim_; i++)
258 projVecs = MVT::CloneViewNonConst(*this->V_,curind);
260 if(this->Op_ != Teuchos::null)
262 lclKV = MVT::CloneViewNonConst(*this->KV_,curind);
267 for (
int i=0; i<this->blockSize_; ++i)
268 newind[i] = this->curDim_ + i;
269 lclV = MVT::CloneViewNonConst(*this->V_,newind);
272 MVT::SetBlock(*Delta,newind,*this->V_);
273 this->curDim_ += this->blockSize_;
276 if(this->auxVecs_.size() > 0)
277 this->orthman_->projectMat(*lclV,this->auxVecs_);
280 if(this->Op_ != Teuchos::null)
282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 283 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
285 this->count_ApplyOp_+= this->blockSize_;
287 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
288 OPT::Apply(*this->Op_,*lclV,*lclKV);
294 int nauxVecs = MVT::GetNumberVecs(*projVecs);
295 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(nauxVecs,this->blockSize_));
297 this->orthman_->innerProdMat(*KprojVecs,*lclKV,*gamma);
300 MVT::MvTimesMatAddMv(-this->ONE,*KprojVecs,*gamma,this->ONE,*lclKV);
303 MVT::MvTimesMatAddMv(-this->ONE,*projVecs,*gamma,this->ONE,*lclV);
306 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma2 = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
307 rank = this->orthman_->normalizeMat(*lclKV,gamma2);
310 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
311 SDsolver.setMatrix(gamma2);
313 RCP<MV> tempMV = MVT::CloneCopy(*lclV);
314 MVT::MvTimesMatAddMv(this->ONE,*tempMV,*gamma2,this->ZERO,*lclV);
316 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
317 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
322 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 323 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
325 this->count_ApplyM_+= this->blockSize_;
327 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
328 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.