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