dune-istl  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 
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/timer.hh>
11 namespace 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};
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};
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};
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};
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:
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:
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:
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:
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:
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  {
514  };
515 
516  template<typename T, typename A, typename A1, int n, int k, int m>
518  {
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  {
539  };
540 
541  template<typename T, typename A, typename A1, int n, int k, int m>
543  {
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.
FieldMatrix< T, n, m > type
Definition: matrixmatrix.hh:513
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
Matrix::RowIterator Row
Definition: matrixmatrix.hh:324
Dune::BCRSMatrix< Dune::FieldMatrix< T, n, m >, TA > Matrix
Definition: matrixmatrix.hh:226
Matrix::size_type size_type
Definition: matrixmatrix.hh:358
Matrix::size_type size_type
Definition: matrixmatrix.hh:378
Col col
Definition: matrixmatrix.hh:349
Matrix::ColIterator Col
Definition: matrixmatrix.hh:325
Matrix::size_type size_type
Definition: matrixmatrix.hh:186
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
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, [[maybe_unused]] bool tryHard=false)
Calculate product of a sparse matrix with a transposed sparse matrices ( ).
Definition: matrixmatrix.hh:558
Matrix & mat
Definition: matrixmatrix.hh:345
Matrix::size_type size_type
Definition: matrixmatrix.hh:228
Dune::BCRSMatrix< FieldMatrix< T, n, m >, TA > Matrix
Definition: matrixmatrix.hh:184
BCRSMatrix< typename MatMultMatResult< FieldMatrix< T, n, k >, FieldMatrix< T, k, m > >::type, std::allocator< typename MatMultMatResult< FieldMatrix< T, n, k >, FieldMatrix< T, k, m > >::type > > type
Definition: matrixmatrix.hh:520
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, [[maybe_unused]] bool tryHard=false)
Calculate product of a transposed sparse matrix with another sparse matrices ( ).
Definition: matrixmatrix.hh:588
Matrix::CreateIterator CreateIterator
Definition: matrixmatrix.hh:227
Matrix::size_type size_type
Definition: matrixmatrix.hh:397
Matrix::CreateIterator CreateIterator
Definition: matrixmatrix.hh:185
BCRSMatrix< typename MatMultMatResult< FieldMatrix< T, n, k >, FieldMatrix< T, k, m > >::type, std::allocator< typename MatMultMatResult< FieldMatrix< T, n, k >, FieldMatrix< T, k, m > >::type > > type
Definition: matrixmatrix.hh:545
Matrix::size_type size_type
Definition: matrixmatrix.hh:416
BCRSMatrix< FieldMatrix< T, n, m >, A >::size_type size_type
Definition: matrixmatrix.hh:155
Definition: allocator.hh:9
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
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
Iterator access to matrix rows
Definition: bcrsmatrix.hh:577
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:955
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1062
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
Definition: matrixutils.hh:25