Dune Core Modules (2.7.1)

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, size_type i)
424 {
426 assert(this->col.index()==i);
427 addMatMultTransposeMat(*this->col,t1,t2);
428 }
429 };
430
431
432 template<int transpose>
433 struct SizeSelector
434 {};
435
436 template<>
437 struct SizeSelector<0>
438 {
439 template<class M1, class M2>
440 static std::tuple<typename M1::size_type, typename M2::size_type>
441 size(const M1& m1, const M2& m2)
442 {
443 return std::make_tuple(m1.N(), m2.M());
444 }
445 };
446
447 template<>
448 struct SizeSelector<1>
449 {
450 template<class M1, class M2>
451 static std::tuple<typename M1::size_type, typename M2::size_type>
452 size(const M1& m1, const M2& m2)
453 {
454 return std::make_tuple(m1.M(), m2.M());
455 }
456 };
457
458
459 template<>
460 struct SizeSelector<2>
461 {
462 template<class M1, class M2>
463 static std::tuple<typename M1::size_type, typename M2::size_type>
464 size(const M1& m1, const M2& m2)
465 {
466 return std::make_tuple(m1.N(), m2.N());
467 }
468 };
469
470 template<int transpose, class T, class A, class A1, class A2, int n1, int m1, int n2, int m2, int n3, int m3>
471 void matMultMat(BCRSMatrix<FieldMatrix<T,n1,m1>,A>& res, const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
472 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
473 {
474 // First step is to count the number of nonzeros
475 typename BCRSMatrix<FieldMatrix<T,n1,m1>,A>::size_type rows, cols;
476 std::tie(rows,cols)=SizeSelector<transpose>::size(mat1, mat2);
477 MatrixInitializer<transpose,T,A,n1,m1> patternInit(res, rows);
478 Timer timer;
479 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,patternInit);
480 res.setSize(rows, cols, patternInit.nonzeros());
481 res.setBuildMode(BCRSMatrix<FieldMatrix<T,n1,m1>,A>::row_wise);
482
483 //std::cout<<"Counting nonzeros took "<<timer.elapsed()<<std::endl;
484 timer.reset();
485
486 // Second step is to allocate the storage for the result and initialize the nonzero pattern
487 patternInit.initPattern(mat1, mat2);
488
489 //std::cout<<"Setting up sparsity pattern took "<<timer.elapsed()<<std::endl;
490 timer.reset();
491 // As a last step calculate the entries
492 res = 0.0;
493 EntryAccumulator<T,A,n1,m1, transpose> entriesAccu(res);
494 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,entriesAccu);
495 //std::cout<<"Calculating entries took "<<timer.elapsed()<<std::endl;
496 }
497
498 }
499
507 template<typename M1, typename M2>
509 {};
510
511 template<typename T, int n, int k, int m>
512 struct MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >
513 {
514 typedef FieldMatrix<T,n,m> type;
515 };
516
517 template<typename T, typename A, typename A1, int n, int k, int m>
518 struct MatMultMatResult<BCRSMatrix<FieldMatrix<T,n,k>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
519 {
520 typedef BCRSMatrix<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type,
521 std::allocator<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type> > type;
522 };
523
524
532 template<typename M1, typename M2>
534 {};
535
536 template<typename T, int n, int k, int m>
538 {
539 typedef FieldMatrix<T,n,m> type;
540 };
541
542 template<typename T, typename A, typename A1, int n, int k, int m>
543 struct TransposedMatMultMatResult<BCRSMatrix<FieldMatrix<T,k,n>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
544 {
545 typedef BCRSMatrix<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type,
546 std::allocator<typename MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >::type> > type;
547 };
548
549
558 template<class T, class A, class A1, class A2, int n, int m, int k>
560 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, bool tryHard=false)
561 {
562 DUNE_UNUSED_PARAMETER(tryHard);
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, bool tryHard=false)
592 {
593 DUNE_UNUSED_PARAMETER(tryHard);
594 matMultMat<1>(res,mat, matt);
595 }
596
597}
598#endif
Implementation of the BCRSMatrix class.
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1023
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:426
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1053
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:660
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:640
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:663
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:460
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
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_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#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: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:559
Dune namespace.
Definition: alignedallocator.hh:14
Helper TMP to get the result type of a sparse matrix matrix multiplication ( )
Definition: matrixmatrix.hh:509
Helper TMP to get the result type of a sparse matrix matrix multiplication ( )
Definition: matrixmatrix.hh:534
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)