dune-fem  2.8-git
ldlsolver.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_LDLSOLVER_HH
2 #define DUNE_FEM_LDLSOLVER_HH
3 
4 #include <algorithm>
5 #include <iostream>
6 #include <tuple>
7 #include <vector>
8 #include <utility>
9 
10 #include <dune/common/exceptions.hh>
11 #include <dune/common/hybridutilities.hh>
15 #include <dune/fem/io/parameter.hh>
19 
20 #if HAVE_SUITESPARSE_LDL
21 
22 #ifdef __cplusplus
23 extern "C"
24 {
25 #include "ldl.h"
26 #include "amd.h"
27 }
28 #endif
29 
30 namespace Dune
31 {
32 namespace Fem
33 {
34 
45 template< class DiscreteFunction,
46  class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
47 class LDLInverseOperator;
48 
49 template< class DiscreteFunction, class Matrix >
50 struct LDLInverseOperatorTraits
51 {
52  typedef DiscreteFunction DiscreteFunctionType;
53  typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
54 
56  typedef OperatorType PreconditionerType;
57 
58  typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
59 
60  typedef LDLInverseOperator< DiscreteFunction, Matrix > InverseOperatorType;
61  typedef SolverParameter SolverParameterType;
62 };
63 
64 
65 
73 template< class DF, class Matrix >
74 class LDLInverseOperator : public InverseOperatorInterface< LDLInverseOperatorTraits< DF, Matrix > >
75 {
76  typedef LDLInverseOperatorTraits< DF, Matrix > Traits;
77  typedef InverseOperatorInterface< Traits > BaseType;
78 
79  friend class InverseOperatorInterface< Traits >;
80 public:
81  typedef LDLInverseOperator< DF, Matrix > ThisType;
82 
83  typedef typename BaseType :: SolverDiscreteFunctionType
84  SolverDiscreteFunctionType;
85 
86  typedef typename BaseType :: OperatorType OperatorType;
87  typedef typename BaseType :: AssembledOperatorType AssembledOperatorType;
88 
89  // \brief The column-compressed matrix type.
90  typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,int> CCSMatrixType;
91 
92  typedef typename SolverDiscreteFunctionType::DofType DofType;
93  typedef typename SolverDiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
94 
101  LDLInverseOperator(const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
102  const ParameterReader &parameter = Parameter::container() ) :
103  verbose_(verbose && Dune::Fem::Parameter::verbose()), ccsmat_()
104  {}
105 
111  LDLInverseOperator(const double& redEps, const double& absLimit, const int& maxIter,
112  const ParameterReader &parameter = Parameter::container() ) :
113  verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_()
114  {}
115 
116  LDLInverseOperator(const SolverParameter &parameter = SolverParameter(Parameter::container()) )
117  : BaseType(parameter), verbose_(BaseType::verbose())
118  {}
119 
127  LDLInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
128  const ParameterReader &parameter = Parameter::container() ) :
129  verbose_(verbose), ccsmat_(), isloaded_(false)
130  {
131  bind(op);
132  }
133 
140  LDLInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter,
141  const ParameterReader &parameter = Parameter::container() ) :
142  verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
143  {
144  bind(op);
145  }
146 
147  LDLInverseOperator(const OperatorType& op, const ParameterReader &parameter = Parameter::container() ) :
148  verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
149  {
150  bind(op);
151  }
152 
153  // \brief Destructor.
154  ~LDLInverseOperator()
155  {
156  unbind();
157  }
158 
159  void bind( const OperatorType& op )
160  {
161  // clear old storage
162  finalize();
163  BaseType::bind( op );
164  }
165 
166  void unbind ()
167  {
168  finalize();
169  BaseType::unbind();
170  }
171 
172  // \brief Decompose matrix.
173  template<typename... A>
174  void prepare(A... ) const
175  {
176  if( ! assembledOperator_ )
177  DUNE_THROW(NotImplemented,"LDLInverseOperator only works for assembled systems!");
178 
179  if(!isloaded_)
180  {
181  ccsmat_ = assembledOperator_->exportMatrix();
182  decompose();
183  isloaded_ = true;
184  }
185  }
186 
187  // \brief Free allocated memory.
188  virtual void finalize()
189  {
190  if(isloaded_)
191  {
192  ccsmat_.free();
193  delete [] D_;
194  delete [] Y_;
195  delete [] Lp_;
196  delete [] Lx_;
197  delete [] Li_;
198  delete [] P_;
199  delete [] Pinv_;
200  isloaded_ = false;
201  }
202  }
203 
210  template<typename... DFs>
211  void apply(const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest) const
212  {
213  // copy DOF's arg into a consecutive vector
214  std::vector<DofType> vecArg(arg.size());
215  auto vecArgIt(vecArg.begin());
216  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},
217  [&](auto i){vecArgIt=std::copy(std::get<i>(arg).dbegin(),std::get<i>(arg).dend(),vecArgIt);});
218  std::vector<DofType> vecDest(dest.size());
219  // apply operator
220  apply(vecArg.data(),vecDest.data());
221  // copy back solution into dest
222  auto vecDestIt(vecDest.begin());
223  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},[&](auto i){for(auto& dof:dofs(std::get<i>(dest))) dof=(*(vecDestIt++));});
224  }
225 
226 protected:
227  void printTexInfo(std::ostream& out) const
228  {
229  out<<"Solver: LDL direct solver"<<std::endl;
230  }
231 
232  // \brief Print some statistics about the LDL decomposition.
233  void printDecompositionInfo() const
234  {
235  amd_info(info_);
236  }
237 
241  DofType* getD()
242  {
243  return D_;
244  }
245 
249  int* getLp()
250  {
251  return Lp_;
252  }
253 
257  int* getLi()
258  {
259  return Li_;
260  }
261 
265  DofType* getLx()
266  {
267  return Lx_;
268  }
269 
273  CCSMatrixType& getCCSMatrix()
274  {
275  return ccsmat_;
276  }
277 
278 protected:
285  void apply(const DofType* arg, DofType* dest) const
286  {
287  prepare();
288 
289  // apply part of the call
290  const std::size_t dimMat(ccsmat_.N());
291  ldl_perm(dimMat, Y_, const_cast<DofType*>(arg), P_);
292  ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
293  ldl_dsolve(dimMat, Y_, D_);
294  ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
295  ldl_permt(dimMat, dest, Y_, P_);
296 
297  const_cast<ThisType*>(this)->finalize();
298  }
299 
306  int apply(const SolverDiscreteFunctionType& arg,
307  SolverDiscreteFunctionType& dest) const
308  {
309  apply(arg.leakPointer(),dest.leakPointer());
310  return 0;
311  }
312 
313 
314 protected:
315  using BaseType :: assembledOperator_;
316  const bool verbose_;
317 
318  mutable CCSMatrixType ccsmat_;
319  mutable bool isloaded_ = false;
320  mutable int* Lp_;
321  mutable int* Parent_;
322  mutable int* Lnz_;
323  mutable int* Flag_;
324  mutable int* Pattern_;
325  mutable int* P_;
326  mutable int* Pinv_;
327  mutable DofType* D_;
328  mutable DofType* Y_;
329  mutable DofType* Lx_;
330  mutable int* Li_;
331  mutable double info_[AMD_INFO];
332 
333  // \brief Computes the LDL decomposition.
334  void decompose() const
335  {
336  // allocate vectors
337  const std::size_t dimMat(ccsmat_.N());
338  D_ = new DofType [dimMat];
339  Y_ = new DofType [dimMat];
340  Lp_ = new int [dimMat + 1];
341  Parent_ = new int [dimMat];
342  Lnz_ = new int [dimMat];
343  Flag_ = new int [dimMat];
344  Pattern_ = new int [dimMat];
345  P_ = new int [dimMat];
346  Pinv_ = new int [dimMat];
347 
348  if(amd_order (dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), P_, (DofType *) NULL, info_) < AMD_OK)
349  DUNE_THROW(InvalidStateException,"LDL Error: AMD failed!");
350  if(verbose_)
351  printDecompositionInfo();
352  // compute the symbolic factorisation
353  ldl_symbolic(dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
354  // initialise those entries of additionalVectors_ whose dimension is known only now
355  Lx_ = new DofType [Lp_[dimMat]];
356  Li_ = new int [Lp_[dimMat]];
357  // compute the numeric factorisation
358  const std::size_t k(ldl_numeric(dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), ccsmat_.getValues(),
359  Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
360  // free temporary vectors
361  delete [] Flag_;
362  delete [] Pattern_;
363  delete [] Parent_;
364  delete [] Lnz_;
365 
366  if(k!=dimMat)
367  {
368  std::cerr<<"LDL Error: D("<<k<<","<<k<<") is zero!"<<std::endl;
369  DUNE_THROW(InvalidStateException,"LDL Error: factorisation failed!");
370  }
371  }
372 };
373 
374 // deprecated old type
375 template<class DF, class Op, bool symmetric=false>
376 using LDLOp = LDLInverseOperator< DF, typename Op::MatrixType >;
377 
378 } // end namespace Fem
379 } // end namespace Dune
380 
381 #endif // #if HAVE_SUITESPARSE_LDL
382 
383 #endif // #ifndef DUNE_FEM_LDLSOLVER_HH
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static ParameterContainer & container()
Definition: io/parameter.hh:193