dune-pdelab  2.4-dev
solvers.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_EIGEN_SOLVERS_HH
4 #define DUNE_PDELAB_EIGEN_SOLVERS_HH
5 
6 #include <dune/common/deprecated.hh>
7 #include <dune/common/timer.hh>
8 #include <dune/common/parallel/mpihelper.hh>
9 
11 
12 #include "../solver.hh"
13 
14 #if HAVE_EIGEN
15 
16 #include <Eigen/Eigen>
17 #include <Eigen/Sparse>
18 
19 namespace Dune {
20  namespace PDELab {
21 
22  //==============================================================================
23  // Here we add some standard linear solvers conforming to the linear solver
24  // interface required to solve linear and nonlinear problems.
25  //==============================================================================
26 
27  template<class PreconditionerImp>
28  class EigenBackend_BiCGSTAB_Base
29  : public LinearResultStorage
30  {
31  public:
37  explicit EigenBackend_BiCGSTAB_Base(unsigned maxiter_=5000)
38  : maxiter(maxiter_)
39  {}
40 
48  template<class M, class V, class W>
49  void apply(M& A, V& z, W& r, typename W::field_type reduction)
50  {
51  using Backend::Native;
52  using Backend::native;
53  // PreconditionerImp preconditioner;
54  using Mat = Native<M>;
55  Eigen::BiCGSTAB<Mat, PreconditionerImp> solver;
56  solver.setMaxIterations(maxiter);
57  solver.setTolerance(reduction);
58  Dune::Timer watch;
59  watch.reset();
60  solver.compute(native(A));
61  native(z) = solver.solve(native(r));
62  double elapsed = watch.elapsed();
63 
64  res.converged = solver.info() == Eigen::ComputationInfo::Success;
65  res.iterations = solver.iterations();
66  res.elapsed = elapsed;
67  res.reduction = solver.error();
68  res.conv_rate = 0;
69  }
70 
71  public:
72  template<class V>
73  typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(const V& v) const
74  {
75  return native(v).norm();
76  }
77 
78  private:
79  unsigned maxiter;
80  int verbose;
81  };
82 
83  class EigenBackend_BiCGSTAB_IILU
84  : public EigenBackend_BiCGSTAB_Base<Eigen::IncompleteLUT<double> >
85  {
86  public:
87  explicit EigenBackend_BiCGSTAB_IILU(unsigned maxiter_=5000)
88  : EigenBackend_BiCGSTAB_Base(maxiter_)
89  {
90  }
91  };
92 
93  class EigenBackend_BiCGSTAB_Diagonal
94  : public EigenBackend_BiCGSTAB_Base<Eigen::DiagonalPreconditioner<double> >
95  {
96  public:
97  explicit EigenBackend_BiCGSTAB_Diagonal(unsigned maxiter_=5000)
98  : EigenBackend_BiCGSTAB_Base(maxiter_)
99  {}
100  };
101 
102  template< class Preconditioner, int UpLo >
103  class EigenBackend_CG_Base
104  : public SequentialNorm, public LinearResultStorage
105  {
106  public:
112  explicit EigenBackend_CG_Base(unsigned maxiter_=5000)
113  : maxiter(maxiter_)
114  {}
115 
123  template<class M, class V, class W>
124  void apply(M& A, V& z, W& r, typename W::field_type reduction)
125  {
126  using Backend::native;
127  using Mat = Backend::Native<M>;
128  Eigen::ConjugateGradient<Mat, UpLo, Preconditioner> solver;
129  solver.setMaxIterations(maxiter);
130  solver.setTolerance(reduction);
131  Dune::Timer watch;
132  watch.reset();
133  solver.compute(native(A));
134  native(z) = solver.solve(native(r));
135  double elapsed = watch.elapsed();
136 
137 
138  res.converged = solver.info() == Eigen::ComputationInfo::Success;
139  res.iterations = solver.iterations();
140  res.elapsed = elapsed;
141  res.reduction = solver.error();
142  res.conv_rate = 0;
143  }
144 
145  private:
146  unsigned maxiter;
147  int verbose;
148  };
149 
150 
151  class EigenBackend_CG_IILU_Up
152  : public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Upper >
153  {
154  public:
155  explicit EigenBackend_CG_IILU_Up(unsigned maxiter_=5000)
156  : EigenBackend_CG_Base(maxiter_)
157  {}
158  };
159 
160  class EigenBackend_CG_Diagonal_Up
161  : public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Upper >
162  {
163  public:
164  explicit EigenBackend_CG_Diagonal_Up(unsigned maxiter_=5000)
165  : EigenBackend_CG_Base(maxiter_)
166  {}
167  };
168 
169  class EigenBackend_CG_IILU_Lo
170  : public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Lower >
171  {
172  public:
173  explicit EigenBackend_CG_IILU_Lo(unsigned maxiter_=5000)
174  : EigenBackend_CG_Base(maxiter_)
175  {}
176  };
177 
178  class EigenBackend_CG_Diagonal_Lo
179  : public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Lower >
180  {
181  public:
182  explicit EigenBackend_CG_Diagonal_Lo(unsigned maxiter_=5000)
183  : EigenBackend_CG_Base(maxiter_)
184  {}
185  };
186 
187 #if EIGEN_VERSION_AT_LEAST(3,2,2)
188  template<template<class, int, class> class Solver, int UpLo>
189 #else
190  template<template<class, int> class Solver, int UpLo>
191 #endif
192  class EigenBackend_SPD_Base
193  : public SequentialNorm, public LinearResultStorage
194  {
195  public:
201  explicit EigenBackend_SPD_Base()
202  {}
203 
211  template<class M, class V, class W>
212  void apply(M& A, V& z, W& r, typename W::field_type reduction)
213  {
214  using Backend::native;
215  using Mat = Backend::Native<M>;
216 #if EIGEN_VERSION_AT_LEAST(3,2,2)
217  // use the approximate minimum degree algorithm for the ordering in
218  // the solver. This reproduces the default ordering for the Cholesky
219  // type solvers.
220  Solver<Mat, UpLo, Eigen::AMDOrdering<typename Mat::Index> > solver;
221 #else
222  Solver<Mat, UpLo> solver;
223 #endif
224  Dune::Timer watch;
225  watch.reset();
226  solver.compute(native(A));
227  native(z) = solver.solve(native(r));
228  double elapsed = watch.elapsed();
229 
230  res.converged = solver.info() == Eigen::ComputationInfo::Success;
231  res.iterations = solver.iterations();
232  res.elapsed = elapsed;
233  res.reduction = solver.error();
234  res.conv_rate = 0;
235  }
236  private:
237  };
238 
239  class EigenBackend_SimplicialCholesky_Up
240  : public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Upper >
241  {
242  public:
243  explicit EigenBackend_SimplicialCholesky_Up()
244  {}
245  };
246 
247  class EigenBackend_SimplicialCholesky_Lo
248  : public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Lower >
249  {
250  public:
251  explicit EigenBackend_SimplicialCholesky_Lo()
252  {}
253  };
254 
255 /* class EigenBackend_SimplicialLDLt_Up
256  * : public EigenBackend_SPD_Base<Eigen::SimplicialLDLt, Eigen::Upper >
257  * {
258  * public:
259  * explicit EigenBackend_SimplicialLDLt_Up()
260  * {}
261  * };
262 
263  * class EigenBackend_SimplicialLDLt_Lo
264  * : public EigenBackend_SPD_Base<Eigen::SimplicialLDLt, Eigen::Lower >
265  * {
266  * public:
267  * explicit EigenBackend_SimplicialLDLt_Lo()
268  * {}
269  * };
270 
271  * class EigenBackend_SimplicialLLt_Up
272  * : public EigenBackend_SPD_Base<Eigen::SimplicialLLt, Eigen::Upper >
273  * {
274  * public:
275  * explicit EigenBackend_SimplicialLLt_Up()
276  * {}
277  * };
278 
279  * class EigenBackend_SimplicialLLt_Lo
280  * : public EigenBackend_SPD_Base<Eigen::SimplicialLLt, Eigen::Lower >
281  * {
282  * public:
283  * explicit EigenBackend_SimplicialLLt_Lo()
284  * {}
285  * };
286 
287  * class EigenBackend_Cholmod_Up
288  * : public EigenBackend_SPD_Base<Eigen::CholmodDecomposition, Eigen::Upper >
289  * {
290  * public:
291  * explicit EigenBackend_Cholmod_Up()
292  * {}
293  * };
294 
295  * class EigenBackend_Cholmod_Lo
296  * : public EigenBackend_SPD_Base<Eigen::CholmodDecomposition, Eigen::Lower >
297  * {
298  * public:
299  * explicit EigenBackend_Cholmod_Lo()
300  * {}
301  * };*/
302 
303  class EigenBackend_LeastSquares
304  : public SequentialNorm, public LinearResultStorage
305  {
306  private:
307  const static unsigned int defaultFlag = Eigen::ComputeThinU | Eigen::ComputeThinV;
308  public:
314  explicit EigenBackend_LeastSquares(unsigned int flags = defaultFlag)
315  : flags_(flags)
316  {}
317 
325  template<class M, class V, class W>
326  void apply(M& A, V& z, W& r, typename W::field_type reduction)
327  {
328  Dune::Timer watch;
329  watch.reset();
330  using Backend::native;
331  using Mat = Backend::Native<M>;
332  Eigen::JacobiSVD<Mat, Eigen::ColPivHouseholderQRPreconditioner> solver(A, flags_);
333  native(z) = solver.solve(native(r));
334  double elapsed = watch.elapsed();
335 
336  res.converged = solver.info() == Eigen::ComputationInfo::Success;
337  res.iterations = solver.iterations();
338  res.elapsed = elapsed;
339  res.reduction = solver.error();
340  res.conv_rate = 0;
341  }
342  private:
343  unsigned int flags_;
344  };
345 
346  } // namespace PDELab
347 } // namespace Dune
348 
349 #endif // HAVE_EIGEN
350 
351 #endif // DUNE_PDELAB_EIGENSOLVERBACKEND
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:198
Definition: adaptivity.hh:27
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper...
Definition: backend/interface.hh:182