Dune Core Modules (2.8.0)

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