Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTraceMin.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//
42// Anasazi: Block Eigensolvers Package
43//
44
49#ifndef ANASAZI_TRACEMIN_HPP
50#define ANASAZI_TRACEMIN_HPP
51
52#include "AnasaziTypes.hpp"
53#include "AnasaziBasicSort.hpp"
55
56#ifdef HAVE_ANASAZI_EPETRA
57 #include "Epetra_Operator.h"
58#endif
59
63#include "Teuchos_ScalarTraits.hpp"
64
67
68#include "Teuchos_LAPACK.hpp"
69#include "Teuchos_BLAS.hpp"
70#include "Teuchos_SerialDenseMatrix.hpp"
71#include "Teuchos_SerialDenseSolver.hpp"
72#include "Teuchos_ParameterList.hpp"
73#include "Teuchos_TimeMonitor.hpp"
74
75// TODO: TraceMin performs some unnecessary computations upon restarting. Fix it!
76
77namespace Anasazi {
78namespace Experimental {
130 template <class ScalarType, class MV, class OP>
131 class TraceMin : public TraceMinBase<ScalarType,MV,OP> {
132 public:
134
135
176 TraceMin( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
177 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
178 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
179 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
180 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
181 Teuchos::ParameterList &params
182 );
183
184 private:
185 //
186 // Convenience typedefs
187 //
191 typedef Teuchos::ScalarTraits<ScalarType> SCT;
192 typedef typename SCT::magnitudeType MagnitudeType;
193 const MagnitudeType ONE;
194 const MagnitudeType ZERO;
195 const MagnitudeType NANVAL;
196
197 // TraceMin specific methods
198 void addToBasis(const Teuchos::RCP<const MV> Delta);
199
200 void harmonicAddToBasis(const Teuchos::RCP<const MV> Delta);
201 };
202
205 //
206 // Implementations
207 //
210
211
213 // Constructor
214 template <class ScalarType, class MV, class OP>
216 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
217 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
218 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
219 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
220 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
221 Teuchos::ParameterList &params
222 ) :
223 TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params),
224 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
225 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
226 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan())
227 {
228 }
229
230
231 template <class ScalarType, class MV, class OP>
232 void TraceMin<ScalarType,MV,OP>::addToBasis(const Teuchos::RCP<const MV> Delta)
233 {
234 MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);
235
236 if(this->hasM_)
237 {
238#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
239 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
240#endif
241 this->count_ApplyM_+= this->blockSize_;
242
243 OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
244 }
245
246 int rank;
247 {
248#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
249 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
250#endif
251
252 if(this->numAuxVecs_ > 0)
253 {
254 rank = this->orthman_->projectAndNormalizeMat(*this->V_,this->auxVecs_,
255 Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
256 Teuchos::null,this->MV_,this->MauxVecs_);
257 }
258 else
259 {
260 rank = this->orthman_->normalizeMat(*this->V_,Teuchos::null,this->MV_);
261 }
262 }
263
264 // FIXME (mfh 07 Oct 2014) This variable is currently unused, but
265 // it would make sense to use it to check whether the block is
266 // rank deficient.
267 (void) rank;
268
269 if(this->Op_ != Teuchos::null)
270 {
271#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
272 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
273#endif
274 this->count_ApplyOp_+= this->blockSize_;
275 OPT::Apply(*this->Op_,*this->V_,*this->KV_);
276 }
277 }
278
279
280
281 template <class ScalarType, class MV, class OP>
282 void TraceMin<ScalarType,MV,OP>::harmonicAddToBasis(const Teuchos::RCP<const MV> Delta)
283 {
284 // V = X - Delta
285 MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);
286
287 // Project out auxVecs
288 if(this->numAuxVecs_ > 0)
289 {
290#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
291 Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
292#endif
293 this->orthman_->projectMat(*this->V_,this->auxVecs_);
294 }
295
296 // Compute KV
297 if(this->Op_ != Teuchos::null)
298 {
299#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
300 Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
301#endif
302 this->count_ApplyOp_+= this->blockSize_;
303
304 OPT::Apply(*this->Op_,*this->V_,*this->KV_);
305 }
306
307 // Normalize lclKV
308 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
309 int rank = this->orthman_->normalizeMat(*this->KV_,gamma);
310 // FIXME (mfh 18 Feb 2015) It would make sense to check the rank.
311 (void) rank;
312
313 // lclV = lclV/gamma
314 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
315 SDsolver.setMatrix(gamma);
316 SDsolver.invert();
317 RCP<MV> tempMV = MVT::CloneCopy(*this->V_);
318 MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*this->V_);
319
320 // Compute MV
321 if(this->hasM_)
322 {
323#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
324 Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
325#endif
326 this->count_ApplyM_+= this->blockSize_;
327
328 OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
329 }
330 }
331
332}} // End of namespace Anasazi
333
334#endif
335
336// End of file AnasaziTraceMin.hpp
Basic implementation of the Anasazi::SortManager class.
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.
Class which provides internal utilities for the Anasazi solvers.
Abstract base class for trace minimization eigensolvers.
Types and exceptions used within Anasazi solvers and interfaces.
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 iteration, a preconditioned iteration for solving linear symmetric p...
TraceMin(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)
TraceMin 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, static class providing utilities for the solvers.
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.