Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
solvereigen.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2007-07-04
7 
8  Copyright (C) 2007-2011 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but 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 Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __SolverEigen_H
30 #define __SolverEigen_H 1
31 
32 #include <boost/tuple/tuple.hpp>
33 #include <boost/parameter.hpp>
34 #include <feel/feelcore/parameter.hpp>
35 
36 #include <feel/feelalg/enums.hpp>
37 #include <feel/feelcore/traits.hpp>
38 #include <feel/feelalg/vector.hpp>
40 
41 namespace Feel
42 {
49 template<typename T>
51 {
52 public:
53 
54 
58 
59  typedef T value_type;
60  typedef typename type_traits<T>::real_type real_type;
61 
63  typedef boost::shared_ptr<solvereigen_type> solvereigen_ptrtype;
64 
65  typedef boost::tuple<size_type, size_type, std::vector<double> > solve_return_type;
66 
67 
69  typedef boost::shared_ptr<vector_type> vector_ptrtype;
71  typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
72 
73  typedef boost::tuple<real_type, real_type, vector_ptrtype> eigenpair_type;
74  typedef std::map<real_type, eigenpair_type> eigenmodes_type;
75 
77 
81 
85  SolverEigen();
86 
102  SolverEigen( po::variables_map const& vm, std::string const& prefix = "" );
103 
107  SolverEigen( SolverEigen const & );
108 
112  virtual ~SolverEigen();
113 
118  static boost::shared_ptr<SolverEigen<value_type> > build( const SolverPackage solver_package = SOLVERS_SLEPC );
119 
126  static boost::shared_ptr<SolverEigen<value_type> > build( po::variables_map const& vm,
127  std::string const& prefix = std::string() );
128 
133  virtual eigenpair_type eigenPair ( unsigned int i ) = 0;
134 
135  /*
136  * Returns the eigen modes in a map
137  */
138  virtual eigenmodes_type eigenModes () = 0;
139 
141 
145 
146 
148 
152 
156  std::string const& prefix() const
157  {
158  return M_prefix;
159  }
160 
165  {
166  return M_eigen_solver_type;
167  }
168 
173  {
174  return M_eigen_problem_type;
175  }
176 
181  {
182  return M_position_of_spectrum;
183  }
184 
189  {
190  return M_spectral_transform;
191  }
192 
197  {
198  return M_nev;
199  }
200 
205  {
206  return M_ncv;
207  }
208 
212  value_type tolerance() const
213  {
214  return M_tolerance;
215  }
216 
221  {
222  return M_maxit;
223  }
224 
226 
230 
231 
236  {
237  M_eigen_solver_type = est;
238  }
239 
244  {
245  M_eigen_problem_type = ept;
246  }
247 
252  {
254  }
255 
260  {
262  }
263 
268  void setTolerance( value_type tol )
269  {
270  M_tolerance = tol;
271  }
272 
277  void setMaxIterations( size_type maxiter )
278  {
279  M_maxit = maxiter;
280  }
281 
286  {
287  M_nev = nev;
288  }
289 
294  {
295  M_ncv = ncv;
296  }
297 
299 
303 
308  bool initialized () const
309  {
310  return M_is_initialized;
311  }
312 
316  virtual void clear () {}
317 
321  virtual void init () = 0;
322 
323 
324  BOOST_PARAMETER_MEMBER_FUNCTION( ( solve_return_type ),
325  solve,
326  tag,
327  ( required
328  ( matrixA,( sparse_matrix_ptrtype ) )
329  ( matrixB,( sparse_matrix_ptrtype ) ) )
330  ( optional
331  ( maxit,( size_type ), 1000 )
332  ( tolerance,( double ), 1e-11 )
333  )
334  )
335  {
336  this->setMaxIterations( maxit );
337  this->setTolerance( tolerance );
338  solve_return_type ret = solve( matrixA, matrixB, M_nev, M_ncv, this->tolerance(), this->maxIterations() );
339  return ret;
340  }
341 
347  virtual solve_return_type solve ( MatrixSparse<value_type> &matrix_A,
348  int nev,
349  int ncv,
350  const double tol,
351  const unsigned int m_its ) = 0;
352 
358  solve_return_type solve ( boost::shared_ptr<MatrixSparse<value_type> > &matrix_A,
359  int nev,
360  int ncv,
361  const double tol,
362  const unsigned int m_its )
363  {
364  return this->solve( *matrix_A, nev, ncv, tol, m_its );
365  }
366 
367 
373  virtual solve_return_type solve ( MatrixSparse<value_type> &matrix_A,
374  MatrixSparse<value_type> &matrix_B,
375  int nev,
376  int ncv,
377  const double tol,
378  const unsigned int m_its ) = 0;
379 
385  solve_return_type solve ( boost::shared_ptr<MatrixSparse<value_type> > &matrix_A,
386  boost::shared_ptr<MatrixSparse<value_type> > &matrix_B,
387  int nev,
388  int ncv,
389  const double tol,
390  const unsigned int m_its )
391  {
392  return this->solve( *matrix_A, *matrix_B, nev, ncv, tol, m_its );
393  }
394 
396 
397 
398 
399 
400 
401 protected:
402 
404  std::string M_prefix;
405 
410 
415 
420 
425 
430 
433 
436 
439 
441  value_type M_tolerance;
442 };
443 
462 BOOST_PARAMETER_MEMBER_FUNCTION( ( typename SolverEigen<double>::eigenmodes_type ),
463  eigs,
464  tag,
465  ( required
466  ( matrixA,( d_sparse_matrix_ptrtype ) )
467  ( matrixB,( d_sparse_matrix_ptrtype ) ) )
468  ( optional
469  ( nev, ( int ), 1 )
470  ( ncv, ( int ), 3 )
471  ( backend,( BackendType ), BACKEND_PETSC )
472  ( solver,( EigenSolverType ), KRYLOVSCHUR )
473  ( problem,( EigenProblemType ), GHEP )
474  ( transform,( SpectralTransformType ), SHIFT )
475  ( spectrum,( PositionOfSpectrum ), LARGEST_MAGNITUDE )
476  ( maxit,( size_type ), 1000 )
477  ( tolerance,( double ), 1e-11 )
478  ( verbose,( bool ), false )
479  )
480  )
481 {
482  typedef boost::shared_ptr<Vector<double> > vector_ptrtype;
483  //boost::shared_ptr<SolverEigen<double> > eigen = SolverEigen<double>::build( backend );
484  boost::shared_ptr<SolverEigen<double> > eigen = SolverEigen<double>::build();
485  eigen->setEigenSolverType( solver );
486  eigen->setEigenProblemType( problem );
487  eigen->setPositionOfSpectrum( spectrum );
488  eigen->setNumberOfEigenValues( nev );
489  eigen->setNumberOfEigenValuesConverged( ncv );
490  eigen->setMaxIterations( maxit );
491  eigen->setSpectralTransform( transform );
492  eigen->setTolerance( tolerance );
493 
494  LOG(INFO) << "number of eigen values = " << nev << "\n";
495  LOG(INFO) << "number of eigen values converged = " << ncv << "\n";
496  LOG(INFO) << "number of eigen value solver iterations = " << maxit << "\n";
497  LOG(INFO) << "eigenvalue tolerance = " << tolerance << "\n";
498 
499  unsigned int nconv, nits;
500  std::vector<double> err( ncv );
501  boost::tie( nconv, nits, err ) = eigen->solve( _matrixA=matrixA,
502  _matrixB=matrixB,
503  _maxit=maxit,
504  _tolerance=tolerance );
505 
506  if ( verbose )
507  {
508  std::for_each( err.begin(), err.end(), []( double e )
509  {
510  std::cout << "||A x - lambda B x ||/||x|| = " << e << "\n";
511  } );
512  }
513 
514  return eigen->eigenModes();
515 }
516 
517 } // Feel
518 #endif /* __SolverEigen_H */

Generated on Sun Oct 20 2013 08:25:05 for Feel++ by doxygen 1.8.4