dune-istl  2.5.1
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/unused.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 
40  namespace
41  {
42 
43  template<int i>
44  struct NonZeroCounter
45  {
46  template<class M>
47  static typename M::size_type count(const M& matrix)
48  {
49  typedef typename M::ConstRowIterator RowIterator;
50 
51  RowIterator endRow = matrix.end();
52  typename M::size_type nonZeros = 0;
53 
54  for(RowIterator row = matrix.begin(); row != endRow; ++row) {
55  typedef typename M::ConstColIterator Entry;
56  Entry endEntry = row->end();
57  for(Entry entry = row->begin(); entry != endEntry; ++entry) {
58  nonZeros += NonZeroCounter<i-1>::count(*entry);
59  }
60  }
61  return nonZeros;
62  }
63  };
64 
65  template<>
66  struct NonZeroCounter<1>
67  {
68  template<class M>
69  static typename M::size_type count(const M& matrix)
70  {
71  return matrix.N()*matrix.M();
72  }
73  };
74 
75  }
76 
81  template<class Matrix, std::size_t blocklevel, std::size_t l=blocklevel>
83  {
88  static void check(const Matrix& mat)
89  {
90  DUNE_UNUSED_PARAMETER(mat);
91 #ifdef DUNE_ISTL_WITH_CHECKING
92  typedef typename Matrix::ConstRowIterator Row;
93  typedef typename Matrix::ConstColIterator Entry;
94  for(Row row = mat.begin(); row!=mat.end(); ++row) {
95  Entry diagonal = row->find(row.index());
96  if(diagonal==row->end())
97  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
98  <<" at block recursion level "<<l-blocklevel);
99  else
101  }
102 #endif
103  }
104  };
105 
106  template<class Matrix, std::size_t l>
108  {
109  static void check(const Matrix& mat)
110  {
111  typedef typename Matrix::ConstRowIterator Row;
112  for(Row row = mat.begin(); row!=mat.end(); ++row) {
113  if(row->find(row.index())==row->end())
114  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
115  <<" at block recursion level "<<l);
116  }
117  }
118  };
119 
120  template<typename FirstRow, typename... Args>
122 
123  template<std::size_t blocklevel, std::size_t l, typename T1, typename... Args>
125  blocklevel,l>
126  {
127  typedef MultiTypeBlockMatrix<T1,Args...> Matrix;
128 
133  static void check(const Matrix& /* mat */)
134  {
135 #ifdef DUNE_ISTL_WITH_CHECKING
136  // TODO Implement check
137 #endif
138  }
139  };
140 
152  template<class M>
153  inline int countNonZeros(const M& matrix)
154  {
155  return NonZeroCounter<M::blocklevel>::count(matrix);
156  }
157  /*
158  template<class M>
159  struct ProcessOnFieldsOfMatrix
160  */
161 
163  namespace
164  {
165  struct CompPair {
166  template<class G,class M>
167  bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2)
168  {
169  return p1.first<p2.first;
170  }
171  };
172 
173  }
174  template<class M, class C>
175  void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os)
176  {
177  typedef typename C::ParallelIndexSet::const_iterator IIter;
178  typedef typename C::OwnerSet OwnerSet;
179  typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
180 
181  GlobalIndex gmax=0;
182 
183  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
184  idx!=eidx; ++idx)
185  gmax=std::max(gmax,idx->global());
186 
187  gmax=ooc.communicator().max(gmax);
188  ooc.buildGlobalLookup();
189 
190  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
191  idx!=eidx; ++idx) {
192  if(OwnerSet::contains(idx->local().attribute()))
193  {
194  typedef typename M::block_type Block;
195 
196  std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
197 
198  // sort rows
199  typedef typename M::ConstColIterator CIter;
200  for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
201  c!=cend; ++c) {
202  const typename C::ParallelIndexSet::IndexPair* pair
203  =ooc.globalLookup().pair(c.index());
204  assert(pair);
205  entries.insert(std::make_pair(pair->global(), *c));
206  }
207 
208  //wait until its the rows turn.
209  GlobalIndex rowidx = idx->global();
210  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
211  while(cur!=rowidx)
212  cur=ooc.communicator().min(rowidx);
213 
214  // print rows
215  typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
216  for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
217  os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl;
218 
219 
220  }
221  }
222 
223  ooc.freeGlobalLookup();
224  // Wait until everybody is finished
225  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
226  while(cur!=ooc.communicator().min(cur)) ;
227  }
228 
229  template<typename M>
230  struct MatrixDimension
231  {};
232 
233 
234  template<typename B, typename TA>
236  {
238  typedef typename Matrix::block_type block_type;
239  typedef typename Matrix::size_type size_type;
240 
241  static size_type rowdim (const Matrix& A, size_type i)
242  {
243  const B* row = A.r[i].getptr();
244  if(row)
246  else
247  return 0;
248  }
249 
250  static size_type coldim (const Matrix& A, size_type c)
251  {
252  // find an entry in column c
253  if (A.nnz_ > 0)
254  {
255  for (size_type k=0; k<A.nnz_; k++) {
256  if (A.j_.get()[k] == c) {
258  }
259  }
260  }
261  else
262  {
263  for (size_type i=0; i<A.N(); i++)
264  {
265  size_type* j = A.r[i].getindexptr();
266  B* a = A.r[i].getptr();
267  for (size_type k=0; k<A.r[i].getsize(); k++)
268  if (j[k]==c) {
270  }
271  }
272  }
273 
274  // not found
275  return 0;
276  }
277 
278  static size_type rowdim (const Matrix& A){
279  size_type nn=0;
280  for (size_type i=0; i<A.N(); i++)
281  nn += rowdim(A,i);
282  return nn;
283  }
284 
285  static size_type coldim (const Matrix& A){
286  typedef typename Matrix::ConstRowIterator ConstRowIterator;
287  typedef typename Matrix::ConstColIterator ConstColIterator;
288 
289  // The following code has a complexity of nnz, and
290  // typically a very small constant.
291  //
292  std::vector<size_type> coldims(A.M(),
293  std::numeric_limits<size_type>::max());
294 
295  for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
296  for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
297  // only compute blocksizes we don't already have
298  if (coldims[col.index()]==std::numeric_limits<size_type>::max())
299  coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
300 
301  size_type sum = 0;
302  for (typename std::vector<size_type>::iterator it=coldims.begin();
303  it!=coldims.end(); ++it)
304  // skip rows for which no coldim could be determined
305  if ((*it)>=0)
306  sum += *it;
307 
308  return sum;
309  }
310  };
311 
312 
313  template<typename B, int n, int m, typename TA>
315  {
317  typedef typename Matrix::size_type size_type;
318 
319  static size_type rowdim (const Matrix& /*A*/, size_type /*i*/)
320  {
321  return n;
322  }
323 
324  static size_type coldim (const Matrix& /*A*/, size_type /*c*/)
325  {
326  return m;
327  }
328 
329  static size_type rowdim (const Matrix& A) {
330  return A.N()*n;
331  }
332 
333  static size_type coldim (const Matrix& A) {
334  return A.M()*m;
335  }
336  };
337 
338  template<typename K, int n, int m>
339  struct MatrixDimension<FieldMatrix<K,n,m> >
340  {
342  typedef typename Matrix::size_type size_type;
343 
344  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
345  {
346  return 1;
347  }
348 
349  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
350  {
351  return 1;
352  }
353 
354  static size_type rowdim(const Matrix& /*A*/)
355  {
356  return n;
357  }
358 
359  static size_type coldim(const Matrix& /*A*/)
360  {
361  return m;
362  }
363  };
364 
365  template <class T>
366  struct MatrixDimension<Dune::DynamicMatrix<T> >
367  {
368  typedef Dune::DynamicMatrix<T> MatrixType;
369  typedef typename MatrixType::size_type size_type;
370 
371  static size_type rowdim(const MatrixType& /*A*/, size_type /*r*/)
372  {
373  return 1;
374  }
375 
376  static size_type coldim(const MatrixType& /*A*/, size_type /*r*/)
377  {
378  return 1;
379  }
380 
381  static size_type rowdim(const MatrixType& A)
382  {
383  return A.N();
384  }
385 
386  static size_type coldim(const MatrixType& A)
387  {
388  return A.M();
389  }
390  };
391 
392  template<typename K, int n, int m, typename TA>
393  struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
394  {
397 
398  static size_type rowdim(const ThisMatrix& /*A*/, size_type /*r*/)
399  {
400  return n;
401  }
402 
403  static size_type coldim(const ThisMatrix& /*A*/, size_type /*r*/)
404  {
405  return m;
406  }
407 
408  static size_type rowdim(const ThisMatrix& A)
409  {
410  return A.N()*n;
411  }
412 
413  static size_type coldim(const ThisMatrix& A)
414  {
415  return A.M()*m;
416  }
417  };
418 
419  template<typename K, int n>
420  struct MatrixDimension<DiagonalMatrix<K,n> >
421  {
422  typedef DiagonalMatrix<K,n> Matrix;
423  typedef typename Matrix::size_type size_type;
424 
425  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
426  {
427  return 1;
428  }
429 
430  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
431  {
432  return 1;
433  }
434 
435  static size_type rowdim(const Matrix& /*A*/)
436  {
437  return n;
438  }
439 
440  static size_type coldim(const Matrix& /*A*/)
441  {
442  return n;
443  }
444  };
445 
446  template<typename K, int n>
448  {
450  typedef typename Matrix::size_type size_type;
451 
452  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
453  {
454  return 1;
455  }
456 
457  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
458  {
459  return 1;
460  }
461 
462  static size_type rowdim(const Matrix& /*A*/)
463  {
464  return n;
465  }
466 
467  static size_type coldim(const Matrix& /*A*/)
468  {
469  return n;
470  }
471  };
472 
476  template<typename T>
477  struct IsMatrix
478  {
479  enum {
483  value = false
484  };
485  };
486 
487  template<typename T>
488  struct IsMatrix<DenseMatrix<T> >
489  {
490  enum {
494  value = true
495  };
496  };
497 
498 
499  template<typename T, typename A>
500  struct IsMatrix<BCRSMatrix<T,A> >
501  {
502  enum {
506  value = true
507  };
508  };
509 
510  template<typename T>
512  {
513  bool operator()(const T* l, const T* r)
514  {
515  return *l < *r;
516  }
517  };
518 
519 }
520 #endif
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:329
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:88
static size_type coldim(const MatrixType &, size_type)
Definition: matrixutils.hh:376
static size_type coldim(const MatrixType &A)
Definition: matrixutils.hh:386
A Matrix class to support different block types.
Definition: matrixutils.hh:121
Col col
Definition: matrixmatrix.hh:349
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:435
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:153
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:462
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:571
static size_type rowdim(const ThisMatrix &A)
Definition: matrixutils.hh:408
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
A generic dynamic dense matrix.
Definition: matrix.hh:554
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:27
static size_type rowdim(const MatrixType &, size_type)
Definition: matrixutils.hh:371
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:640
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition: matrixutils.hh:316
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Definition: matrixutils.hh:511
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:634
ScaledIdentityMatrix< K, n > Matrix
Definition: matrixutils.hh:449
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:477
Definition: bcrsmatrix.hh:70
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:349
static void check(const Matrix &mat)
Definition: matrixutils.hh:109
MatrixType::size_type size_type
Definition: matrixutils.hh:369
Matrix::size_type size_type
Definition: matrixutils.hh:423
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:359
size_type M() const
Return the number of columns.
Definition: matrix.hh:695
static void check(const Matrix &)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:133
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:82
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition: matrixutils.hh:175
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:467
Definition: matrixutils.hh:25
B * a
Definition: bcrsmatrix.hh:1968
size_type N() const
Return the number of rows.
Definition: matrix.hh:690
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:430
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:333
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:319
Matrix::block_type block_type
Definition: matrixutils.hh:238
size_type getsize() const
get size
Definition: bvector.hh:1153
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41
static size_type coldim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:403
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:448
bool operator()(const T *l, const T *r)
Definition: matrixutils.hh:513
Matrix & mat
Definition: matrixmatrix.hh:345
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:440
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:452
Definition: basearray.hh:19
Matrix::size_type size_type
Definition: matrixutils.hh:239
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:425
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition: matrixutils.hh:395
size_type * getindexptr()
get pointer
Definition: bvector.hh:1136
static size_type coldim(const ThisMatrix &A)
Definition: matrixutils.hh:413
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:354
static size_type rowdim(const MatrixType &A)
Definition: matrixutils.hh:381
BCRSMatrix< B, TA > Matrix
Definition: matrixutils.hh:237
B * getptr()
get pointer
Definition: bvector.hh:1130
ConstIterator class for sequential access.
Definition: matrix.hh:397
std::shared_ptr< size_type > j_
Definition: bcrsmatrix.hh:1971
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1903
static size_type rowdim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:398
row_type * r
Definition: bcrsmatrix.hh:1965
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:344
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:457
static size_type coldim(const Matrix &A, size_type c)
Definition: matrixutils.hh:250
MultiTypeBlockMatrix< T1, Args... > Matrix
Definition: matrixutils.hh:127
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:324
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1909
DiagonalMatrix< K, n > Matrix
Definition: matrixutils.hh:422
static size_type rowdim(const Matrix &A, size_type i)
Definition: matrixutils.hh:241
Matrix::size_type size_type
Definition: matrixutils.hh:450
iterator class for sequential access
Definition: basearray.hh:563
derive error class from the base class in common
Definition: istlexception.hh:16
Iterator access to matrix rows
Definition: bcrsmatrix.hh:536
Matrix::size_type size_type
Definition: matrixutils.hh:342
ThisMatrix::size_type size_type
Definition: matrixutils.hh:396
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:278
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
size_type nnz_
Definition: bcrsmatrix.hh:1960
FieldMatrix< K, n, m > Matrix
Definition: matrixutils.hh:341
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:285
Dune::DynamicMatrix< T > MatrixType
Definition: matrixutils.hh:368
Matrix::size_type size_type
Definition: matrixutils.hh:317
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583