dune-pdelab  2.0.0
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/parallel/mpihelper.hh>
8 
10 
11 #include "solver.hh"
12 
13 #if HAVE_EIGEN
14 
15 #include <Eigen/Eigen>
16 
17 namespace Dune {
18  namespace PDELab {
19 
20  //==============================================================================
21  // Here we add some standard linear solvers conforming to the linear solver
22  // interface required to solve linear and nonlinear problems.
23  //==============================================================================
24 
25  template<class PreconditionerImp>
26  class EigenBackend_BiCGSTAB_Base
27  : public LinearResultStorage
28  {
29  public:
35  explicit EigenBackend_BiCGSTAB_Base(unsigned maxiter_=5000)
36  : maxiter(maxiter_)
37  {}
38 
46  template<class M, class V, class W>
47  void apply(M& A, V& z, W& r, typename W::field_type reduction)
48  {
49  // PreconditionerImp preconditioner;
50  typedef typename M::Container Mat;
51  Eigen::BiCGSTAB<Mat, PreconditionerImp> solver;
52  solver.setMaxIterations(maxiter);
53  solver.setTolerance(reduction);
54  Dune::Timer watch;
55  watch.reset();
56  solver.compute(A.base());
57  z.base() = solver.solve(r.base());
58  double elapsed = watch.elapsed();
59 
60  res.converged = solver.info() == Eigen::ComputationInfo::Success;
61  res.iterations = solver.iterations();
62  res.elapsed = elapsed;
63  res.reduction = solver.error();
64  res.conv_rate = 0;
65  }
66 
67  public:
68  template<class V>
69  typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(const V& v) const
70  {
71  return v.base().norm();
72  }
73 
74  private:
75  unsigned maxiter;
76  int verbose;
77  };
78 
79  class EigenBackend_BiCGSTAB_IILU
80  : public EigenBackend_BiCGSTAB_Base<Eigen::IncompleteLUT<double> >
81  {
82  public:
83  explicit EigenBackend_BiCGSTAB_IILU(unsigned maxiter_=5000)
84  : EigenBackend_BiCGSTAB_Base(maxiter_)
85  {
86  }
87  };
88 
89  class EigenBackend_BiCGSTAB_Diagonal
90  : public EigenBackend_BiCGSTAB_Base<Eigen::DiagonalPreconditioner<double> >
91  {
92  public:
93  explicit EigenBackend_BiCGSTAB_Diagonal(unsigned maxiter_=5000)
94  : EigenBackend_BiCGSTAB_Base(maxiter_)
95  {}
96  };
97 
98  template< class Preconditioner, int UpLo >
99  class EigenBackend_CG_Base
100  : public SequentialNorm, public LinearResultStorage
101  {
102  public:
108  explicit EigenBackend_CG_Base(unsigned maxiter_=5000)
109  : maxiter(maxiter_)
110  {}
111 
119  template<class M, class V, class W>
120  void apply(M& A, V& z, W& r, typename W::field_type reduction)
121  {
122  typedef typename M::Container Mat;
123  Eigen::ConjugateGradient<Mat, UpLo, Preconditioner> solver;
124  solver.setMaxIterations(maxiter);
125  solver.setTolerance(reduction);
126  Dune::Timer watch;
127  watch.reset();
128  solver.compute(A.base());
129  z.base() = solver.solve(r.base());
130  double elapsed = watch.elapsed();
131 
132 
133  res.converged = solver.info() == Eigen::ComputationInfo::Success;
134  res.iterations = solver.iterations();
135  res.elapsed = elapsed;
136  res.reduction = solver.error();
137  res.conv_rate = 0;
138  }
139 
140  private:
141  unsigned maxiter;
142  int verbose;
143  };
144 
145 
146  class EigenBackend_CG_IILU_Up
147  : public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Upper >
148  {
149  public:
150  explicit EigenBackend_CG_IILU_Up(unsigned maxiter_=5000)
151  : EigenBackend_CG_Base(maxiter_)
152  {}
153  };
154 
155  class EigenBackend_CG_Diagonal_Up
156  : public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Upper >
157  {
158  public:
159  explicit EigenBackend_CG_Diagonal_Up(unsigned maxiter_=5000)
160  : EigenBackend_CG_Base(maxiter_)
161  {}
162  };
163 
164  class EigenBackend_CG_IILU_Lo
165  : public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Lower >
166  {
167  public:
168  explicit EigenBackend_CG_IILU_Lo(unsigned maxiter_=5000)
169  : EigenBackend_CG_Base(maxiter_)
170  {}
171  };
172 
173  class EigenBackend_CG_Diagonal_Lo
174  : public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Lower >
175  {
176  public:
177  explicit EigenBackend_CG_Diagonal_Lo(unsigned maxiter_=5000)
178  : EigenBackend_CG_Base(maxiter_)
179  {}
180  };
181 
182  template<template<class, int> class Solver, int UpLo>
183  class EigenBackend_SPD_Base
184  : public SequentialNorm, public LinearResultStorage
185  {
186  public:
192  explicit EigenBackend_SPD_Base()
193  {}
194 
202  template<class M, class V, class W>
203  void apply(M& A, V& z, W& r, typename W::field_type reduction)
204  {
205  typedef typename M::Container Mat;
206  Solver<Mat, UpLo> solver;
207  Dune::Timer watch;
208  watch.reset();
209  solver.compute(A.base());
210  z.base() = solver.solve(r.base());
211  double elapsed = watch.elapsed();
212 
213  res.converged = solver.info() == Eigen::ComputationInfo::Success;
214  res.iterations = solver.iterations();
215  res.elapsed = elapsed;
216  res.reduction = solver.error();
217  res.conv_rate = 0;
218  }
219  private:
220  };
221 
222  class EigenBackend_SimplicialCholesky_Up
223  : public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Upper >
224  {
225  public:
226  explicit EigenBackend_SimplicialCholesky_Up()
227  {}
228  };
229 
230  class EigenBackend_SimplicialCholesky_Lo
231  : public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Lower >
232  {
233  public:
234  explicit EigenBackend_SimplicialCholesky_Lo()
235  {}
236  };
237 
238 /* class EigenBackend_SimplicialLDLt_Up
239  * : public EigenBackend_SPD_Base<Eigen::SimplicialLDLt, Eigen::Upper >
240  * {
241  * public:
242  * explicit EigenBackend_SimplicialLDLt_Up()
243  * {}
244  * };
245 
246  * class EigenBackend_SimplicialLDLt_Lo
247  * : public EigenBackend_SPD_Base<Eigen::SimplicialLDLt, Eigen::Lower >
248  * {
249  * public:
250  * explicit EigenBackend_SimplicialLDLt_Lo()
251  * {}
252  * };
253 
254  * class EigenBackend_SimplicialLLt_Up
255  * : public EigenBackend_SPD_Base<Eigen::SimplicialLLt, Eigen::Upper >
256  * {
257  * public:
258  * explicit EigenBackend_SimplicialLLt_Up()
259  * {}
260  * };
261 
262  * class EigenBackend_SimplicialLLt_Lo
263  * : public EigenBackend_SPD_Base<Eigen::SimplicialLLt, Eigen::Lower >
264  * {
265  * public:
266  * explicit EigenBackend_SimplicialLLt_Lo()
267  * {}
268  * };
269 
270  * class EigenBackend_Cholmod_Up
271  * : public EigenBackend_SPD_Base<Eigen::CholmodDecomposition, Eigen::Upper >
272  * {
273  * public:
274  * explicit EigenBackend_Cholmod_Up()
275  * {}
276  * };
277 
278  * class EigenBackend_Cholmod_Lo
279  * : public EigenBackend_SPD_Base<Eigen::CholmodDecomposition, Eigen::Lower >
280  * {
281  * public:
282  * explicit EigenBackend_Cholmod_Lo()
283  * {}
284  * };*/
285 
286  class EigenBackend_LeastSquares
287  : public SequentialNorm, public LinearResultStorage
288  {
289  private:
290  const static unsigned int defaultFlag = Eigen::ComputeThinU | Eigen::ComputeThinV;
291  public:
297  explicit EigenBackend_LeastSquares(unsigned int flags = defaultFlag)
298  : flags_(flags)
299  {}
300 
308  template<class M, class V, class W>
309  void apply(M& A, V& z, W& r, typename W::field_type reduction)
310  {
311  Dune::Timer watch;
312  watch.reset();
313  typedef typename M::Container Mat;
314  Eigen::JacobiSVD<Mat, Eigen::ColPivHouseholderQRPreconditioner> solver(A, flags_);
315  z.base() = solver.solve(r.base());
316  double elapsed = watch.elapsed();
317 
318  res.converged = solver.info() == Eigen::ComputationInfo::Success;
319  res.iterations = solver.iterations();
320  res.elapsed = elapsed;
321  res.reduction = solver.error();
322  res.conv_rate = 0;
323  }
324  private:
325  unsigned int flags_;
326  };
327 
328  } // namespace PDELab
329 } // namespace Dune
330 
331 #endif // HAVE_EIGEN
332 
333 #endif // DUNE_PDELAB_EIGENSOLVERBACKEND