Anasazi Version of the Day
Loading...
Searching...
No Matches
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 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
46#ifndef ANASAZI_TRACEMIN_DAVIDSON_HPP
47#define ANASAZI_TRACEMIN_DAVIDSON_HPP
48
49#include "AnasaziConfigDefs.hpp"
55
56#include "Teuchos_ScalarTraits.hpp"
57#include "Teuchos_SerialDenseMatrix.hpp"
58#include "Teuchos_ParameterList.hpp"
59#include "Teuchos_TimeMonitor.hpp"
60
61
62namespace Anasazi {
63namespace Experimental {
64
79 template <class ScalarType, class MV, class OP>
80 class TraceMinDavidson : public TraceMinBase<ScalarType,MV,OP> {
81 public:
82
91 TraceMinDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
92 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
93 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
94 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
95 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
96 Teuchos::ParameterList &params
97 );
98
99 private:
100 //
101 // Convenience typedefs
102 //
105 typedef Teuchos::ScalarTraits<ScalarType> SCT;
106 typedef typename SCT::magnitudeType MagnitudeType;
107
108 // TraceMin specific methods
109 void addToBasis(const Teuchos::RCP<const MV> Delta);
110
111 void harmonicAddToBasis(const Teuchos::RCP<const MV> Delta);
112 };
113
116 //
117 // Implementations
118 //
121
122
123
125 // Constructor
126 template <class ScalarType, class MV, class OP>
128 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
129 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
130 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
131 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
132 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
133 Teuchos::ParameterList &params
134 ) :
135 TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params)
136 {
137 }
138
139
141 // 1. Project Delta so that V' M Delta = 0 and Q' M Delta = 0
142 // 2. Normalize Delta so that Delta' M Delta = I
143 // 3. Add Delta to the end of V: V = [V Delta]
144 // 4. Update KV and MV
145 template <class ScalarType, class MV, class OP>
146 void TraceMinDavidson<ScalarType,MV,OP>::addToBasis(Teuchos::RCP<const MV> Delta)
147 {
148 // TODO: We should also test the row length and map, etc...
149 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
150 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
151
152 int rank;
153 // Vector of indices
154 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
155 // Pointer to the meaningful parts of V, KV, and MV
156 Teuchos::RCP<MV> lclV, lclKV, lclMV;
157 // Holds the vectors we project against
158 Teuchos::Array< Teuchos::RCP<const MV> > projVecs = this->auxVecs_;
159
160 // Get the existing parts of the basis and add them to the list of things we project against
161 for(int i=0; i<this->curDim_; i++)
162 curind[i] = i;
163 lclV = MVT::CloneViewNonConst(*this->V_,curind);
164 projVecs.push_back(lclV);
165
166 // Get the new part of the basis (where we're going to insert Delta)
167 for (int i=0; i<this->blockSize_; ++i)
168 newind[i] = this->curDim_ + i;
169 lclV = MVT::CloneViewNonConst(*this->V_,newind);
170
171 // Insert Delta at the end of V
172 MVT::SetBlock(*Delta,newind,*this->V_);
173 this->curDim_ += this->blockSize_;
174
175 // Project out the components of Delta in the direction of V
176 if(this->hasM_)
177 {
178 // It is more efficient to provide the orthomanager with MV
179 Teuchos::Array< Teuchos::RCP<const MV> > MprojVecs = this->MauxVecs_;
180 lclMV = MVT::CloneViewNonConst(*this->MV_,curind);
181 MprojVecs.push_back(lclMV);
182
183 // Compute M*Delta
184 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
185 {
186 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
187 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
188 #endif
189 this->count_ApplyM_+= this->blockSize_;
190 OPT::Apply(*this->MOp_,*lclV,*lclMV);
191 }
192
193 {
194 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
195 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
196 #endif
197
198 // Project and normalize Delta
199 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs,
200 Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
201 Teuchos::null,lclMV,MprojVecs);
202 }
203
204 MprojVecs.pop_back();
205 }
206 else
207 {
208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
209 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
210 #endif
211
212 // Project and normalize Delta
213 rank = this->orthman_->projectAndNormalizeMat(*lclV,projVecs);
214 }
215
216 projVecs.pop_back();
217
218 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
219 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
220
221 // Update KV
222 if(this->Op_ != Teuchos::null)
223 {
224 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
225 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
226 #endif
227 this->count_ApplyOp_+= this->blockSize_;
228
229 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
230 OPT::Apply(*this->Op_,*lclV,*lclKV);
231 }
232 }
233
234
235
237 // 1. Project Delta so that V' M Delta = 0 and Q' M Delta = 0
238 // 2. Normalize Delta so that Delta' M Delta = I
239 // 3. Add Delta to the end of V: V = [V Delta]
240 // 4. Update KV and MV
241 template <class ScalarType, class MV, class OP>
242 void TraceMinDavidson<ScalarType,MV,OP>::harmonicAddToBasis(Teuchos::RCP<const MV> Delta)
243 {
244 // TODO: We should also test the row length and map, etc...
245 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*Delta) != this->blockSize_, std::invalid_argument,
246 "Anasazi::TraceMinDavidson::addToBasis(): Delta does not have blockSize_ columns");
247
248 int rank;
249 // Vector of indices
250 std::vector<int> curind(this->curDim_), newind(this->blockSize_);
251 // Pointer to the meaningful parts of V, KV, and MV
252 Teuchos::RCP<MV> lclV, lclKV, lclMV, projVecs, KprojVecs;
253 // Array of things we project against
254
255 // Get the existing parts of the basis and add them to the list of things we project against
256 for(int i=0; i<this->curDim_; i++)
257 curind[i] = i;
258 projVecs = MVT::CloneViewNonConst(*this->V_,curind);
259
260 if(this->Op_ != Teuchos::null)
261 {
262 lclKV = MVT::CloneViewNonConst(*this->KV_,curind);
263 KprojVecs = lclKV;
264 }
265
266 // Get the new part of the basis (where we're going to insert Delta)
267 for (int i=0; i<this->blockSize_; ++i)
268 newind[i] = this->curDim_ + i;
269 lclV = MVT::CloneViewNonConst(*this->V_,newind);
270
271 // Insert Delta at the end of V
272 MVT::SetBlock(*Delta,newind,*this->V_);
273 this->curDim_ += this->blockSize_;
274
275 // Project the auxVecs out of Delta
276 if(this->auxVecs_.size() > 0)
277 this->orthman_->projectMat(*lclV,this->auxVecs_);
278
279 // Update KV
280 if(this->Op_ != Teuchos::null)
281 {
282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
283 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
284 #endif
285 this->count_ApplyOp_+= this->blockSize_;
286
287 lclKV = MVT::CloneViewNonConst(*this->KV_,newind);
288 OPT::Apply(*this->Op_,*lclV,*lclKV);
289 }
290
291 // Project out the components of Delta in the direction of V
292
293 // gamma = KauxVecs' lclKV
294 int nauxVecs = MVT::GetNumberVecs(*projVecs);
295 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(nauxVecs,this->blockSize_));
296
297 this->orthman_->innerProdMat(*KprojVecs,*lclKV,*gamma);
298
299 // lclKV = lclKV - KauxVecs gamma
300 MVT::MvTimesMatAddMv(-this->ONE,*KprojVecs,*gamma,this->ONE,*lclKV);
301
302 // lclV = lclV - auxVecs gamma
303 MVT::MvTimesMatAddMv(-this->ONE,*projVecs,*gamma,this->ONE,*lclV);
304
305 // Normalize lclKV
306 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma2 = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
307 rank = this->orthman_->normalizeMat(*lclKV,gamma2);
308
309 // lclV = lclV/gamma
310 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
311 SDsolver.setMatrix(gamma2);
312 SDsolver.invert();
313 RCP<MV> tempMV = MVT::CloneCopy(*lclV);
314 MVT::MvTimesMatAddMv(this->ONE,*tempMV,*gamma2,this->ZERO,*lclV);
315
316 TEUCHOS_TEST_FOR_EXCEPTION(rank != this->blockSize_,TraceMinBaseOrthoFailure,
317 "Anasazi::TraceMinDavidson::addToBasis(): Couldn't generate basis of full rank.");
318
319 // Update MV
320 if(this->hasM_)
321 {
322 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
323 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
324 #endif
325 this->count_ApplyM_+= this->blockSize_;
326
327 lclMV = MVT::CloneViewNonConst(*this->MV_,newind);
328 OPT::Apply(*this->MOp_,*lclV,*lclMV);
329 }
330 }
331
332}} // End of namespace Anasazi
333
334#endif
335
336// End of file AnasaziTraceMinDavidson.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Abstract base class for trace minimization eigensolvers.
This class defines the interface required by an eigensolver and status test class to compute solution...
This is an abstract base class for the trace minimization eigensolvers.
This class implements a TraceMin-Davidson iteration for solving symmetric generalized eigenvalue prob...
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.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Output managers remove the need for the eigensolver to know any information about the required output...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.