Dune Core Modules (2.8.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_ISTL_MATRIXUTILS_HH
4#define DUNE_ISTL_MATRIXUTILS_HH
5
6#include <set>
7#include <vector>
8#include <limits>
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>
25 class FieldMatrix;
26
27 template<class T, class A>
28 class Matrix;
29#endif
30
44 template<class Matrix, std::size_t blocklevel, std::size_t l=blocklevel>
46 {
51 static void check([[maybe_unused]] const Matrix& mat)
52 {
53#ifdef DUNE_ISTL_WITH_CHECKING
54 typedef typename Matrix::ConstRowIterator Row;
55 typedef typename Matrix::ConstColIterator Entry;
56 for(Row row = mat.begin(); row!=mat.end(); ++row) {
57 Entry diagonal = row->find(row.index());
58 if(diagonal==row->end())
59 DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
60 <<" at block recursion level "<<l-blocklevel);
61 else{
62 auto m = Impl::asMatrix(*diagonal);
63 CheckIfDiagonalPresent<decltype(m),blocklevel-1,l>::check(m);
64 }
65 }
66#endif
67 }
68 };
69
70 template<class Matrix, std::size_t l>
71 struct CheckIfDiagonalPresent<Matrix,0,l>
72 {
73 static void check(const Matrix& mat)
74 {
75 typedef typename Matrix::ConstRowIterator Row;
76 for(Row row = mat.begin(); row!=mat.end(); ++row) {
77 if(row->find(row.index())==row->end())
78 DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index()
79 <<" at block recursion level "<<l);
80 }
81 }
82 };
83
84 template<typename FirstRow, typename... Args>
85 class MultiTypeBlockMatrix;
86
87 template<std::size_t blocklevel, std::size_t l, typename T1, typename... Args>
88 struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,Args...>,
89 blocklevel,l>
90 {
91 typedef MultiTypeBlockMatrix<T1,Args...> Matrix;
92
97 static void check(const Matrix& /* mat */)
98 {
99#ifdef DUNE_ISTL_WITH_CHECKING
100 // TODO Implement check
101#endif
102 }
103 };
104
116 template<class M>
117 inline auto countNonZeros(const M&,
118 [[maybe_unused]] typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae = nullptr)
119 {
120 return 1;
121 }
122
123 template<class M>
124 inline auto countNonZeros(const M& matrix,
125 [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
126 {
127 typename M::size_type nonZeros = 0;
128 for(auto&& row : matrix)
129 for(auto&& entry : row)
130 nonZeros += countNonZeros(entry);
131 return nonZeros;
132 }
133
134 /*
135 template<class M>
136 struct ProcessOnFieldsOfMatrix
137 */
138
140 namespace
141 {
142 struct CompPair {
143 template<class G,class M>
144 bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2) const
145 {
146 return p1.first<p2.first;
147 }
148 };
149
150 }
151 template<class M, class C>
152 void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os)
153 {
154 typedef typename C::ParallelIndexSet::const_iterator IIter;
155 typedef typename C::OwnerSet OwnerSet;
156 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
157
158 GlobalIndex gmax=0;
159
160 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
161 idx!=eidx; ++idx)
162 gmax=std::max(gmax,idx->global());
163
164 gmax=ooc.communicator().max(gmax);
165 ooc.buildGlobalLookup();
166
167 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
168 idx!=eidx; ++idx) {
169 if(OwnerSet::contains(idx->local().attribute()))
170 {
171 typedef typename M::block_type Block;
172
173 std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
174
175 // sort rows
176 typedef typename M::ConstColIterator CIter;
177 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
178 c!=cend; ++c) {
179 const typename C::ParallelIndexSet::IndexPair* pair
180 =ooc.globalLookup().pair(c.index());
181 assert(pair);
182 entries.insert(std::make_pair(pair->global(), *c));
183 }
184
185 //wait until its the rows turn.
186 GlobalIndex rowidx = idx->global();
188 while(cur!=rowidx)
189 cur=ooc.communicator().min(rowidx);
190
191 // print rows
192 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
193 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
194 os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl;
195
196
197 }
198 }
199
200 ooc.freeGlobalLookup();
201 // Wait until everybody is finished
203 while(cur!=ooc.communicator().min(cur)) ;
204 }
205
206 // Default implementation for scalar types
207 template<typename M>
208 struct MatrixDimension
209 {
210 static_assert(IsNumber<M>::value, "MatrixDimension is not implemented for this type!");
211
212 static auto rowdim(const M& A)
213 {
214 return 1;
215 }
216
217 static auto coldim(const M& A)
218 {
219 return 1;
220 }
221 };
222
223 // Default implementation for scalar types
224 template<typename B, typename TA>
225 struct MatrixDimension<Matrix<B,TA> >
226 {
227 using block_type = typename Matrix<B,TA>::block_type;
228 using size_type = typename Matrix<B,TA>::size_type;
229
230 static size_type rowdim (const Matrix<B,TA>& A, size_type i)
231 {
232 return MatrixDimension<block_type>::rowdim(A[i][0]);
233 }
234
235 static size_type coldim (const Matrix<B,TA>& A, size_type c)
236 {
237 return MatrixDimension<block_type>::coldim(A[0][c]);
238 }
239
240 static size_type rowdim (const Matrix<B,TA>& A)
241 {
242 size_type nn=0;
243 for (size_type i=0; i<A.N(); i++)
244 nn += rowdim(A,i);
245 return nn;
246 }
247
248 static size_type coldim (const Matrix<B,TA>& A)
249 {
250 size_type nn=0;
251 for (size_type i=0; i<A.M(); i++)
252 nn += coldim(A,i);
253 return nn;
254 }
255 };
256
257
258 template<typename B, typename TA>
259 struct MatrixDimension<BCRSMatrix<B,TA> >
260 {
261 typedef BCRSMatrix<B,TA> Matrix;
262 typedef typename Matrix::block_type block_type;
263 typedef typename Matrix::size_type size_type;
264
265 static size_type rowdim (const Matrix& A, size_type i)
266 {
267 const B* row = A.r[i].getptr();
268 if(row)
269 return MatrixDimension<block_type>::rowdim(*row);
270 else
271 return 0;
272 }
273
274 static size_type coldim (const Matrix& A, size_type c)
275 {
276 // find an entry in column c
277 if (A.nnz_ > 0)
278 {
279 for (size_type k=0; k<A.nnz_; k++) {
280 if (A.j_.get()[k] == c) {
281 return MatrixDimension<block_type>::coldim(A.a[k]);
282 }
283 }
284 }
285 else
286 {
287 for (size_type i=0; i<A.N(); i++)
288 {
289 size_type* j = A.r[i].getindexptr();
290 B* a = A.r[i].getptr();
291 for (size_type k=0; k<A.r[i].getsize(); k++)
292 if (j[k]==c) {
293 return MatrixDimension<block_type>::coldim(a[k]);
294 }
295 }
296 }
297
298 // not found
299 return 0;
300 }
301
302 static size_type rowdim (const Matrix& A){
303 size_type nn=0;
304 for (size_type i=0; i<A.N(); i++)
305 nn += rowdim(A,i);
306 return nn;
307 }
308
309 static size_type coldim (const Matrix& A){
310 typedef typename Matrix::ConstRowIterator ConstRowIterator;
311 typedef typename Matrix::ConstColIterator ConstColIterator;
312
313 // The following code has a complexity of nnz, and
314 // typically a very small constant.
315 //
316 std::vector<size_type> coldims(A.M(),
318
319 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
320 for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
321 // only compute blocksizes we don't already have
322 if (coldims[col.index()]==std::numeric_limits<size_type>::max())
323 coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
324
325 size_type sum = 0;
326 for (typename std::vector<size_type>::iterator it=coldims.begin();
327 it!=coldims.end(); ++it)
328 // skip rows for which no coldim could be determined
329 if ((*it)>=0)
330 sum += *it;
331
332 return sum;
333 }
334 };
335
336
337 template<typename B, int n, int m, typename TA>
338 struct MatrixDimension<BCRSMatrix<FieldMatrix<B,n,m> ,TA> >
339 {
340 typedef BCRSMatrix<FieldMatrix<B,n,m> ,TA> Matrix;
341 typedef typename Matrix::size_type size_type;
342
343 static size_type rowdim (const Matrix& /*A*/, size_type /*i*/)
344 {
345 return n;
346 }
347
348 static size_type coldim (const Matrix& /*A*/, size_type /*c*/)
349 {
350 return m;
351 }
352
353 static size_type rowdim (const Matrix& A) {
354 return A.N()*n;
355 }
356
357 static size_type coldim (const Matrix& A) {
358 return A.M()*m;
359 }
360 };
361
362 template<typename K, int n, int m>
363 struct MatrixDimension<FieldMatrix<K,n,m> >
364 {
365 typedef FieldMatrix<K,n,m> Matrix;
366 typedef typename Matrix::size_type size_type;
367
368 static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
369 {
370 return 1;
371 }
372
373 static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
374 {
375 return 1;
376 }
377
378 static size_type rowdim(const Matrix& /*A*/)
379 {
380 return n;
381 }
382
383 static size_type coldim(const Matrix& /*A*/)
384 {
385 return m;
386 }
387 };
388
389 template <class T>
390 struct MatrixDimension<Dune::DynamicMatrix<T> >
391 {
392 typedef Dune::DynamicMatrix<T> MatrixType;
393 typedef typename MatrixType::size_type size_type;
394
395 static size_type rowdim(const MatrixType& /*A*/, size_type /*r*/)
396 {
397 return 1;
398 }
399
400 static size_type coldim(const MatrixType& /*A*/, size_type /*r*/)
401 {
402 return 1;
403 }
404
405 static size_type rowdim(const MatrixType& A)
406 {
407 return A.N();
408 }
409
410 static size_type coldim(const MatrixType& A)
411 {
412 return A.M();
413 }
414 };
415
416 template<typename K, int n, int m, typename TA>
417 struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
418 {
419 typedef Matrix<FieldMatrix<K,n,m>, TA> ThisMatrix;
420 typedef typename ThisMatrix::size_type size_type;
421
422 static size_type rowdim(const ThisMatrix& /*A*/, size_type /*r*/)
423 {
424 return n;
425 }
426
427 static size_type coldim(const ThisMatrix& /*A*/, size_type /*r*/)
428 {
429 return m;
430 }
431
432 static size_type rowdim(const ThisMatrix& A)
433 {
434 return A.N()*n;
435 }
436
437 static size_type coldim(const ThisMatrix& A)
438 {
439 return A.M()*m;
440 }
441 };
442
443 template<typename K, int n>
444 struct MatrixDimension<DiagonalMatrix<K,n> >
445 {
446 typedef DiagonalMatrix<K,n> Matrix;
447 typedef typename Matrix::size_type size_type;
448
449 static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
450 {
451 return 1;
452 }
453
454 static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
455 {
456 return 1;
457 }
458
459 static size_type rowdim(const Matrix& /*A*/)
460 {
461 return n;
462 }
463
464 static size_type coldim(const Matrix& /*A*/)
465 {
466 return n;
467 }
468 };
469
470 template<typename K, int n>
471 struct MatrixDimension<ScaledIdentityMatrix<K,n> >
472 {
473 typedef ScaledIdentityMatrix<K,n> Matrix;
474 typedef typename Matrix::size_type size_type;
475
476 static size_type rowdim(const Matrix& /*A*/, size_type /*r*/)
477 {
478 return 1;
479 }
480
481 static size_type coldim(const Matrix& /*A*/, size_type /*r*/)
482 {
483 return 1;
484 }
485
486 static size_type rowdim(const Matrix& /*A*/)
487 {
488 return n;
489 }
490
491 static size_type coldim(const Matrix& /*A*/)
492 {
493 return n;
494 }
495 };
496
500 template<typename T>
501 struct IsMatrix
502 {
503 enum {
507 value = false
508 };
509 };
510
511 template<typename T>
512 struct IsMatrix<DenseMatrix<T> >
513 {
514 enum {
518 value = true
519 };
520 };
521
522
523 template<typename T, typename A>
524 struct IsMatrix<BCRSMatrix<T,A> >
525 {
526 enum {
530 value = true
531 };
532 };
533
534 template<typename T>
535 struct PointerCompare
536 {
537 bool operator()(const T* l, const T* r)
538 {
539 return *l < *r;
540 }
541 };
542
543}
544#endif
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:59
derive error class from the base class in common
Definition: istlexception.hh:17
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:584
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
This file implements a quadratic diagonal matrix of fixed size.
Traits for type conversions and type information.
This file implements a dense matrix with dynamic numbers of rows and columns.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto countNonZeros(const M &, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
Implements a scalar matrix view wrapper around an existing scalar.
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:46
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:51
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:507
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)