dune-istl  2.2.1
matrixutils.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 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>
13 #include"istlexception.hh"
14 
15 namespace Dune
16 {
17 
18 #ifndef DOYXGEN
19  template<typename B, typename A>
20  class BCRSMatrix;
21 
22  template<typename K, int n, int m>
23  class FieldMatrix;
24 
25  template<class T, class A>
26  class Matrix;
27 #endif
28 
38  namespace
39  {
40 
41  template<int i>
42  struct NonZeroCounter
43  {
44  template<class M>
45  static typename M::size_type count(const M& matrix)
46  {
47  typedef typename M::ConstRowIterator RowIterator;
48 
49  RowIterator endRow = matrix.end();
50  typename M::size_type nonZeros = 0;
51 
52  for(RowIterator row = matrix.begin(); row != endRow; ++row){
53  typedef typename M::ConstColIterator Entry;
54  Entry endEntry = row->end();
55  for(Entry entry = row->begin(); entry != endEntry; ++entry){
56  nonZeros += NonZeroCounter<i-1>::count(*entry);
57  }
58  }
59  return nonZeros;
60  }
61  };
62 
63  template<>
64  struct NonZeroCounter<1>
65  {
66  template<class M>
67  static typename M::size_type count(const M& matrix)
68  {
69  return matrix.N()*matrix.M();
70  }
71  };
72 
73  }
74 
79  template<class Matrix, std::size_t blocklevel, std::size_t l=blocklevel>
81  {
86  static void check(const Matrix& mat)
87  {
88 #ifdef DUNE_ISTL_WITH_CHECKING
89  typedef typename Matrix::ConstRowIterator Row;
90  typedef typename Matrix::ConstColIterator Entry;
91  for(Row row = mat.begin(); row!=mat.end(); ++row){
92  Entry diagonal = row->find(row.index());
93  if(diagonal==row->end())
94  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
95  <<" at block recursion level "<<l-blocklevel);
96  else
98  }
99 #endif
100  }
101  };
102 
103  template<class Matrix, std::size_t l>
105  {
106  static void check(const Matrix& mat)
107  {
108  typedef typename Matrix::ConstRowIterator Row;
109  for(Row row = mat.begin(); row!=mat.end(); ++row){
110  if(row->find(row.index())==row->end())
111  DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
112  <<" at block recursion level "<<l);
113  }
114  }
115  };
116 
117  template<typename T1, typename T2, typename T3, typename T4, typename T5,
118  typename T6, typename T7, typename T8, typename T9>
119  class MultiTypeBlockMatrix;
120 
121  template<typename T1, typename T2, typename T3, typename T4, typename T5,
122  typename T6, typename T7, typename T8, typename T9, std::size_t blocklevel, std::size_t l>
123  struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>,
124  blocklevel,l>
125  {
126  typedef MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> Matrix;
127 
132  static void check(const Matrix& mat)
133  {
134 #ifdef DUNE_ISTL_WITH_CHECKING
135  // TODO Implement check
136 #endif
137  }
138  };
139 
151  template<class M>
152  inline int countNonZeros(const M& matrix)
153  {
155  }
156  /*
157  template<class M>
158  struct ProcessOnFieldsOfMatrix
159  */
160 
162  namespace
163  {
164  struct CompPair{
165  template<class G,class M>
166  bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2)
167  {
168  return p1.first<p2.first;
169  }
170  };
171 
172  }
173  template<class M, class C>
174  void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os)
175  {
176  typedef typename C::ParallelIndexSet::const_iterator IIter;
177  typedef typename C::OwnerSet OwnerSet;
178  typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
179 
180  GlobalIndex gmax=0;
181 
182  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
183  idx!=eidx; ++idx)
184  gmax=std::max(gmax,idx->global());
185 
186  gmax=ooc.communicator().max(gmax);
187  ooc.buildGlobalLookup();
188 
189  for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
190  idx!=eidx; ++idx){
191  if(OwnerSet::contains(idx->local().attribute()))
192  {
193  typedef typename M::block_type Block;
194 
195  std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
196 
197  // sort rows
198  typedef typename M::ConstColIterator CIter;
199  for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
200  c!=cend; ++c){
201  const typename C::ParallelIndexSet::IndexPair* pair
202  =ooc.globalLookup().pair(c.index());
203  assert(pair);
204  entries.insert(std::make_pair(pair->global(), *c));
205  }
206 
207  //wait until its the rows turn.
208  GlobalIndex rowidx = idx->global();
209  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
210  while(cur!=rowidx)
211  cur=ooc.communicator().min(rowidx);
212 
213  // print rows
214  typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
215  for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
216  os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl;
217 
218 
219  }
220  }
221 
222  ooc.freeGlobalLookup();
223  // Wait until everybody is finished
224  GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
225  while(cur!=ooc.communicator().min(cur));
226  }
227 
228  template<typename M>
230  {
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>
314  struct MatrixDimension<BCRSMatrix<FieldMatrix<B,n,m> ,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  {
341  typedef FieldMatrix<K,n,m> Matrix;
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<typename K, int n, int m, typename TA>
366  struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
367  {
370 
371  static size_type rowdim(const ThisMatrix& A, size_type r)
372  {
373  return n;
374  }
375 
376  static size_type coldim(const ThisMatrix& A, size_type r)
377  {
378  return m;
379  }
380 
381  static size_type rowdim(const ThisMatrix& A)
382  {
383  return A.N()*n;
384  }
385 
386  static size_type coldim(const ThisMatrix& A)
387  {
388  return A.M()*m;
389  }
390  };
391 
395  template<typename T>
396  struct IsMatrix
397  {
398  enum{
402  value = false
403  };
404  };
405 
406  template<typename T>
407  struct IsMatrix<DenseMatrix<T> >
408  {
409  enum{
413  value = true
414  };
415  };
416 
417 
418  template<typename T, typename A>
419  struct IsMatrix<BCRSMatrix<T,A> >
420  {
421  enum{
425  value = true
426  };
427  };
428 
429 }
430 #endif