dune-istl  2.8.0
ldl.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_LDL_HH
4 #define DUNE_ISTL_LDL_HH
5 
6 #if HAVE_SUITESPARSE_LDL || defined DOXYGEN
7 
8 #include <iostream>
9 #include <memory>
10 #include <type_traits>
11 
12 #ifdef __cplusplus
13 extern "C"
14 {
15 #include "ldl.h"
16 #include "amd.h"
17 }
18 #endif
19 
20 #include <dune/common/exceptions.hh>
21 
23 #include <dune/istl/solvers.hh>
24 #include <dune/istl/solvertype.hh>
26 
27 namespace Dune {
39  // forward declarations
40  template<class M, class T, class TM, class TD, class TA>
41  class SeqOverlappingSchwarz;
42 
43  template<class T, bool tag>
44  struct SeqOverlappingSchwarzAssemblerHelper;
45 
52  template<class Matrix>
53  class LDL
54  {};
55 
69  template<typename T, typename A, int n, int m>
70  class LDL<BCRSMatrix<FieldMatrix<T,n,m>,A > >
71  : public InverseOperator<BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >,
72  BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > >
73  {
74  public:
79  typedef Dune::ISTL::Impl::BCCSMatrix<T,int> LDLMatrix;
81  typedef ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A>, int> MatrixInitializer;
83  typedef Dune::BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > > domain_type;
85  typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;
86 
89  {
90  return SolverCategory::Category::sequential;
91  }
92 
102  LDL(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
103  {
104  //check whether T is a supported type
105  static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
106  setMatrix(matrix);
107  }
108 
118  LDL(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false), verbose_(verbose)
119  {
120  //check whether T is a supported type
121  static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
122  setMatrix(matrix);
123  }
124 
134  LDL(const Matrix& matrix, const ParameterTree& config)
135  : LDL(matrix, config.get<int>("verbose", 0))
136  {}
137 
139  LDL() : matrixIsLoaded_(false), verbose_(0)
140  {}
141 
143  virtual ~LDL()
144  {
145  if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
146  free();
147  }
148 
151  {
152  const int dimMat(ldlMatrix_.N());
153  ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
154  ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
155  ldl_dsolve(dimMat, Y_, D_);
156  ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
157  ldl_permt(dimMat, reinterpret_cast<double*>(&x[0]), Y_, P_);
158  // this is a direct solver
159  res.iterations = 1;
160  res.converged = true;
161  }
162 
164  virtual void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
165  {
166  apply(x,b,res);
167  }
168 
174  void apply(T* x, T* b)
175  {
176  const int dimMat(ldlMatrix_.N());
177  ldl_perm(dimMat, Y_, b, P_);
178  ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
179  ldl_dsolve(dimMat, Y_, D_);
180  ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
181  ldl_permt(dimMat, x, Y_, P_);
182  }
183 
184  void setOption([[maybe_unused]] unsigned int option, [[maybe_unused]] double value)
185  {}
186 
188  void setMatrix(const Matrix& matrix)
189  {
190  if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
191  free();
192 
193  if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
194  ldlMatrix_.free();
195  ldlMatrix_.setSize(MatrixDimension<Matrix>::rowdim(matrix),
197  ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
198 
199  copyToBCCSMatrix(initializer, matrix);
200 
201  decompose();
202  }
203 
204  template<class S>
205  void setSubMatrix(const Matrix& matrix, const S& rowIndexSet)
206  {
207  if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
208  free();
209 
210  if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
211  ldlMatrix_.free();
212 
213  ldlMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(matrix) / matrix.N(),
214  rowIndexSet.size()*MatrixDimension<Matrix>::coldim(matrix) / matrix.M());
215  ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
216 
217  copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(matrix,rowIndexSet));
218 
219  decompose();
220  }
221 
226  inline void setVerbosity(int v)
227  {
228  verbose_=v;
229  }
230 
236  {
237  return ldlMatrix_;
238  }
239 
244  void free()
245  {
246  delete [] D_;
247  delete [] Y_;
248  delete [] Lp_;
249  delete [] Lx_;
250  delete [] Li_;
251  delete [] P_;
252  delete [] Pinv_;
253  ldlMatrix_.free();
254  matrixIsLoaded_ = false;
255  }
256 
258  inline const char* name()
259  {
260  return "LDL";
261  }
262 
267  inline double* getD()
268  {
269  return D_;
270  }
271 
276  inline int* getLp()
277  {
278  return Lp_;
279  }
280 
285  inline int* getLi()
286  {
287  return Li_;
288  }
289 
294  inline double* getLx()
295  {
296  return Lx_;
297  }
298 
299  private:
300  template<class M,class X, class TM, class TD, class T1>
301  friend class SeqOverlappingSchwarz;
302 
303  friend struct SeqOverlappingSchwarzAssemblerHelper<LDL<Matrix>,true>;
304 
306  void decompose()
307  {
308  // allocate vectors
309  const int dimMat(ldlMatrix_.N());
310  D_ = new double [dimMat];
311  Y_ = new double [dimMat];
312  Lp_ = new int [dimMat + 1];
313  Parent_ = new int [dimMat];
314  Lnz_ = new int [dimMat];
315  Flag_ = new int [dimMat];
316  Pattern_ = new int [dimMat];
317  P_ = new int [dimMat];
318  Pinv_ = new int [dimMat];
319 
320  double Info [AMD_INFO];
321  if(amd_order (dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), P_, (double *) NULL, Info) < AMD_OK)
322  DUNE_THROW(InvalidStateException,"Error: AMD failed!");
323  if(verbose_ > 0)
324  amd_info (Info);
325  // compute the symbolic factorisation
326  ldl_symbolic(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
327  // initialise those entries of additionalVectors_ whose dimension is known only now
328  Lx_ = new double [Lp_[dimMat]];
329  Li_ = new int [Lp_[dimMat]];
330  // compute the numeric factorisation
331  const int rank(ldl_numeric(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), ldlMatrix_.getValues(),
332  Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
333  // free temporary vectors
334  delete [] Flag_;
335  delete [] Pattern_;
336  delete [] Parent_;
337  delete [] Lnz_;
338 
339  if(rank!=dimMat)
340  DUNE_THROW(InvalidStateException,"Error: LDL factorisation failed!");
341  }
342 
343  LDLMatrix ldlMatrix_;
344  bool matrixIsLoaded_;
345  int verbose_;
346  int* Lp_;
347  int* Parent_;
348  int* Lnz_;
349  int* Flag_;
350  int* Pattern_;
351  int* P_;
352  int* Pinv_;
353  double* D_;
354  double* Y_;
355  double* Lx_;
356  int* Li_;
357  };
358 
359  template<typename T, typename A, int n, int m>
361  {
362  enum {value = true};
363  };
364 
365  template<typename T, typename A, int n, int m>
367  {
368  enum {value = true};
369  };
370 
371  struct LDLCreator {
372  template<class F> struct isValidBlock : std::false_type{};
373  template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
374 
375  template<typename TL, typename M>
376  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
377  typename Dune::TypeListElement<2, TL>::type>>
378  operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
379  std::enable_if_t<
380  isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
381  {
382  int verbose = config.get("verbose", 0);
383  return std::make_shared<Dune::LDL<M>>(mat,verbose);
384  }
385 
386  // second version with SFINAE to validate the template parameters of LDL
387  template<typename TL, typename M>
388  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
389  typename Dune::TypeListElement<2, TL>::type>>
390  operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
391  std::enable_if_t<
392  !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
393  {
394  DUNE_THROW(UnsupportedType,
395  "Unsupported Type in LDL (only double and std::complex<double> supported)");
396  }
397  };
399 
400 } // end namespace Dune
401 
402 
403 #endif //HAVE_SUITESPARSE_LDL
404 #endif //DUNE_ISTL_LDL_HH
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Dune::BlockVector< FieldVector< T, m >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, m > > > domain_type
The type of the domain of the solver.
Definition: ldl.hh:83
void setOption([[maybe_unused]] unsigned int option, [[maybe_unused]] double value)
Definition: ldl.hh:184
std::shared_ptr< Dune::InverseOperator< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL, const M &mat, const Dune::ParameterTree &config, std::enable_if_t< isValidBlock< typename Dune::TypeListElement< 1, TL >::type::block_type >::value, int >=0) const
Definition: ldl.hh:378
Dune::ISTL::Impl::BCCSMatrix< T, int > LDLMatrix
The corresponding SuperLU Matrix type.
Definition: ldl.hh:79
DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator())
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition: ldl.hh:226
void apply(T *x, T *b)
Additional apply method with c-arrays in analogy to superlu.
Definition: ldl.hh:174
virtual void apply(domain_type &x, range_type &b, [[maybe_unused]] double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: ldl.hh:164
double * getLx()
Get factorization Lx.
Definition: ldl.hh:294
Dune::BlockVector< FieldVector< T, n >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, n > > > range_type
The type of the range of the solver.
Definition: ldl.hh:85
virtual ~LDL()
Default constructor.
Definition: ldl.hh:143
void free()
Free allocated space.
Definition: ldl.hh:244
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: ldl.hh:150
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: ldl.hh:188
const char * name()
Get method name.
Definition: ldl.hh:258
int * getLi()
Get factorization Li.
Definition: ldl.hh:285
LDL(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: ldl.hh:118
LDL(const Matrix &matrix, const ParameterTree &config)
Constructs the LDL solver.
Definition: ldl.hh:134
int * getLp()
Get factorization Lp.
Definition: ldl.hh:276
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: ldl.hh:88
ISTL::Impl::BCCSMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A >, int > MatrixInitializer
Type of an associated initializer class.
Definition: ldl.hh:81
void setSubMatrix(const Matrix &matrix, const S &rowIndexSet)
Definition: ldl.hh:205
LDL()
Default constructor.
Definition: ldl.hh:139
double * getD()
Get factorization diagonal matrix D.
Definition: ldl.hh:267
LDLMatrix & getInternalMatrix()
Return the column compress matrix.
Definition: ldl.hh:235
LDL(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition: ldl.hh:102
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
Definition: matrixutils.hh:209
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
A vector of blocks with memory management.
Definition: bvector.hh:393
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
Definition: overlappingschwarz.hh:692
Use the LDL package to directly solve linear systems – empty default class.
Definition: ldl.hh:54
Definition: ldl.hh:371
Definition: ldl.hh:372
Definition: matrixutils.hh:25
Statistics about the application of an inverse operator.
Definition: solver.hh:46
int iterations
Number of iterations.
Definition: solver.hh:65
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Abstract base class for all solvers.
Definition: solver.hh:97
Category
Definition: solvercategory.hh:21
Definition: solverregistry.hh:75
Definition: solvertype.hh:14
@ value
Whether this is a direct solver.
Definition: solvertype.hh:22
Definition: solvertype.hh:28
@ value
whether the solver internally uses column compressed storage
Definition: solvertype.hh:34