Anasazi  Version of the Day
AnasaziTraceMinDavidson.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_TRACEMIN_DAVIDSON_HPP
34 #define ANASAZI_TRACEMIN_DAVIDSON_HPP
35 
36 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziEigensolver.hpp"
41 #include "AnasaziTraceMinBase.hpp"
42 
43 #include "Teuchos_ScalarTraits.hpp"
44 #include "Teuchos_SerialDenseMatrix.hpp"
45 #include "Teuchos_ParameterList.hpp"
46 #include "Teuchos_TimeMonitor.hpp"
47 
48 
49 namespace Anasazi {
50 namespace Experimental {
51 
66  template <class ScalarType, class MV, class OP>
67  class TraceMinDavidson : public TraceMinBase<ScalarType,MV,OP> {
68  public:
69 
78  TraceMinDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
79  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
80  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
81  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
82  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
83  Teuchos::ParameterList &params
84  );
85 
86  private:
87  //
88  // Convenience typedefs
89  //
92  typedef Teuchos::ScalarTraits<ScalarType> SCT;
93  typedef typename SCT::magnitudeType MagnitudeType;
94 
95  // TraceMin specific methods
96  void addToBasis(const Teuchos::RCP<const MV> Delta);
97 
98  void harmonicAddToBasis(const Teuchos::RCP<const MV> Delta);
99  };
100 
103  //
104  // Implementations
105  //
108 
109 
110 
112  // Constructor
113  template <class ScalarType, class MV, class OP>
115  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
116  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
117  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
118  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
119  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
120  Teuchos::ParameterList &params
121  ) :
122  TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params)
123  {
124  }
125 
126 
128  // 1. Project Delta so that V' M Delta = 0 and Q' M Delta = 0
129  // 2. Normalize Delta so that Delta' M Delta = I
130  // 3. Add Delta to the end of V: V = [V Delta]
131  // 4. Update KV and MV
132  template <class ScalarType, class MV, class OP>
133  void TraceMinDavidson<ScalarType,MV,OP>::addToBasis(Teuchos::RCP<const MV> Delta)
134  {
135  // TODO: We should also test the row length and map, etc...
136  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
137  "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
138 
139  int rank;
140  // Vector of indices
141  std::vector<int> curind(this->curDim_), newind(this->blockSize_);
142  // Pointer to the meaningful parts of V, KV, and MV
143  Teuchos::RCP<MV> lclV, lclKV, lclMV;
144  // Holds the vectors we project against
145  Teuchos::Array< Teuchos::RCP<const MV> > projVecs = this->auxVecs_;
146 
147  // Get the existing parts of the basis and add them to the list of things we project against
148  for(int i=0; i<this->curDim_; i++)
149  curind[i] = i;
150  lclV = MVT::CloneViewNonConst(*this->V_,curind);
151  projVecs.push_back(lclV);
152 
153  // Get the new part of the basis (where we're going to insert Delta)
154  for (int i=0; i<this->blockSize_; ++i)
155  newind[i] = this->curDim_ + i;
156  lclV = MVT::CloneViewNonConst(*this->V_,newind);
157 
158  // Insert Delta at the end of V
159  MVT::SetBlock(*Delta,newind,*this->V_);
160  this->curDim_ += this->blockSize_;
161 
162  // Project out the components of Delta in the direction of V
163  if(this->hasM_)
164  {
165  // It is more efficient to provide the orthomanager with MV
166  Teuchos::Array< Teuchos::RCP<const MV> > MprojVecs = this->MauxVecs_;
167  lclMV = MVT::CloneViewNonConst(*this->MV_,curind);
168  MprojVecs.push_back(lclMV);
169 
170  // Compute M*Delta
171  lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
172  {
173  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
174  Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
175  #endif
176  this->count_ApplyM_+= this->blockSize_;
177  OPT::Apply(*this->MOp_,*lclV,*lclMV);
178  }
179 
180  {
181  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
182  Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
183  #endif
184 
185  // Project and normalize Delta
186  rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs,
187  Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
188  Teuchos::null,lclMV,MprojVecs);
189  }
190 
191  MprojVecs.pop_back();
192  }
193  else
194  {
195  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
196  Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
197  #endif
198 
199  // Project and normalize Delta
200  rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs);
201  }
202 
203  projVecs.pop_back();
204 
205  TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
206  "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
207 
208  // Update KV
209  if(this->Op_ != Teuchos::null)
210  {
211  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
212  Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
213  #endif
214  this->count_ApplyOp_+= this->blockSize_;
215 
216  lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
217  OPT::Apply(*this->Op_,*lclV,*lclKV);
218  }
219  }
220 
221 
222 
224  // 1. Project Delta so that V' M Delta = 0 and Q' M Delta = 0
225  // 2. Normalize Delta so that Delta' M Delta = I
226  // 3. Add Delta to the end of V: V = [V Delta]
227  // 4. Update KV and MV
228  template <class ScalarType, class MV, class OP>
229  void TraceMinDavidson<ScalarType,MV,OP>::harmonicAddToBasis(Teuchos::RCP<const MV> Delta)
230  {
231  // TODO: We should also test the row length and map, etc...
232  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
233  "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
234 
235  int rank;
236  // Vector of indices
237  std::vector<int> curind(this->curDim_), newind(this->blockSize_);
238  // Pointer to the meaningful parts of V, KV, and MV
239  Teuchos::RCP<MV> lclV, lclKV, lclMV, projVecs, KprojVecs;
240  // Array of things we project against
241 
242  // Get the existing parts of the basis and add them to the list of things we project against
243  for(int i=0; i<this->curDim_; i++)
244  curind[i] = i;
245  projVecs = MVT::CloneViewNonConst(*this->V_,curind);
246 
247  if(this->Op_ != Teuchos::null)
248  {
249  lclKV = MVT::CloneViewNonConst(*this->KV_,curind);
250  KprojVecs = lclKV;
251  }
252 
253  // Get the new part of the basis (where we're going to insert Delta)
254  for (int i=0; i<this->blockSize_; ++i)
255  newind[i] = this->curDim_ + i;
256  lclV = MVT::CloneViewNonConst(*this->V_,newind);
257 
258  // Insert Delta at the end of V
259  MVT::SetBlock(*Delta,newind,*this->V_);
260  this->curDim_ += this->blockSize_;
261 
262  // Project the auxVecs out of Delta
263  if(this->auxVecs_.size() > 0)
264  this->orthman_->projectMat(*lclV,this->auxVecs_);
265 
266  // Update KV
267  if(this->Op_ != Teuchos::null)
268  {
269  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
270  Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
271  #endif
272  this->count_ApplyOp_+= this->blockSize_;
273 
274  lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
275  OPT::Apply(*this->Op_,*lclV,*lclKV);
276  }
277 
278  // Project out the components of Delta in the direction of V
279 
280  // gamma = KauxVecs' lclKV
281  int nauxVecs = MVT::GetNumberVecs(*projVecs);
282  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(nauxVecs,this->blockSize_));
283 
284  this->orthman_->innerProdMat(*KprojVecs,*lclKV,*gamma);
285 
286  // lclKV = lclKV - KauxVecs gamma
287  MVT::MvTimesMatAddMv(-this->ONE,*KprojVecs,*gamma,this->ONE,*lclKV);
288 
289  // lclV = lclV - auxVecs gamma
290  MVT::MvTimesMatAddMv(-this->ONE,*projVecs,*gamma,this->ONE,*lclV);
291 
292  // Normalize lclKV
293  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma2 = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
294  rank = this->orthman_->normalizeMat(*lclKV,gamma2);
295 
296  // lclV = lclV/gamma
297  Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
298  SDsolver.setMatrix(gamma2);
299  SDsolver.invert();
300  RCP<MV> tempMV = MVT::CloneCopy(*lclV);
301  MVT::MvTimesMatAddMv(this->ONE,*tempMV,*gamma2,this->ZERO,*lclV);
302 
303  TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
304  "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
305 
306  // Update MV
307  if(this->hasM_)
308  {
309  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
310  Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
311  #endif
312  this->count_ApplyM_+= this->blockSize_;
313 
314  lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
315  OPT::Apply(*this->MOp_,*lclV,*lclMV);
316  }
317  }
318 
319 }} // End of namespace Anasazi
320 
321 #endif
322 
323 // End of file AnasaziTraceMinDavidson.hpp
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&#39;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 &params)
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&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi&#39;s solvers.