dune-istl  2.3.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_MATRIX_UTILS_HH
4 #define DUNE_MATRIX_UTILS_HH
5 
6 #include <set>
7 #include <vector>
8 #include <limits>
9 #include <dune/common/typetraits.hh>
10 #include <dune/common/static_assert.hh>
11 #include <dune/common/fmatrix.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 T1, typename T2, typename T3, typename T4, typename T5,
121  typename T6, typename T7, typename T8, typename T9>
123 
124  template<typename T1, typename T2, typename T3, typename T4, typename T5,
125  typename T6, typename T7, typename T8, typename T9, std::size_t blocklevel, std::size_t l>
126  struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>,
127  blocklevel,l>
128  {
130 
135  static void check(const Matrix& mat)
136  {
137 #ifdef DUNE_ISTL_WITH_CHECKING
138  // TODO Implement check
139 #endif
140  }
141  };
142 
154  template<class M>
155  inline int countNonZeros(const M& matrix)
156  {
158  }
159  /*
160  template<class M>
161  struct ProcessOnFieldsOfMatrix
162  */
163 
165  namespace
166  {
167  struct CompPair {
168  template<class G,class M>
169  bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2)
170  {
171  return p1.first<p2.first;
172  }
173  };
174 
175  }
176  template<class M, class C>
177  void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os)
178  {
179  typedef typename C::ParallelIndexSet::const_iterator IIter;
180  typedef typename C::OwnerSet OwnerSet;
181  typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
182 
183  GlobalIndex gmax=0;
184 
185  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
186  idx!=eidx; ++idx)
187  gmax=std::max(gmax,idx->global());
188 
189  gmax=ooc.communicator().max(gmax);
190  ooc.buildGlobalLookup();
191 
192  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
193  idx!=eidx; ++idx) {
194  if(OwnerSet::contains(idx->local().attribute()))
195  {
196  typedef typename M::block_type Block;
197 
198  std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
199 
200  // sort rows
201  typedef typename M::ConstColIterator CIter;
202  for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
203  c!=cend; ++c) {
204  const typename C::ParallelIndexSet::IndexPair* pair
205  =ooc.globalLookup().pair(c.index());
206  assert(pair);
207  entries.insert(std::make_pair(pair->global(), *c));
208  }
209 
210  //wait until its the rows turn.
211  GlobalIndex rowidx = idx->global();
212  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
213  while(cur!=rowidx)
214  cur=ooc.communicator().min(rowidx);
215 
216  // print rows
217  typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
218  for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
219  os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl;
220 
221 
222  }
223  }
224 
225  ooc.freeGlobalLookup();
226  // Wait until everybody is finished
227  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
228  while(cur!=ooc.communicator().min(cur)) ;
229  }
230 
231  template<typename M>
232  struct MatrixDimension
233  {};
234 
235 
236  template<typename B, typename TA>
238  {
240  typedef typename Matrix::block_type block_type;
241  typedef typename Matrix::size_type size_type;
242 
243  static size_type rowdim (const Matrix& A, size_type i)
244  {
245  const B* row = A.r[i].getptr();
246  if(row)
248  else
249  return 0;
250  }
251 
252  static size_type coldim (const Matrix& A, size_type c)
253  {
254  // find an entry in column c
255  if (A.nnz>0)
256  {
257  for (size_type k=0; k<A.nnz; k++) {
258  if (A.j.get()[k]==c) {
260  }
261  }
262  }
263  else
264  {
265  for (size_type i=0; i<A.N(); i++)
266  {
267  size_type* j = A.r[i].getindexptr();
268  B* a = A.r[i].getptr();
269  for (size_type k=0; k<A.r[i].getsize(); k++)
270  if (j[k]==c) {
272  }
273  }
274  }
275 
276  // not found
277  return 0;
278  }
279 
280  static size_type rowdim (const Matrix& A){
281  size_type nn=0;
282  for (size_type i=0; i<A.N(); i++)
283  nn += rowdim(A,i);
284  return nn;
285  }
286 
287  static size_type coldim (const Matrix& A){
288  typedef typename Matrix::ConstRowIterator ConstRowIterator;
289  typedef typename Matrix::ConstColIterator ConstColIterator;
290 
291  // The following code has a complexity of nnz, and
292  // typically a very small constant.
293  //
294  std::vector<size_type> coldims(A.M(),
295  std::numeric_limits<size_type>::max());
296 
297  for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
298  for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
299  // only compute blocksizes we don't already have
300  if (coldims[col.index()]==std::numeric_limits<size_type>::max())
301  coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
302 
303  size_type sum = 0;
304  for (typename std::vector<size_type>::iterator it=coldims.begin();
305  it!=coldims.end(); ++it)
306  // skip rows for which no coldim could be determined
307  if ((*it)>=0)
308  sum += *it;
309 
310  return sum;
311  }
312  };
313 
314 
315  template<typename B, int n, int m, typename TA>
317  {
319  typedef typename Matrix::size_type size_type;
320 
321  static size_type rowdim (const Matrix& /*A*/, size_type /*i*/)
322  {
323  return n;
324  }
325 
326  static size_type coldim (const Matrix& /*A*/, size_type /*c*/)
327  {
328  return m;
329  }
330 
331  static size_type rowdim (const Matrix& A) {
332  return A.N()*n;
333  }
334 
335  static size_type coldim (const Matrix& A) {
336  return A.M()*m;
337  }
338  };
339 
340  template<typename K, int n, int m>
341  struct MatrixDimension<FieldMatrix<K,n,m> >
342  {
344  typedef typename Matrix::size_type size_type;
345 
346  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
347  {
348  return 1;
349  }
350 
351  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
352  {
353  return 1;
354  }
355 
356  static size_type rowdim(const Matrix& /*A*/)
357  {
358  return n;
359  }
360 
361  static size_type coldim(const Matrix& /*A*/)
362  {
363  return m;
364  }
365  };
366 
367  template<typename K, int n, int m, typename TA>
368  struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
369  {
372 
373  static size_type rowdim(const ThisMatrix& /*A*/, size_type /*r*/)
374  {
375  return n;
376  }
377 
378  static size_type coldim(const ThisMatrix& /*A*/, size_type /*r*/)
379  {
380  return m;
381  }
382 
383  static size_type rowdim(const ThisMatrix& A)
384  {
385  return A.N()*n;
386  }
387 
388  static size_type coldim(const ThisMatrix& A)
389  {
390  return A.M()*m;
391  }
392  };
393 
394  template<typename K, int n>
395  struct MatrixDimension<DiagonalMatrix<K,n> >
396  {
397  typedef DiagonalMatrix<K,n> Matrix;
398  typedef typename Matrix::size_type size_type;
399 
400  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
401  {
402  return 1;
403  }
404 
405  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
406  {
407  return 1;
408  }
409 
410  static size_type rowdim(const Matrix& /*A*/)
411  {
412  return n;
413  }
414 
415  static size_type coldim(const Matrix& /*A*/)
416  {
417  return n;
418  }
419  };
420 
421  template<typename K, int n>
423  {
425  typedef typename Matrix::size_type size_type;
426 
427  static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
428  {
429  return 1;
430  }
431 
432  static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
433  {
434  return 1;
435  }
436 
437  static size_type rowdim(const Matrix& /*A*/)
438  {
439  return n;
440  }
441 
442  static size_type coldim(const Matrix& /*A*/)
443  {
444  return n;
445  }
446  };
447 
451  template<typename T>
452  struct IsMatrix
453  {
454  enum {
458  value = false
459  };
460  };
461 
462  template<typename T>
463  struct IsMatrix<DenseMatrix<T> >
464  {
465  enum {
469  value = true
470  };
471  };
472 
473 
474  template<typename T, typename A>
475  struct IsMatrix<BCRSMatrix<T,A> >
476  {
477  enum {
481  value = true
482  };
483  };
484 
485  template<typename T>
487  {
488  bool operator()(const T* l, const T* r)
489  {
490  return *l < *r;
491  }
492  };
493 
494 }
495 #endif
static size_type coldim(const Matrix &A, size_type c)
Definition: matrixutils.hh:252
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:280
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:41
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:405
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:155
FieldMatrix< K, n, m > Matrix
Definition: matrixutils.hh:343
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition: matrixutils.hh:177
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1856
static size_type coldim(const ThisMatrix &A)
Definition: matrixutils.hh:388
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:321
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:448
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:346
Matrix::size_type size_type
Definition: matrixutils.hh:241
ScaledIdentityMatrix< K, n > Matrix
Definition: matrixutils.hh:424
Definition: matrixutils.hh:486
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:335
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:427
Matrix::size_type size_type
Definition: matrixutils.hh:344
static void check(const Matrix &mat)
Definition: matrixutils.hh:109
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:82
B * getptr()
get pointer
Definition: bvector.hh:1051
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition: matrixutils.hh:370
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:632
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41
A generic dynamic dense matrix.
Definition: matrix.hh:24
static size_type rowdim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:373
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:400
static size_type rowdim(const ThisMatrix &A)
Definition: matrixutils.hh:383
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
Definition: bcrsmatrix.hh:72
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:85
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1850
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:437
Col col
Definition: matrixmatrix.hh:347
Matrix::size_type size_type
Definition: matrixutils.hh:398
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition: matrixutils.hh:318
MultiTypeBlockMatrix< T1, T2, T3, T4, T5, T6, T7, T8, T9 > Matrix
Definition: matrixutils.hh:129
static size_type rowdim(const Matrix &A, size_type i)
Definition: matrixutils.hh:243
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:415
Matrix & A
Definition: matrixmatrix.hh:216
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:287
iterator class for sequential access
Definition: basearray.hh:594
bool operator()(const T *l, const T *r)
Definition: matrixutils.hh:488
VariableBlockVector< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:50
size_type * getindexptr()
get pointer
Definition: bvector.hh:1057
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:361
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Matrix::block_type block_type
Definition: matrixutils.hh:240
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:27
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:452
Row row
Definition: matrixmatrix.hh:345
Matrix::size_type size_type
Definition: matrixutils.hh:319
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
BCRSMatrix< B, TA > Matrix
Definition: matrixutils.hh:239
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:626
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
std::size_t count
Definition: matrixmatrix.hh:215
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:439
derive error class from the base class in common
Definition: istlexception.hh:16
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:442
True if T is an ISTL matrix.
Definition: matrixutils.hh:458
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:351
Matrix::size_type size_type
Definition: matrixutils.hh:425
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:326
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:331
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:410
DiagonalMatrix< K, n > Matrix
Definition: matrixutils.hh:397
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:88
Matrix & mat
Definition: matrixmatrix.hh:343
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:135
Definition: matrixutils.hh:122
static size_type coldim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:378
size_type getsize() const
get size
Definition: bvector.hh:1074
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:356
ThisMatrix::size_type size_type
Definition: matrixutils.hh:371
Definition: matrixutils.hh:25
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:432
Iterator access to matrix rows
Definition: bcrsmatrix.hh:528