DUNE PDELab (git)

matrixmatrix.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_MATRIXMATRIX_HH
6#define DUNE_ISTL_MATRIXMATRIX_HH
7
8#include <tuple>
9
12#include <dune/common/timer.hh>
13namespace Dune
14{
15
26 namespace
27 {
28
37 template<int b>
38 struct NonzeroPatternTraverser
39 {};
40
41
42 template<>
43 struct NonzeroPatternTraverser<0>
44 {
45 template<class T,class A1, class A2, class F, int n, int m, int k>
46 static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>& A,
48 F& func)
49 {
50 if(A.M()!=B.N())
51 DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<A.M()<<"!="<<B.N());
52
53 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>::ConstRowIterator Row;
54 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>::ConstColIterator Col;
55 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;
56 for(Row row= A.begin(); row != A.end(); ++row) {
57 // Loop over all column entries
58 for(Col col = row->begin(); col != row->end(); ++col) {
59 // entry at i,k
60 // search for all nonzeros in row k
61 for(BCol bcol = B[col.index()].begin(); bcol != B[col.index()].end(); ++bcol) {
62 func(*col, *bcol, row.index(), bcol.index());
63 }
64 }
65 }
66 }
67
68 };
69
70 template<>
71 struct NonzeroPatternTraverser<1>
72 {
73 template<class T, class A1, class A2, class F, int n, int m, int k>
74 static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>& A,
76 F& func)
77 {
78
79 if(A.N()!=B.N())
80 DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<A.N()<<"!="<<B.N());
81
82 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstRowIterator Row;
83 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstColIterator Col;
84 typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;
85
86 for(Row row=A.begin(); row!=A.end(); ++row) {
87 for(Col col=row->begin(); col!=row->end(); ++col) {
88 for(BCol bcol = B[row.index()].begin(); bcol != B[row.index()].end(); ++bcol) {
89 func(*col, *bcol, col.index(), bcol.index());
90 }
91 }
92 }
93 }
94 };
95
96 template<>
97 struct NonzeroPatternTraverser<2>
98 {
99 template<class T, class A1, class A2, class F, int n, int m, int k>
100 static void traverse(const BCRSMatrix<FieldMatrix<T,n,m>,A1>& mat,
101 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt,
102 F& func)
103 {
104 if(mat.M()!=matt.M())
105 DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<mat.M()<<"!="<<matt.M());
106
107 typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstRowIterator row_iterator;
108 typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstColIterator col_iterator;
109 typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstRowIterator row_iterator_t;
110 typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstColIterator col_iterator_t;
111
112 for(row_iterator mrow=mat.begin(); mrow != mat.end(); ++mrow) {
113 //iterate over the column entries
114 // mt is a transposed matrix crs therefore it is treated as a ccs matrix
115 // and the row_iterator iterates over the columns of the transposed matrix.
116 // search the row of the transposed matrix for an entry with the same index
117 // as the mcol iterator
118
119 for(row_iterator_t mtcol=matt.begin(); mtcol != matt.end(); ++mtcol) {
120 //Search for col entries in mat that have a corresponding row index in matt
121 // (i.e. corresponding col index in the as this is the transposed matrix
122 col_iterator_t mtrow=mtcol->begin();
123 bool funcCalled = false;
124 for(col_iterator mcol=mrow->begin(); mcol != mrow->end(); ++mcol) {
125 // search
126 // TODO: This should probably be substituted by a binary search
127 for( ; mtrow != mtcol->end(); ++mtrow)
128 if(mtrow.index()>=mcol.index())
129 break;
130 if(mtrow != mtcol->end() && mtrow.index()==mcol.index()) {
131 func(*mcol, *mtrow, mtcol.index());
132 funcCalled = true;
133 // In some cases we only search for one pair, then we break here
134 // and continue with the next column.
135 if(F::do_break)
136 break;
137 }
138 }
139 // move on with func only if func was called, otherwise they might
140 // get out of sync
141 if (funcCalled)
142 func.nextCol();
143 }
144 func.nextRow();
145 }
146 }
147 };
148
149
150
151 template<class T, class A, int n, int m>
152 class SparsityPatternInitializer
153 {
154 public:
155 enum {do_break=true};
156 typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::CreateIterator CreateIterator;
157 typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::size_type size_type;
158
159 SparsityPatternInitializer(CreateIterator iter)
160 : rowiter(iter)
161 {}
162
163 template<class T1, class T2>
164 void operator()(const T1&, const T2&, size_type j)
165 {
166 rowiter.insert(j);
167 }
168
169 void nextRow()
170 {
171 ++rowiter;
172 }
173 void nextCol()
174 {}
175
176 private:
177 CreateIterator rowiter;
178 };
179
180
181 template<int transpose, class T, class TA, int n, int m>
182 class MatrixInitializer
183 {
184 public:
185 enum {do_break=true};
186 typedef typename Dune::BCRSMatrix<FieldMatrix<T,n,m>,TA> Matrix;
187 typedef typename Matrix::CreateIterator CreateIterator;
188 typedef typename Matrix::size_type size_type;
189
190 MatrixInitializer(Matrix& A_, size_type)
191 : count(0), A(A_)
192 {}
193 template<class T1, class T2>
194 void operator()(const T1&, const T2&, int)
195 {
196 ++count;
197 }
198
199 void nextCol()
200 {}
201
202 void nextRow()
203 {}
204
205 std::size_t nonzeros()
206 {
207 return count;
208 }
209
210 template<class A1, class A2, int n2, int m2, int n3, int m3>
211 void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
212 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
213 {
214 SparsityPatternInitializer<T, TA, n, m> sparsity(A.createbegin());
215 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,sparsity);
216 }
217
218 private:
219 std::size_t count;
220 Matrix& A;
221 };
222
223 template<class T, class TA, int n, int m>
224 class MatrixInitializer<1,T,TA,n,m>
225 {
226 public:
227 enum {do_break=false};
229 typedef typename Matrix::CreateIterator CreateIterator;
230 typedef typename Matrix::size_type size_type;
231
232 MatrixInitializer(Matrix& A_, size_type rows)
233 : A(A_), entries(rows)
234 {}
235
236 template<class T1, class T2>
237 void operator()(const T1&, const T2&, size_type i, size_type j)
238 {
239 entries[i].insert(j);
240 }
241
242 void nextCol()
243 {}
244
245 size_type nonzeros()
246 {
247 size_type nnz=0;
248 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
249 for(Iter iter = entries.begin(); iter != entries.end(); ++iter)
250 nnz+=(*iter).size();
251 return nnz;
252 }
253 template<class A1, class A2, int n2, int m2, int n3, int m3>
254 void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>&,
255 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>&)
256 {
257 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
258 CreateIterator citer = A.createbegin();
259 for(Iter iter = entries.begin(); iter != entries.end(); ++iter, ++citer) {
260 typedef std::set<size_t>::const_iterator SetIter;
261 for(SetIter index=iter->begin(); index != iter->end(); ++index)
262 citer.insert(*index);
263 }
264 }
265
266 private:
267 Matrix& A;
268 std::vector<std::set<size_t> > entries;
269 };
270
271 template<class T, class TA, int n, int m>
272 struct MatrixInitializer<0,T,TA,n,m>
273 : public MatrixInitializer<1,T,TA,n,m>
274 {
275 MatrixInitializer(Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>& A_,
276 typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>::size_type rows)
277 : MatrixInitializer<1,T,TA,n,m>(A_,rows)
278 {}
279 };
280
281
282 template<class T, class T1, class T2, int n, int m, int k>
283 void addMatMultTransposeMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,n,m>& mat,
284 const FieldMatrix<T2,k,m>& matt)
285 {
286 typedef typename FieldMatrix<T,n,k>::size_type size_type;
287
288 for(size_type row=0; row<n; ++row)
289 for(size_type col=0; col<k; ++col) {
290 for(size_type i=0; i < m; ++i)
291 res[row][col]+=mat[row][i]*matt[col][i];
292 }
293 }
294
295 template<class T, class T1, class T2, int n, int m, int k>
296 void addTransposeMatMultMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,m,n>& mat,
297 const FieldMatrix<T2,m,k>& matt)
298 {
299 typedef typename FieldMatrix<T,n,k>::size_type size_type;
300 for(size_type i=0; i<m; ++i)
301 for(size_type row=0; row<n; ++row) {
302 for(size_type col=0; col < k; ++col)
303 res[row][col]+=mat[i][row]*matt[i][col];
304 }
305 }
306
307 template<class T, class T1, class T2, int n, int m, int k>
308 void addMatMultMat(FieldMatrix<T,n,m>& res, const FieldMatrix<T1,n,k>& mat,
309 const FieldMatrix<T2,k,m>& matt)
310 {
311 typedef typename FieldMatrix<T,n,k>::size_type size_type;
312 for(size_type row=0; row<n; ++row)
313 for(size_type col=0; col<m; ++col) {
314 for(size_type i=0; i < k; ++i)
315 res[row][col]+=mat[row][i]*matt[i][col];
316 }
317 }
318
319
320 template<class T, class A, int n, int m>
321 class EntryAccumulatorFather
322 {
323 public:
324 enum {do_break=false};
325 typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
326 typedef typename Matrix::RowIterator Row;
327 typedef typename Matrix::ColIterator Col;
328
329 EntryAccumulatorFather(Matrix& mat_)
330 : mat(mat_), row(mat.begin())
331 {
332 mat=0;
333 col=row->begin();
334 }
335 void nextRow()
336 {
337 ++row;
338 if(row!=mat.end())
339 col=row->begin();
340 }
341
342 void nextCol()
343 {
344 ++this->col;
345 }
346 protected:
347 Matrix& mat;
348 private:
349 Row row;
350 protected:
351 Col col;
352 };
353
354 template<class T, class A, int n, int m, int transpose>
355 class EntryAccumulator
356 : public EntryAccumulatorFather<T,A,n,m>
357 {
358 public:
359 typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
360 typedef typename Matrix::size_type size_type;
361
362 EntryAccumulator(Matrix& mat_)
363 : EntryAccumulatorFather<T,A,n,m>(mat_)
364 {}
365
366 template<class T1, class T2>
367 void operator()(const T1& t1, const T2& t2, size_type i)
368 {
369 assert(this->col.index()==i);
370 addMatMultMat(*(this->col),t1,t2);
371 }
372 };
373
374 template<class T, class A, int n, int m>
375 class EntryAccumulator<T,A,n,m,0>
376 : public EntryAccumulatorFather<T,A,n,m>
377 {
378 public:
379 typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
380 typedef typename Matrix::size_type size_type;
381
382 EntryAccumulator(Matrix& mat_)
383 : EntryAccumulatorFather<T,A,n,m>(mat_)
384 {}
385
386 template<class T1, class T2>
387 void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
388 {
389 addMatMultMat(this->mat[i][j], t1, t2);
390 }
391 };
392
393 template<class T, class A, int n, int m>
394 class EntryAccumulator<T,A,n,m,1>
395 : public EntryAccumulatorFather<T,A,n,m>
396 {
397 public:
398 typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
399 typedef typename Matrix::size_type size_type;
400
401 EntryAccumulator(Matrix& mat_)
402 : EntryAccumulatorFather<T,A,n,m>(mat_)
403 {}
404
405 template<class T1, class T2>
406 void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
407 {
408 addTransposeMatMultMat(this->mat[i][j], t1, t2);
409 }
410 };
411
412 template<class T, class A, int n, int m>
413 class EntryAccumulator<T,A,n,m,2>
414 : public EntryAccumulatorFather<T,A,n,m>
415 {
416 public:
417 typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
418 typedef typename Matrix::size_type size_type;
419
420 EntryAccumulator(Matrix& mat_)
421 : EntryAccumulatorFather<T,A,n,m>(mat_)
422 {}
423
424 template<class T1, class T2>
425 void operator()(const T1& t1, const T2& t2, [[maybe_unused]] size_type i)
426 {
427 assert(this->col.index()==i);
428 addMatMultTransposeMat(*this->col,t1,t2);
429 }
430 };
431
432
433 template<int transpose>
434 struct SizeSelector
435 {};
436
437 template<>
438 struct SizeSelector<0>
439 {
440 template<class M1, class M2>
441 static std::tuple<typename M1::size_type, typename M2::size_type>
442 size(const M1& m1, const M2& m2)
443 {
444 return std::make_tuple(m1.N(), m2.M());
445 }
446 };
447
448 template<>
449 struct SizeSelector<1>
450 {
451 template<class M1, class M2>
452 static std::tuple<typename M1::size_type, typename M2::size_type>
453 size(const M1& m1, const M2& m2)
454 {
455 return std::make_tuple(m1.M(), m2.M());
456 }
457 };
458
459
460 template<>
461 struct SizeSelector<2>
462 {
463 template<class M1, class M2>
464 static std::tuple<typename M1::size_type, typename M2::size_type>
465 size(const M1& m1, const M2& m2)
466 {
467 return std::make_tuple(m1.N(), m2.N());
468 }
469 };
470
471 template<int transpose, class T, class A, class A1, class A2, int n1, int m1, int n2, int m2, int n3, int m3>
472 void matMultMat(BCRSMatrix<FieldMatrix<T,n1,m1>,A>& res, const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
473 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
474 {
475 // First step is to count the number of nonzeros
476 typename BCRSMatrix<FieldMatrix<T,n1,m1>,A>::size_type rows, cols;
477 std::tie(rows,cols)=SizeSelector<transpose>::size(mat1, mat2);
478 MatrixInitializer<transpose,T,A,n1,m1> patternInit(res, rows);
479 Timer timer;
480 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,patternInit);
481 res.setSize(rows, cols, patternInit.nonzeros());
482 res.setBuildMode(BCRSMatrix<FieldMatrix<T,n1,m1>,A>::row_wise);
483
484 //std::cout<<"Counting nonzeros took "<<timer.elapsed()<<std::endl;
485 timer.reset();
486
487 // Second step is to allocate the storage for the result and initialize the nonzero pattern
488 patternInit.initPattern(mat1, mat2);
489
490 //std::cout<<"Setting up sparsity pattern took "<<timer.elapsed()<<std::endl;
491 timer.reset();
492 // As a last step calculate the entries
493 res = 0.0;
494 EntryAccumulator<T,A,n1,m1, transpose> entriesAccu(res);
495 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,entriesAccu);
496 //std::cout<<"Calculating entries took "<<timer.elapsed()<<std::endl;
497 }
498
499 }
500
508 template<typename M1, typename M2>
510 {};
511
512 template<typename T, int n, int k, int m>
513 struct MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >
514 {
515 typedef FieldMatrix<T,n,m> type;
516 };
517
518 template<typename T, typename A, typename A1, int n, int k, int m>
519 struct MatMultMatResult<BCRSMatrix<FieldMatrix<T,n,k>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
520 {
521 typedef BCRSMatrix<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type,
522 std::allocator<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type> > type;
523 };
524
525
533 template<typename M1, typename M2>
535 {};
536
537 template<typename T, int n, int k, int m>
539 {
540 typedef FieldMatrix<T,n,m> type;
541 };
542
543 template<typename T, typename A, typename A1, int n, int k, int m>
544 struct TransposedMatMultMatResult<BCRSMatrix<FieldMatrix<T,k,n>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
545 {
546 typedef BCRSMatrix<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type,
547 std::allocator<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type> > type;
548 };
549
550
559 template<class T, class A, class A1, class A2, int n, int m, int k>
561 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, [[maybe_unused]] bool tryHard=false)
562 {
563 matMultMat<2>(res,mat, matt);
564 }
565
574 template<class T, class A, class A1, class A2, int n, int m, int k>
576 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, bool tryHard=false)
577 {
578 matMultMat<0>(res,mat, matt);
579 }
580
589 template<class T, class A, class A1, class A2, int n, int m, int k>
591 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, [[maybe_unused]] bool tryHard=false)
592 {
593 matMultMat<1>(res,mat, matt);
594 }
595
596}
597#endif
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1061
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1091
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:697
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:677
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:700
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:497
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1094
A dense n x m matrix.
Definition: fmatrix.hh:117
Implementation of the BCRSMatrix class.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void transposeMatMultMat(BCRSMatrix< FieldMatrix< T, n, m >, A > &res, const BCRSMatrix< FieldMatrix< T, k, n >, A1 > &mat, const BCRSMatrix< FieldMatrix< T, k, m >, A2 > &matt, bool tryHard=false)
Calculate product of a transposed sparse matrix with another sparse matrices ( ).
Definition: matrixmatrix.hh:590
void matMultMat(BCRSMatrix< FieldMatrix< T, n, m >, A > &res, const BCRSMatrix< FieldMatrix< T, n, k >, A1 > &mat, const BCRSMatrix< FieldMatrix< T, k, m >, A2 > &matt, bool tryHard=false)
Calculate product of two sparse matrices ( ).
Definition: matrixmatrix.hh:575
void matMultTransposeMat(BCRSMatrix< FieldMatrix< T, n, k >, A > &res, const BCRSMatrix< FieldMatrix< T, n, m >, A1 > &mat, const BCRSMatrix< FieldMatrix< T, k, m >, A2 > &matt, bool tryHard=false)
Calculate product of a sparse matrix with a transposed sparse matrices ( ).
Definition: matrixmatrix.hh:560
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Helper TMP to get the result type of a sparse matrix matrix multiplication ( )
Definition: matrixmatrix.hh:510
Helper TMP to get the result type of a sparse matrix matrix multiplication ( )
Definition: matrixmatrix.hh:535
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)