dune-fem  2.8-git
eigen.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_EIGEN_HH
2 #define DUNE_FEM_SOLVER_EIGEN_HH
3 
4 #ifdef HAVE_EIGEN
5 
6 #include <Eigen/IterativeLinearSolvers>
7 
11 
12 namespace Dune
13 {
14 
15  namespace Fem
16  {
17  template <class DiscreteFunction>
18  struct EigenInverseOperatorTraits
19  {
20  typedef Fem::EigenLinearOperator< DiscreteFunction, DiscreteFunction> OperatorType;
21  typedef typename OperatorType::MatrixType::MatrixStorageType Matrix;
22  };
23 
24  template< class DiscreteFunction, class EigenOp >
25  class EigenInverseOperator
26  : public Fem::Operator< DiscreteFunction, DiscreteFunction >
27  {
28  typedef Fem::Operator< DiscreteFunction, DiscreteFunction > Base;
29 
30  public:
31  typedef typename Base::DomainFunctionType DomainFunction;
32  typedef typename Base::RangeFunctionType RangeFunction;
33 
34  typedef Fem::EigenLinearOperator< RangeFunction, DomainFunction > OperatorType;
35 
36  private:
37  typedef typename OperatorType::MatrixType::MatrixStorageType Matrix;
38 
39  public:
40  EigenInverseOperator ( double redEps, double absLimit, unsigned int maxIter, bool verbose,
41  const ParameterReader &parameter = Parameter::container() )
42  : solver_(std::make_unique<EigenOp>()), absLimit_( absLimit )
43  {}
44 
45  EigenInverseOperator ( double redEps, double absLimit, unsigned int maxIter = std::numeric_limits<int>::max(),
46  const ParameterReader &parameter = Parameter::container() )
47  : solver_(std::make_unique<EigenOp>()), absLimit_( absLimit )
48  {}
49 
50  EigenInverseOperator ( const OperatorType &op, double redEps, double absLimit, unsigned int maxIter, bool verbose,
51  const ParameterReader &parameter = Parameter::container() )
52  : solver_(std::make_unique<EigenOp>()), absLimit_( absLimit )
53  {
54  bind( op );
55  }
56 
57  EigenInverseOperator ( const OperatorType &op, double redEps, double absLimit, unsigned int maxIter,
58  const ParameterReader &parameter = Parameter::container() )
59  : solver_(std::make_unique<EigenOp>()), absLimit_( absLimit )
60  {
61  bind( op );
62  }
63 
64  EigenInverseOperator ( const OperatorType &op,
65  double reduction, double absLimit,
66  const SolverParameter &parameter )
67  : EigenInverseOperator( reduction, absLimit, std::numeric_limits<int>::max(), parameter.parameter() )
68  {}
69  EigenInverseOperator ( double reduction, double absLimit,
70  const SolverParameter &parameter )
71  : EigenInverseOperator( reduction, absLimit, std::numeric_limits<int>::max(), parameter.parameter() )
72  {}
73 
74  void bind ( const OperatorType &op )
75  {
76  op_ = &op;
77  setup();
78  }
79  void unbind() { op_ = nullptr; }
80 
81  virtual void operator() ( const DomainFunction &u, RangeFunction &w ) const
82  {
83  if ( op_ )
84  w.dofVector().array().coefficients() = solver_->solve( u.dofVector().array().coefficients() );
85  }
86 
87  unsigned int iterations () const { return solver_->iterations(); }
88  void setMaxIterations ( unsigned int ) {}
89 
90  protected:
91  void setup ()
92  {
93  assert( op_ );
94  solver_->setTolerance( absLimit_ );
95  solver_->analyzePattern( op_->matrix().data() );
96  solver_->factorize( op_->matrix().data() );
97  }
98 
99  const OperatorType *op_ = nullptr;
100  std::unique_ptr<EigenOp> solver_;
101  double absLimit_;
102  };
103 
104  template< class DiscreteFunction>
105  using EigenCGInverseOperator = EigenInverseOperator<DiscreteFunction,
106  Eigen::ConjugateGradient<typename EigenInverseOperatorTraits<DiscreteFunction>::Matrix> >;
107 
108  template< class DiscreteFunction>
109  using EigenBiCGStabInverseOperator = EigenInverseOperator<DiscreteFunction,
110  Eigen::BiCGSTAB<typename EigenInverseOperatorTraits<DiscreteFunction>::Matrix> >;
111 
112  } // namespace Fem
113 } // namespace Dune
114 
115 #endif
116 #endif // #ifndef DUNE_FEM_SOLVER_VIENNACL_HH
Definition: bindguard.hh:11
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static double max(const Double &v, const double p)
Definition: double.hh:398
static constexpr T max(T a)
Definition: utility.hh:77
static ParameterContainer & container()
Definition: io/parameter.hh:193