3 #ifndef DUNE_PDELAB_EIGEN_SOLVERS_HH
4 #define DUNE_PDELAB_EIGEN_SOLVERS_HH
6 #include <dune/common/deprecated.hh>
7 #include <dune/common/parallel/mpihelper.hh>
15 #include <Eigen/Eigen>
25 template<
class PreconditionerImp>
26 class EigenBackend_BiCGSTAB_Base
27 :
public LinearResultStorage
35 explicit EigenBackend_BiCGSTAB_Base(
unsigned maxiter_=5000)
46 template<
class M,
class V,
class W>
47 void apply(M& A, V& z, W& r,
typename W::field_type reduction)
50 typedef typename M::Container Mat;
51 Eigen::BiCGSTAB<Mat, PreconditionerImp> solver;
52 solver.setMaxIterations(maxiter);
53 solver.setTolerance(reduction);
56 solver.compute(A.base());
57 z.base() = solver.solve(r.base());
58 double elapsed = watch.elapsed();
60 res.converged = solver.info() == Eigen::ComputationInfo::Success;
61 res.iterations = solver.iterations();
62 res.elapsed = elapsed;
63 res.reduction = solver.error();
69 typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(
const V& v)
const
71 return v.base().norm();
79 class EigenBackend_BiCGSTAB_IILU
80 :
public EigenBackend_BiCGSTAB_Base<Eigen::IncompleteLUT<double> >
83 explicit EigenBackend_BiCGSTAB_IILU(
unsigned maxiter_=5000)
84 : EigenBackend_BiCGSTAB_Base(maxiter_)
89 class EigenBackend_BiCGSTAB_Diagonal
90 :
public EigenBackend_BiCGSTAB_Base<Eigen::DiagonalPreconditioner<double> >
93 explicit EigenBackend_BiCGSTAB_Diagonal(
unsigned maxiter_=5000)
94 : EigenBackend_BiCGSTAB_Base(maxiter_)
98 template<
class Preconditioner,
int UpLo >
99 class EigenBackend_CG_Base
100 :
public SequentialNorm,
public LinearResultStorage
108 explicit EigenBackend_CG_Base(
unsigned maxiter_=5000)
119 template<
class M,
class V,
class W>
120 void apply(M& A, V& z, W& r,
typename W::field_type reduction)
122 typedef typename M::Container Mat;
123 Eigen::ConjugateGradient<Mat, UpLo, Preconditioner> solver;
124 solver.setMaxIterations(maxiter);
125 solver.setTolerance(reduction);
128 solver.compute(A.base());
129 z.base() = solver.solve(r.base());
130 double elapsed = watch.elapsed();
133 res.converged = solver.info() == Eigen::ComputationInfo::Success;
134 res.iterations = solver.iterations();
135 res.elapsed = elapsed;
136 res.reduction = solver.error();
146 class EigenBackend_CG_IILU_Up
147 :
public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Upper >
150 explicit EigenBackend_CG_IILU_Up(
unsigned maxiter_=5000)
151 : EigenBackend_CG_Base(maxiter_)
155 class EigenBackend_CG_Diagonal_Up
156 :
public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Upper >
159 explicit EigenBackend_CG_Diagonal_Up(
unsigned maxiter_=5000)
160 : EigenBackend_CG_Base(maxiter_)
164 class EigenBackend_CG_IILU_Lo
165 :
public EigenBackend_CG_Base<Eigen::IncompleteLUT<double>, Eigen::Lower >
168 explicit EigenBackend_CG_IILU_Lo(
unsigned maxiter_=5000)
169 : EigenBackend_CG_Base(maxiter_)
173 class EigenBackend_CG_Diagonal_Lo
174 :
public EigenBackend_CG_Base<Eigen::DiagonalPreconditioner<double>, Eigen::Lower >
177 explicit EigenBackend_CG_Diagonal_Lo(
unsigned maxiter_=5000)
178 : EigenBackend_CG_Base(maxiter_)
182 template<
template<
class,
int>
class Solver,
int UpLo>
183 class EigenBackend_SPD_Base
184 :
public SequentialNorm,
public LinearResultStorage
192 explicit EigenBackend_SPD_Base()
202 template<
class M,
class V,
class W>
203 void apply(M& A, V& z, W& r,
typename W::field_type reduction)
205 typedef typename M::Container Mat;
206 Solver<Mat, UpLo> solver;
209 solver.compute(A.base());
210 z.base() = solver.solve(r.base());
211 double elapsed = watch.elapsed();
213 res.converged = solver.info() == Eigen::ComputationInfo::Success;
214 res.iterations = solver.iterations();
215 res.elapsed = elapsed;
216 res.reduction = solver.error();
222 class EigenBackend_SimplicialCholesky_Up
223 :
public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Upper >
226 explicit EigenBackend_SimplicialCholesky_Up()
230 class EigenBackend_SimplicialCholesky_Lo
231 :
public EigenBackend_SPD_Base<Eigen::SimplicialCholesky, Eigen::Lower >
234 explicit EigenBackend_SimplicialCholesky_Lo()
286 class EigenBackend_LeastSquares
287 :
public SequentialNorm,
public LinearResultStorage
290 const static unsigned int defaultFlag = Eigen::ComputeThinU | Eigen::ComputeThinV;
297 explicit EigenBackend_LeastSquares(
unsigned int flags = defaultFlag)
308 template<
class M,
class V,
class W>
309 void apply(M& A, V& z, W& r,
typename W::field_type reduction)
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();
318 res.converged = solver.info() == Eigen::ComputationInfo::Success;
319 res.iterations = solver.iterations();
320 res.elapsed = elapsed;
321 res.reduction = solver.error();
333 #endif // DUNE_PDELAB_EIGENSOLVERBACKEND