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