dune-istl 3.0-git
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
17namespace Dune
18{
19
20#ifndef DOYXGEN
21 template<typename B, typename A>
22 class BCRSMatrix;
23
24 template<typename K, int n, int m>
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 {
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>
121 class MultiTypeBlockMatrix;
122
123 template<std::size_t blocklevel, std::size_t l, typename T1, typename... Args>
125 blocklevel,l>
126 {
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>
231 {};
232
233
234 template<typename B, typename TA>
236 {
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())
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>
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>
367 {
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>
394 {
396 typedef typename ThisMatrix::size_type size_type;
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>
421 {
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>
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
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Col col
Definition matrixmatrix.hh:347
Matrix & mat
Definition matrixmatrix.hh:343
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition matrixutils.hh:153
Definition basearray.hh:19
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition matrixutils.hh:175
Iterator implementation class
Definition basearray.hh:84
Definition matrixutils.hh:231
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:412
derive error class from the base class in common
Definition istlexception.hh:16
ConstIterator class for sequential access.
Definition matrix.hh:398
A generic dynamic dense matrix.
Definition matrix.hh:555
A::size_type size_type
Type for indices and sizes.
Definition matrix.hh:571
RowIterator end()
Get iterator to one beyond last row.
Definition matrix.hh:615
RowIterator begin()
Get iterator to first row.
Definition matrix.hh:609
T block_type
Export the type representing the components.
Definition matrix.hh:562
Definition matrixutils.hh:25
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition matrixutils.hh:83
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition matrixutils.hh:88
static void check(const Matrix &mat)
Definition matrixutils.hh:109
A Matrix class to support different block types.
Definition multitypeblockmatrix.hh:44
static void check(const Matrix &)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition matrixutils.hh:133
MultiTypeBlockMatrix< T1, Args... > Matrix
Definition matrixutils.hh:127
BCRSMatrix< B, TA > Matrix
Definition matrixutils.hh:237
static size_type coldim(const Matrix &A)
Definition matrixutils.hh:285
Matrix::block_type block_type
Definition matrixutils.hh:238
static size_type coldim(const Matrix &A, size_type c)
Definition matrixutils.hh:250
Matrix::size_type size_type
Definition matrixutils.hh:239
static size_type rowdim(const Matrix &A, size_type i)
Definition matrixutils.hh:241
static size_type rowdim(const Matrix &A)
Definition matrixutils.hh:278
static size_type coldim(const Matrix &A)
Definition matrixutils.hh:333
static size_type rowdim(const Matrix &, size_type)
Definition matrixutils.hh:319
static size_type rowdim(const Matrix &A)
Definition matrixutils.hh:329
Matrix::size_type size_type
Definition matrixutils.hh:317
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition matrixutils.hh:316
static size_type coldim(const Matrix &, size_type)
Definition matrixutils.hh:324
static size_type rowdim(const Matrix &, size_type)
Definition matrixutils.hh:344
static size_type coldim(const Matrix &)
Definition matrixutils.hh:359
Matrix::size_type size_type
Definition matrixutils.hh:342
FieldMatrix< K, n, m > Matrix
Definition matrixutils.hh:341
static size_type coldim(const Matrix &, size_type)
Definition matrixutils.hh:349
static size_type rowdim(const Matrix &)
Definition matrixutils.hh:354
static size_type coldim(const MatrixType &A)
Definition matrixutils.hh:386
static size_type rowdim(const MatrixType &A)
Definition matrixutils.hh:381
static size_type rowdim(const MatrixType &, size_type)
Definition matrixutils.hh:371
MatrixType::size_type size_type
Definition matrixutils.hh:369
static size_type coldim(const MatrixType &, size_type)
Definition matrixutils.hh:376
Dune::DynamicMatrix< T > MatrixType
Definition matrixutils.hh:368
static size_type coldim(const ThisMatrix &A)
Definition matrixutils.hh:413
static size_type rowdim(const ThisMatrix &A)
Definition matrixutils.hh:408
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition matrixutils.hh:395
static size_type coldim(const ThisMatrix &, size_type)
Definition matrixutils.hh:403
static size_type rowdim(const ThisMatrix &, size_type)
Definition matrixutils.hh:398
ThisMatrix::size_type size_type
Definition matrixutils.hh:396
static size_type coldim(const Matrix &, size_type)
Definition matrixutils.hh:430
Matrix::size_type size_type
Definition matrixutils.hh:423
static size_type coldim(const Matrix &)
Definition matrixutils.hh:440
static size_type rowdim(const Matrix &)
Definition matrixutils.hh:435
DiagonalMatrix< K, n > Matrix
Definition matrixutils.hh:422
static size_type rowdim(const Matrix &, size_type)
Definition matrixutils.hh:425
static size_type coldim(const Matrix &)
Definition matrixutils.hh:467
static size_type rowdim(const Matrix &, size_type)
Definition matrixutils.hh:452
static size_type coldim(const Matrix &, size_type)
Definition matrixutils.hh:457
Matrix::size_type size_type
Definition matrixutils.hh:450
ScaledIdentityMatrix< K, n > Matrix
Definition matrixutils.hh:449
static size_type rowdim(const Matrix &)
Definition matrixutils.hh:462
Test whether a type is an ISTL Matrix.
Definition matrixutils.hh:478
@ value
True if T is an ISTL matrix.
Definition matrixutils.hh:483
Definition matrixutils.hh:512
bool operator()(const T *l, const T *r)
Definition matrixutils.hh:513
A multiple of the identity matrix of static size.
Definition scaledidmatrix.hh:28