dune-istl  2.9.0
cholmod.hh
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3 #pragma once
4 
5 #if HAVE_SUITESPARSE_CHOLMOD
6 
7 #include <dune/common/fmatrix.hh>
8 #include <dune/common/fvector.hh>
10 #include <dune/istl/bvector.hh>
11 #include<dune/istl/solver.hh>
13 #include <dune/istl/foreach.hh>
14 
15 #include <vector>
16 #include <memory>
17 
18 #include <cholmod.h>
19 
20 namespace Dune {
21 
22 namespace Impl{
23 
32  struct NoIgnore
33  {
34  const NoIgnore& operator[](std::size_t) const { return *this; }
35  explicit operator bool() const { return false; }
36  static constexpr std::size_t size() { return 0; }
37 
38  };
39 
40 
41  template<class BlockedVector, class FlatVector>
42  void copyToFlatVector(const BlockedVector& blockedVector, FlatVector& flatVector)
43  {
44  // traverse the vector once just to compute the size
45  std::size_t len = flatVectorForEach(blockedVector, [&](auto&&, auto...){});
46  flatVector.resize(len);
47 
48  flatVectorForEach(blockedVector, [&](auto&& entry, auto offset){
49  flatVector[offset] = entry;
50  });
51  }
52 
53  // special (dummy) case for NoIgnore
54  template<class FlatVector>
55  void copyToFlatVector(const NoIgnore&, FlatVector&)
56  {
57  // just do nothing
58  return;
59  }
60 
61  template<class FlatVector, class BlockedVector>
62  void copyToBlockedVector(const FlatVector& flatVector, BlockedVector& blockedVector)
63  {
64  flatVectorForEach(blockedVector, [&](auto& entry, auto offset){
65  entry = flatVector[offset];
66  });
67  }
68 
69 
70 } //namespace Impl
71 
76 template<class Vector>
77 class Cholmod : public InverseOperator<Vector, Vector>
78 {
79 public:
80 
86  Cholmod()
87  {
88  cholmod_start(&c_);
89  }
90 
96  ~Cholmod()
97  {
98  if (L_)
99  cholmod_free_factor(&L_, &c_);
100  cholmod_finish(&c_);
101  }
102 
103  // forbid copying to avoid freeing memory twice
104  Cholmod(const Cholmod&) = delete;
105  Cholmod& operator=(const Cholmod&) = delete;
106 
107 
110  void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
111  {
112  apply(x,b,res);
113  }
114 
120  void apply(Vector& x, Vector& b, InverseOperatorResult& res)
121  {
122  // do nothing if N=0
123  if ( nIsZero_ )
124  {
125  return;
126  }
127 
128  if (x.size() != b.size())
129  DUNE_THROW(Exception, "Error in apply(): sizes of x and b do not match!");
130 
131  // cast to double array
132  auto b2 = std::make_unique<double[]>(L_->n);
133  auto x2 = std::make_unique<double[]>(L_->n);
134 
135  // copy to cholmod
136  auto bp = b2.get();
137 
138  flatVectorForEach(b, [&](auto&& entry, auto&& flatIndex){
139  if ( subIndices_.empty() )
140  bp[ flatIndex ] = entry;
141  else
142  if( subIndices_[ flatIndex ] != std::numeric_limits<std::size_t>::max() )
143  bp[ subIndices_[ flatIndex ] ] = entry;
144  });
145 
146  // create a cholmod dense object
147  auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
148  // cast because void-ptr
149  auto b4 = static_cast<double*>(b3->x);
150  std::copy(b2.get(), b2.get() + L_->n, b4);
151 
152  // solve for a cholmod x object
153  auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
154  // cast because void-ptr
155  auto xp = static_cast<double*>(x3->x);
156 
157  // copy into x
158  flatVectorForEach(x, [&](auto&& entry, auto&& flatIndex){
159  if ( subIndices_.empty() )
160  entry = xp[ flatIndex ];
161  else
162  if( subIndices_[ flatIndex ] != std::numeric_limits<std::size_t>::max() )
163  entry = xp[ subIndices_[ flatIndex ] ];
164  });
165 
166  // statistics for a direct solver
167  res.iterations = 1;
168  res.converged = true;
169  }
170 
171 
177  template<class Matrix>
178  void setMatrix(const Matrix& matrix)
179  {
180  const Impl::NoIgnore* noIgnore = nullptr;
181  setMatrix(matrix, noIgnore);
182  }
183 
198  template<class Matrix, class Ignore>
199  void setMatrix(const Matrix& matrix, const Ignore* ignore)
200  {
201  // count the number of entries and diagonal entries
202  int nonZeros = 0;
203  int numberOfIgnoredDofs = 0;
204 
205 
206  auto [flatRows,flatCols] = flatMatrixForEach( matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
207  if( flatRowIndex <= flatColIndex )
208  nonZeros++;
209  });
210 
211  std::vector<bool> flatIgnore;
212 
213  if ( ignore )
214  {
215  Impl::copyToFlatVector(*ignore,flatIgnore);
216  numberOfIgnoredDofs = std::count(flatIgnore.begin(),flatIgnore.end(),true);
217  }
218 
219  // Total number of rows
220  int N = flatRows - numberOfIgnoredDofs;
221 
222  nIsZero_ = (N <= 0);
223 
224  if ( nIsZero_ )
225  {
226  return;
227  }
228 
229  /*
230  * CHOLMOD uses compressed-column sparse matrices, but for symmetric
231  * matrices this is the same as the compressed-row sparse matrix used
232  * by DUNE. So we can just store Mᵀ instead of M (as M = Mᵀ).
233  */
234  const auto deleter = [c = &this->c_](auto* p) {
235  cholmod_free_sparse(&p, c);
236  };
237  auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
238  cholmod_allocate_sparse(N, // # rows
239  N, // # cols
240  nonZeros, // # of nonzeroes
241  1, // indices are sorted ( 1 = true)
242  1, // matrix is "packed" ( 1 = true)
243  -1, // stype of matrix ( -1 = consider the lower part only )
244  CHOLMOD_REAL, // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
245  &c_ // cholmod_common ptr
246  ), deleter);
247 
248  // copy the data of BCRS matrix to Cholmod Sparse matrix
249  int* Ap = static_cast<int*>(M->p);
250  int* Ai = static_cast<int*>(M->i);
251  double* Ax = static_cast<double*>(M->x);
252 
253 
254  if ( ignore )
255  {
256  // init the mapping
257  subIndices_.resize(flatRows,std::numeric_limits<std::size_t>::max());
258 
259  std::size_t subIndexCounter = 0;
260 
261  for ( std::size_t i=0; i<flatRows; i++ )
262  {
263  if ( not flatIgnore[ i ] )
264  {
265  subIndices_[ i ] = subIndexCounter++;
266  }
267  }
268  }
269 
270  // at first, we need to compute the row starts "Ap"
271  // therefore, we count all (not ignored) entries in each row and in the end we accumulate everything
272  flatMatrixForEach(matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
273 
274  // stop if ignored
275  if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
276  return;
277 
278  // stop if in lower half
279  if ( flatRowIndex > flatColIndex )
280  return;
281 
282  // ok, count the entry
283  auto idx = ignore ? subIndices_[flatRowIndex] : flatRowIndex;
284  Ap[idx+1]++;
285 
286  });
287 
288  // now accumulate
289  Ap[0] = 0;
290  for ( int i=0; i<N; i++ )
291  {
292  Ap[i+1] += Ap[i];
293  }
294 
295  // we need a compressed row position counter
296  std::vector<std::size_t> rowPosition(N,0);
297 
298  // now we can set the entries
299  flatMatrixForEach(matrix, [&](auto&& entry, auto&& flatRowIndex, auto&& flatColIndex){
300 
301  // stop if ignored
302  if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
303  return;
304 
305  // stop if in lower half
306  if ( flatRowIndex > flatColIndex )
307  return;
308 
309  // ok, set the entry
310  auto rowIdx = ignore ? subIndices_[flatRowIndex] : flatRowIndex;
311  auto colIdx = ignore ? subIndices_[flatColIndex] : flatColIndex;
312  auto rowStart = Ap[rowIdx];
313  auto rowPos = rowPosition[rowIdx];
314  Ai[ rowStart + rowPos ] = colIdx;
315  Ax[ rowStart + rowPos ] = entry;
316  rowPosition[rowIdx]++;
317 
318  });
319 
320  // Now analyse the pattern and optimal row order
321  L_ = cholmod_analyze(M.get(), &c_);
322 
323  // Do the factorization (this may take some time)
324  cholmod_factorize(M.get(), L_, &c_);
325  }
326 
327  virtual SolverCategory::Category category() const
328  {
329  return SolverCategory::Category::sequential;
330  }
331 
337  cholmod_common& cholmodCommonObject()
338  {
339  return c_;
340  }
341 
347  cholmod_factor& cholmodFactor()
348  {
349  return *L_;
350  }
351 
357  const cholmod_factor& cholmodFactor() const
358  {
359  return *L_;
360  }
361 private:
362 
363  // create a std::unique_ptr to a cholmod_dense object with a deleter
364  // that calls the appropriate cholmod cleanup routine
365  auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
366  {
367  const auto deleter = [c](auto* p) {
368  cholmod_free_dense(&p, c);
369  };
370  return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
371  }
372 
373  cholmod_common c_;
374  cholmod_factor* L_ = nullptr;
375 
376  // indicator for a 0x0 problem (due to ignore dof's)
377  bool nIsZero_ = false;
378 
379  // vector mapping all indices in flat order to the not ignored indices
380  std::vector<std::size_t> subIndices_;
381 };
382 
383  struct CholmodCreator{
384  template<class F> struct isValidBlock : std::false_type{};
385  template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
386  template<int k> struct isValidBlock<FieldVector<float,k>> : std::true_type{};
387 
388  template<class TL, typename M>
389  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
390  typename Dune::TypeListElement<2, TL>::type>>
391  operator()(TL /*tl*/, const M& mat, const Dune::ParameterTree& /*config*/,
392  std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
393  {
394  using D = typename Dune::TypeListElement<1, TL>::type;
395  auto solver = std::make_shared<Dune::Cholmod<D>>();
396  solver->setMatrix(mat);
397  return solver;
398  }
399 
400  // second version with SFINAE to validate the template parameters of Cholmod
401  template<typename TL, typename M>
402  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
403  typename Dune::TypeListElement<2, TL>::type>>
404  operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
405  std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
406  {
407  DUNE_THROW(UnsupportedType, "Unsupported Type in Cholmod");
408  }
409  };
410  DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator());
411 
412 } /* namespace Dune */
413 
414 #endif // HAVE_SUITESPARSE_CHOLMOD
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Define general, extensible interface for inverse operators.
DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator())
Matrix & mat
Definition: matrixmatrix.hh:347
Definition: allocator.hh:11
std::size_t flatVectorForEach(Vector &&vector, F &&f, std::size_t offset=0)
Traverse a blocked vector and call a functor at each scalar entry.
Definition: foreach.hh:95
std::pair< std::size_t, std::size_t > flatMatrixForEach(Matrix &&matrix, F &&f, std::size_t rowOffset=0, std::size_t colOffset=0)
Traverse a blocked matrix and call a functor at each scalar entry.
Definition: foreach.hh:132
Category
Definition: solvercategory.hh:23