dune-istl  2.8.0
matrixutils.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_MATRIXUTILS_HH
4 #define DUNE_ISTL_MATRIXUTILS_HH
5 
6 #include <set>
7 #include <vector>
8 #include <limits>
9 #include <dune/common/typetraits.hh>
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/dynmatrix.hh>
12 #include <dune/common/diagonalmatrix.hh>
13 #include <dune/common/scalarmatrixview.hh>
15 #include "istlexception.hh"
16 
17 namespace Dune
18 {
19 
20 #ifndef DOYXGEN
21  template<typename B, typename A>
22  class BCRSMatrix;
23 
24  template<typename K, int n, int m>
25  class FieldMatrix;
26 
27  template<class T, class A>
28  class Matrix;
29 #endif
30 
44  template<class Matrix, std::size_t blocklevel, std::size_t l=blocklevel>
46  {
51  static void check([[maybe_unused]] const Matrix& mat)
52  {
53 #ifdef DUNE_ISTL_WITH_CHECKING
54  typedef typename Matrix::ConstRowIterator Row;
55  typedef typename Matrix::ConstColIterator Entry;
56  for(Row row = mat.begin(); row!=mat.end(); ++row) {
57  Entry diagonal = row->find(row.index());
58  if(diagonal==row->end())
59  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
60  <<" at block recursion level "<<l-blocklevel);
61  else{
62  auto m = Impl::asMatrix(*diagonal);
63  CheckIfDiagonalPresent<decltype(m),blocklevel-1,l>::check(m);
64  }
65  }
66 #endif
67  }
68  };
69 
70  template<class Matrix, std::size_t l>
72  {
73  static void check(const Matrix& mat)
74  {
75  typedef typename Matrix::ConstRowIterator Row;
76  for(Row row = mat.begin(); row!=mat.end(); ++row) {
77  if(row->find(row.index())==row->end())
78  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
79  <<" at block recursion level "<<l);
80  }
81  }
82  };
83 
84  template<typename FirstRow, typename... Args>
85  class MultiTypeBlockMatrix;
86 
87  template<std::size_t blocklevel, std::size_t l, typename T1, typename... Args>
89  blocklevel,l>
90  {
91  typedef MultiTypeBlockMatrix<T1,Args...> Matrix;
92 
97  static void check(const Matrix& /* mat */)
98  {
99 #ifdef DUNE_ISTL_WITH_CHECKING
100  // TODO Implement check
101 #endif
102  }
103  };
104 
116  template<class M>
117  inline auto countNonZeros(const M&,
118  [[maybe_unused]] typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae = nullptr)
119  {
120  return 1;
121  }
122 
123  template<class M>
124  inline auto countNonZeros(const M& matrix,
125  [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
126  {
127  typename M::size_type nonZeros = 0;
128  for(auto&& row : matrix)
129  for(auto&& entry : row)
130  nonZeros += countNonZeros(entry);
131  return nonZeros;
132  }
133 
134  /*
135  template<class M>
136  struct ProcessOnFieldsOfMatrix
137  */
138 
140  namespace
141  {
142  struct CompPair {
143  template<class G,class M>
144  bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2) const
145  {
146  return p1.first<p2.first;
147  }
148  };
149 
150  }
151  template<class M, class C>
152  void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os)
153  {
154  typedef typename C::ParallelIndexSet::const_iterator IIter;
155  typedef typename C::OwnerSet OwnerSet;
156  typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
157 
158  GlobalIndex gmax=0;
159 
160  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
161  idx!=eidx; ++idx)
162  gmax=std::max(gmax,idx->global());
163 
164  gmax=ooc.communicator().max(gmax);
165  ooc.buildGlobalLookup();
166 
167  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
168  idx!=eidx; ++idx) {
169  if(OwnerSet::contains(idx->local().attribute()))
170  {
171  typedef typename M::block_type Block;
172 
173  std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
174 
175  // sort rows
176  typedef typename M::ConstColIterator CIter;
177  for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
178  c!=cend; ++c) {
179  const typename C::ParallelIndexSet::IndexPair* pair
180  =ooc.globalLookup().pair(c.index());
181  assert(pair);
182  entries.insert(std::make_pair(pair->global(), *c));
183  }
184 
185  //wait until its the rows turn.
186  GlobalIndex rowidx = idx->global();
187  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
188  while(cur!=rowidx)
189  cur=ooc.communicator().min(rowidx);
190 
191  // print rows
192  typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
193  for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
194  os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl;
195 
196 
197  }
198  }
199 
200  ooc.freeGlobalLookup();
201  // Wait until everybody is finished
202  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
203  while(cur!=ooc.communicator().min(cur)) ;
204  }
205 
206  // Default implementation for scalar types
207  template<typename M>
209  {
210  static_assert(IsNumber<M>::value, "MatrixDimension is not implemented for this type!");
211 
212  static auto rowdim(const M& A)
213  {
214  return 1;
215  }
216 
217  static auto coldim(const M& A)
218  {
219  return 1;
220  }
221  };
222 
223  // Default implementation for scalar types
224  template<typename B, typename TA>
225  struct MatrixDimension<Matrix<B,TA> >
226  {
229 
230  static size_type rowdim (const Matrix<B,TA>& A, size_type i)
231  {
232  return MatrixDimension<block_type>::rowdim(A[i][0]);
233  }
234 
235  static size_type coldim (const Matrix<B,TA>& A, size_type c)
236  {
237  return MatrixDimension<block_type>::coldim(A[0][c]);
238  }
239 
240  static size_type rowdim (const Matrix<B,TA>& A)
241  {
242  size_type nn=0;
243  for (size_type i=0; i<A.N(); i++)
244  nn += rowdim(A,i);
245  return nn;
246  }
247 
248  static size_type coldim (const Matrix<B,TA>& A)
249  {
250  size_type nn=0;
251  for (size_type i=0; i<A.M(); i++)
252  nn += coldim(A,i);
253  return nn;
254  }
255  };
256 
257 
258  template<typename B, typename TA>
260  {
262  typedef typename Matrix::block_type block_type;
263  typedef typename Matrix::size_type size_type;
264 
265  static size_type rowdim (const Matrix& A, size_type i)
266  {
267  const B* row = A.r[i].getptr();
268  if(row)
270  else
271  return 0;
272  }
273 
274  static size_type coldim (const Matrix& A, size_type c)
275  {
276  // find an entry in column c
277  if (A.nnz_ > 0)
278  {
279  for (size_type k=0; k<A.nnz_; k++) {
280  if (A.j_.get()[k] == c) {
282  }
283  }
284  }
285  else
286  {
287  for (size_type i=0; i<A.N(); i++)
288  {
289  size_type* j = A.r[i].getindexptr();
290  B* a = A.r[i].getptr();
291  for (size_type k=0; k<A.r[i].getsize(); k++)
292  if (j[k]==c) {
294  }
295  }
296  }
297 
298  // not found
299  return 0;
300  }
301 
302  static size_type rowdim (const Matrix& A){
303  size_type nn=0;
304  for (size_type i=0; i<A.N(); i++)
305  nn += rowdim(A,i);
306  return nn;
307  }
308 
309  static size_type coldim (const Matrix& A){
310  typedef typename Matrix::ConstRowIterator ConstRowIterator;
311  typedef typename Matrix::ConstColIterator ConstColIterator;
312 
313  // The following code has a complexity of nnz, and
314  // typically a very small constant.
315  //
316  std::vector<size_type> coldims(A.M(),
317  std::numeric_limits<size_type>::max());
318 
319  for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
320  for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
321  // only compute blocksizes we don't already have
322  if (coldims[col.index()]==std::numeric_limits<size_type>::max())
323  coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
324 
325  size_type sum = 0;
326  for (typename std::vector<size_type>::iterator it=coldims.begin();
327  it!=coldims.end(); ++it)
328  // skip rows for which no coldim could be determined
329  if ((*it)>=0)
330  sum += *it;
331 
332  return sum;
333  }
334  };
335 
336 
337  template<typename B, int n, int m, typename TA>
339  {
341  typedef typename Matrix::size_type size_type;
342 
343  static size_type rowdim (const Matrix& /*A*/, size_type /*i*/)
344  {
345  return n;
346  }
347 
348  static size_type coldim (const Matrix& /*A*/, size_type /*c*/)
349  {
350  return m;
351  }
352 
353  static size_type rowdim (const Matrix& A) {
354  return A.N()*n;
355  }
356 
357  static size_type coldim (const Matrix& A) {
358  return A.M()*m;
359  }
360  };
361 
362  template<typename K, int n, int m>
363  struct MatrixDimension<FieldMatrix<K,n,m> >
364  {
366  typedef typename Matrix::size_type size_type;
367 
368  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
369  {
370  return 1;
371  }
372 
373  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
374  {
375  return 1;
376  }
377 
378  static size_type rowdim(const Matrix& /*A*/)
379  {
380  return n;
381  }
382 
383  static size_type coldim(const Matrix& /*A*/)
384  {
385  return m;
386  }
387  };
388 
389  template <class T>
390  struct MatrixDimension<Dune::DynamicMatrix<T> >
391  {
392  typedef Dune::DynamicMatrix<T> MatrixType;
393  typedef typename MatrixType::size_type size_type;
394 
395  static size_type rowdim(const MatrixType& /*A*/, size_type /*r*/)
396  {
397  return 1;
398  }
399 
400  static size_type coldim(const MatrixType& /*A*/, size_type /*r*/)
401  {
402  return 1;
403  }
404 
405  static size_type rowdim(const MatrixType& A)
406  {
407  return A.N();
408  }
409 
410  static size_type coldim(const MatrixType& A)
411  {
412  return A.M();
413  }
414  };
415 
416  template<typename K, int n, int m, typename TA>
417  struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
418  {
421 
422  static size_type rowdim(const ThisMatrix& /*A*/, size_type /*r*/)
423  {
424  return n;
425  }
426 
427  static size_type coldim(const ThisMatrix& /*A*/, size_type /*r*/)
428  {
429  return m;
430  }
431 
432  static size_type rowdim(const ThisMatrix& A)
433  {
434  return A.N()*n;
435  }
436 
437  static size_type coldim(const ThisMatrix& A)
438  {
439  return A.M()*m;
440  }
441  };
442 
443  template<typename K, int n>
444  struct MatrixDimension<DiagonalMatrix<K,n> >
445  {
446  typedef DiagonalMatrix<K,n> Matrix;
447  typedef typename Matrix::size_type size_type;
448 
449  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
450  {
451  return 1;
452  }
453 
454  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
455  {
456  return 1;
457  }
458 
459  static size_type rowdim(const Matrix& /*A*/)
460  {
461  return n;
462  }
463 
464  static size_type coldim(const Matrix& /*A*/)
465  {
466  return n;
467  }
468  };
469 
470  template<typename K, int n>
472  {
474  typedef typename Matrix::size_type size_type;
475 
476  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
477  {
478  return 1;
479  }
480 
481  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
482  {
483  return 1;
484  }
485 
486  static size_type rowdim(const Matrix& /*A*/)
487  {
488  return n;
489  }
490 
491  static size_type coldim(const Matrix& /*A*/)
492  {
493  return n;
494  }
495  };
496 
500  template<typename T>
501  struct IsMatrix
502  {
503  enum {
507  value = false
508  };
509  };
510 
511  template<typename T>
512  struct IsMatrix<DenseMatrix<T> >
513  {
514  enum {
518  value = true
519  };
520  };
521 
522 
523  template<typename T, typename A>
524  struct IsMatrix<BCRSMatrix<T,A> >
525  {
526  enum {
530  value = true
531  };
532  };
533 
534  template<typename T>
536  {
537  bool operator()(const T* l, const T* r)
538  {
539  return *l < *r;
540  }
541  };
542 
543 }
544 #endif
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
auto countNonZeros(const M &, [[maybe_unused]] typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition: matrixutils.hh:152
Definition: matrixutils.hh:209
static auto coldim(const M &A)
Definition: matrixutils.hh:217
static auto rowdim(const M &A)
Definition: matrixutils.hh:212
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:498
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:739
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:489
Iterator access to matrix rows
Definition: bcrsmatrix.hh:577
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
derive error class from the base class in common
Definition: istlexception.hh:17
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:618
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:612
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
Definition: matrixutils.hh:25
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:46
static void check([[maybe_unused]] const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:51
static void check(const Matrix &mat)
Definition: matrixutils.hh:73
static void check(const Matrix &)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:97
MultiTypeBlockMatrix< T1, Args... > Matrix
Definition: matrixutils.hh:91
static size_type rowdim(const Matrix< B, TA > &A, size_type i)
Definition: matrixutils.hh:230
static size_type coldim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:248
static size_type rowdim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:240
typename Matrix< B, TA >::size_type size_type
Definition: matrixutils.hh:228
static size_type coldim(const Matrix< B, TA > &A, size_type c)
Definition: matrixutils.hh:235
typename Matrix< B, TA >::block_type block_type
Definition: matrixutils.hh:227
BCRSMatrix< B, TA > Matrix
Definition: matrixutils.hh:261
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:309
Matrix::block_type block_type
Definition: matrixutils.hh:262
static size_type coldim(const Matrix &A, size_type c)
Definition: matrixutils.hh:274
Matrix::size_type size_type
Definition: matrixutils.hh:263
static size_type rowdim(const Matrix &A, size_type i)
Definition: matrixutils.hh:265
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:302
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:357
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:343
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:353
Matrix::size_type size_type
Definition: matrixutils.hh:341
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition: matrixutils.hh:340
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:348
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:368
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:383
Matrix::size_type size_type
Definition: matrixutils.hh:366
FieldMatrix< K, n, m > Matrix
Definition: matrixutils.hh:365
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:373
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:378
static size_type coldim(const MatrixType &A)
Definition: matrixutils.hh:410
static size_type rowdim(const MatrixType &A)
Definition: matrixutils.hh:405
static size_type rowdim(const MatrixType &, size_type)
Definition: matrixutils.hh:395
MatrixType::size_type size_type
Definition: matrixutils.hh:393
static size_type coldim(const MatrixType &, size_type)
Definition: matrixutils.hh:400
Dune::DynamicMatrix< T > MatrixType
Definition: matrixutils.hh:392
static size_type coldim(const ThisMatrix &A)
Definition: matrixutils.hh:437
static size_type rowdim(const ThisMatrix &A)
Definition: matrixutils.hh:432
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition: matrixutils.hh:419
static size_type coldim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:427
static size_type rowdim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:422
ThisMatrix::size_type size_type
Definition: matrixutils.hh:420
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:454
Matrix::size_type size_type
Definition: matrixutils.hh:447
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:464
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:459
DiagonalMatrix< K, n > Matrix
Definition: matrixutils.hh:446
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:449
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:491
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:476
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:481
Matrix::size_type size_type
Definition: matrixutils.hh:474
ScaledIdentityMatrix< K, n > Matrix
Definition: matrixutils.hh:473
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:486
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:507
Definition: matrixutils.hh:536
bool operator()(const T *l, const T *r)
Definition: matrixutils.hh:537
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:28
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41