dune-istl  2.2.1
matrixmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_MATRIXMATRIX_HH
2 #define DUNE_MATRIXMATRIX_HH
4 #include <dune/common/fmatrix.hh>
5 #include <dune/common/tuples.hh>
6 #include <dune/common/timer.hh>
7 namespace Dune
8 {
9 
20  namespace
21  {
22 
31  template<int b>
32  struct NonzeroPatternTraverser
33  {};
34 
35 
36  template<>
37  struct NonzeroPatternTraverser<0>
38  {
39  template<class T,class A1, class A2, class F, int n, int m, int k>
40  static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>& A,
41  const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>& B,
42  F& func)
43  {
44  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::size_type size_type;
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>::ConstRowIterator BRow;
52  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;
53  for(Row row= A.begin(); row != A.end(); ++row){
54  // Loop over all column entries
55  for(Col col = row->begin(); col != row->end(); ++col){
56  // entry at i,k
57  // search for all nonzeros in row k
58  for(BCol bcol = B[col.index()].begin(); bcol != B[col.index()].end(); ++bcol){
59  func(*col, *bcol, row.index(), bcol.index());
60  }
61  }
62  }
63  }
64 
65  };
66 
67  template<>
68  struct NonzeroPatternTraverser<1>
69  {
70  template<class T, class A1, class A2, class F, int n, int m, int k>
71  static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>& A,
72  const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>& B,
73  F& func)
74  {
75 
76  if(A.N()!=B.N())
77  DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<A.N()<<"!="<<B.N());
78 
79  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstRowIterator Row;
80  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstColIterator Col;
81  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;
82  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::size_type size_t1;
83  typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::size_type size_t2;
84 
85  for(Row row=A.begin(); row!=A.end(); ++row){
86  for(Col col=row->begin(); col!=row->end(); ++col){
87  for(BCol bcol = B[row.index()].begin(); bcol != B[row.index()].end(); ++bcol){
88  func(*col, *bcol, col.index(), bcol.index());
89  }
90  }
91  }
92  }
93  };
94 
95  template<>
96  struct NonzeroPatternTraverser<2>
97  {
98  template<class T, class A1, class A2, class F, int n, int m, int k>
99  static void traverse(const BCRSMatrix<FieldMatrix<T,n,m>,A1>& mat,
100  const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt,
101  F& func)
102  {
103  if(mat.M()!=matt.M())
104  DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<mat.N()<<"!="<<matt.N());
105 
106  typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstRowIterator row_iterator;
107  typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstColIterator col_iterator;
108  typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstRowIterator row_iterator_t;
109  typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstColIterator col_iterator_t;
110 
111  for(row_iterator mrow=mat.begin(); mrow != mat.end(); ++mrow){
112  //iterate over the column entries
113  // mt is a transposed matrix crs therefore it is treated as a ccs matrix
114  // and the row_iterator iterates over the columns of the transposed matrix.
115  // search the row of the transposed matrix for an entry with the same index
116  // as the mcol iterator
117 
118  for(row_iterator_t mtcol=matt.begin(); mtcol != matt.end(); ++mtcol){
119  //Search for col entries in mat that have a corrsponding row index in matt
120  // (i.e. corresponding col index in the as this is the transposed matrix
121  col_iterator_t mtrow=mtcol->begin();
122  bool funcCalled = false;
123  for(col_iterator mcol=mrow->begin(); mcol != mrow->end(); ++mcol){
124  // search
125  // TODO: This should probably be substituted by a binary search
126  for( ;mtrow != mtcol->end(); ++mtrow)
127  if(mtrow.index()>=mcol.index())
128  break;
129  if(mtrow != mtcol->end() && mtrow.index()==mcol.index()){
130  func(*mcol, *mtrow, mtcol.index());
131  funcCalled = true;
132  // In some cases we only search for one pair, then we break here
133  // and continue with the next column.
134  if(F::do_break)
135  break;
136  }
137  }
138  // move on with func only if func was called, otherwise they might
139  // get out of sync
140  if (funcCalled)
141  func.nextCol();
142  }
143  func.nextRow();
144  }
145  }
146  };
147 
148 
149 
150  template<class T, class A, int n, int m>
151  class SparsityPatternInitializer
152  {
153  public:
154  enum{do_break=true};
155  typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::CreateIterator CreateIterator;
156  typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::size_type size_type;
157 
158  SparsityPatternInitializer(CreateIterator iter)
159  : rowiter(iter)
160  {}
161 
162  template<class T1, class T2>
163  void operator()(const T1& t1, const T2& t2, size_type j)
164  {
165  rowiter.insert(j);
166  }
167 
168  void nextRow()
169  {
170  ++rowiter;
171  }
172  void nextCol()
173  {}
174 
175  private:
176  CreateIterator rowiter;
177  };
178 
179 
180  template<int transpose, class T, class TA, int n, int m>
181  class MatrixInitializer
182  {
183  public:
184  enum{do_break=true};
187  typedef typename Matrix::size_type size_type;
188 
189  MatrixInitializer(Matrix& A_, size_type rows)
190  : count(0), A(A_)
191  {}
192  template<class T1, class T2>
193  void operator()(const T1& t1, const T2& t2, int j)
194  {
195  ++count;
196  }
197 
198  void nextCol()
199  {}
200 
201  void nextRow()
202  {}
203 
204  std::size_t nonzeros()
205  {
206  return count;
207  }
208 
209  template<class A1, class A2, int n2, int m2, int n3, int m3>
210  void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
211  const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
212  {
213  SparsityPatternInitializer<T, TA, n, m> sparsity(A.createbegin());
214  NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,sparsity);
215  }
216 
217  private:
218  std::size_t count;
219  Matrix& A;
220  };
221 
222  template<class T, class TA, int n, int m>
223  class MatrixInitializer<1,T,TA,n,m>
224  {
225  public:
226  enum{do_break=false};
229  typedef typename Matrix::size_type size_type;
230 
231  MatrixInitializer(Matrix& A_, size_type rows)
232  : A(A_), entries(rows)
233  {}
234 
235  template<class T1, class T2>
236  void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
237  {
238  entries[i].insert(j);
239  }
240 
241  void nextCol()
242  {}
243 
244  size_type nonzeros()
245  {
246  size_type nnz=0;
247  typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
248  for(Iter iter = entries.begin(); iter != entries.end(); ++iter)
249  nnz+=(*iter).size();
250  return nnz;
251  }
252  template<class A1, class A2, int n2, int m2, int n3, int m3>
253  void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
254  const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
255  {
256  typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
257  CreateIterator citer = A.createbegin();
258  for(Iter iter = entries.begin(); iter != entries.end(); ++iter, ++citer){
259  typedef std::set<size_t>::const_iterator SetIter;
260  for(SetIter index=iter->begin(); index != iter->end(); ++index)
261  citer.insert(*index);
262  }
263  }
264 
265  private:
266  Matrix& A;
267  std::vector<std::set<size_t> > entries;
268  };
269 
270  template<class T, class TA, int n, int m>
271  struct MatrixInitializer<0,T,TA,n,m>
272  : public MatrixInitializer<1,T,TA,n,m>
273  {
274  MatrixInitializer(Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>& A_,
275  typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>::size_type rows)
276  :MatrixInitializer<1,T,TA,n,m>(A_,rows)
277  {}
278  };
279 
280 
281  template<class T, class T1, class T2, int n, int m, int k>
282  void addMatMultTransposeMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,n,m>& mat,
283  const FieldMatrix<T2,k,m>& matt)
284  {
285  typedef typename FieldMatrix<T,n,k>::size_type size_type;
286 
287  for(size_type row=0; row<n; ++row)
288  for(size_type col=0; col<k;++col){
289  for(size_type i=0; i < m; ++i)
290  res[row][col]+=mat[row][i]*matt[col][i];
291  }
292  }
293 
294  template<class T, class T1, class T2, int n, int m, int k>
295  void addTransposeMatMultMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,m,n>& mat,
296  const FieldMatrix<T2,m,k>& matt)
297  {
298  typedef typename FieldMatrix<T,n,k>::size_type size_type;
299  for(size_type i=0; i<m; ++i)
300  for(size_type row=0; row<n;++row){
301  for(size_type col=0; col < k; ++col)
302  res[row][col]+=mat[i][row]*matt[i][col];
303  }
304  }
305 
306  template<class T, class T1, class T2, int n, int m, int k>
307  void addMatMultMat(FieldMatrix<T,n,m>& res, const FieldMatrix<T1,n,k>& mat,
308  const FieldMatrix<T2,k,m>& matt)
309  {
310  typedef typename FieldMatrix<T,n,k>::size_type size_type;
311  for(size_type row=0; row<n; ++row)
312  for(size_type col=0; col<m;++col){
313  for(size_type i=0; i < k; ++i)
314  res[row][col]+=mat[row][i]*matt[i][col];
315  }
316  }
317 
318 
319  template<class T, class A, int n, int m>
320  class EntryAccumulatorFather
321  {
322  public:
323  enum{do_break=false};
324  typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
325  typedef typename Matrix::RowIterator Row;
326  typedef typename Matrix::ColIterator Col;
327 
328  EntryAccumulatorFather(Matrix& mat_)
329  :mat(mat_), row(mat.begin())
330  {
331  mat=0;
332  col=row->begin();
333  }
334  void nextRow()
335  {
336  ++row;
337  if(row!=mat.end())
338  col=row->begin();
339  }
340 
341  void nextCol()
342  {
343  ++this->col;
344  }
345  protected:
346  Matrix& mat;
347  private:
348  Row row;
349  protected:
350  Col col;
351  };
352 
353  template<class T, class A, int n, int m, int transpose>
354  class EntryAccumulator
355  : public EntryAccumulatorFather<T,A,n,m>
356  {
357  public:
358  typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
359  typedef typename Matrix::size_type size_type;
360 
361  EntryAccumulator(Matrix& mat_)
362  : EntryAccumulatorFather<T,A,n,m>(mat_)
363  {}
364 
365  template<class T1, class T2>
366  void operator()(const T1& t1, const T2& t2, size_type i)
367  {
368  assert(this->col.index()==i);
369  addMatMultMat(*(this->col),t1,t2);
370  }
371  };
372 
373  template<class T, class A, int n, int m>
374  class EntryAccumulator<T,A,n,m,0>
375  : public EntryAccumulatorFather<T,A,n,m>
376  {
377  public:
378  typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
379  typedef typename Matrix::size_type size_type;
380 
381  EntryAccumulator(Matrix& mat_)
382  : EntryAccumulatorFather<T,A,n,m>(mat_)
383  {}
384 
385  template<class T1, class T2>
386  void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
387  {
388  addMatMultMat(this->mat[i][j], t1, t2);
389  }
390  };
391 
392  template<class T, class A, int n, int m>
393  class EntryAccumulator<T,A,n,m,1>
394  : public EntryAccumulatorFather<T,A,n,m>
395  {
396  public:
397  typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
398  typedef typename Matrix::size_type size_type;
399 
400  EntryAccumulator(Matrix& mat_)
401  : EntryAccumulatorFather<T,A,n,m>(mat_)
402  {}
403 
404  template<class T1, class T2>
405  void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
406  {
407  addTransposeMatMultMat(this->mat[i][j], t1, t2);
408  }
409  };
410 
411  template<class T, class A, int n, int m>
412  class EntryAccumulator<T,A,n,m,2>
413  : public EntryAccumulatorFather<T,A,n,m>
414  {
415  public:
416  typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
417  typedef typename Matrix::size_type size_type;
418 
419  EntryAccumulator(Matrix& mat_)
420  : EntryAccumulatorFather<T,A,n,m>(mat_)
421  {}
422 
423  template<class T1, class T2>
424  void operator()(const T1& t1, const T2& t2, size_type i)
425  {
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 
437  template<>
438  struct SizeSelector<0>
439  {
440  template<class M1, class M2>
441  static tuple<typename M1::size_type, typename M2::size_type>
442  size(const M1& m1, const M2& m2)
443  {
444  return make_tuple(m1.N(), m2.M());
445  }
446  };
447 
448  template<>
449  struct SizeSelector<1>
450  {
451  template<class M1, class M2>
452  static tuple<typename M1::size_type, typename M2::size_type>
453  size(const M1& m1, const M2& m2)
454  {
455  return 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 tuple<typename M1::size_type, typename M2::size_type>
465  size(const M1& m1, const M2& m2)
466  {
467  return 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  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  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 
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  {
522  FieldMatrix<T,k,m> >::type,A> type;
523  };
524 
533  template<class T, class A, class A1, class A2, int n, int m, int k>
534  void matMultTransposeMat(BCRSMatrix<FieldMatrix<T,n,k>,A>& res, const BCRSMatrix<FieldMatrix<T,n,m>,A1>& mat,
535  const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, bool tryHard=false)
536  {
537  matMultMat<2>(res,mat, matt);
538  }
539 
548  template<class T, class A, class A1, class A2, int n, int m, int k>
549  void matMultMat(BCRSMatrix<FieldMatrix<T,n,m>,A>& res, const BCRSMatrix<FieldMatrix<T,n,k>,A1>& mat,
550  const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, bool tryHard=false)
551  {
552  matMultMat<0>(res,mat, matt);
553  }
554 
563  template<class T, class A, class A1, class A2, int n, int m, int k>
564  void transposeMatMultMat(BCRSMatrix<FieldMatrix<T,n,m>,A>& res, const BCRSMatrix<FieldMatrix<T,k,n>,A1>& mat,
565  const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt, bool tryHard=false)
566  {
567  matMultMat<1>(res,mat, matt);
568  }
569 
570 }
571 #endif
572