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