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>
32 struct NonzeroPatternTraverser
37 struct NonzeroPatternTraverser<0>
39 template<
class T,
class A1,
class A2,
class F,
int n,
int m,
int k>
47 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<
A.M()<<
"!="<<B.N());
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;
58 for(BCol bcol = B[
col.index()].begin(); bcol != B[
col.index()].end(); ++bcol){
59 func(*
col, *bcol,
row.index(), bcol.index());
68 struct NonzeroPatternTraverser<1>
70 template<
class T,
class A1,
class A2,
class F,
int n,
int m,
int k>
77 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<
A.N()<<
"!="<<B.N());
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;
87 for(BCol bcol = B[
row.index()].begin(); bcol != B[
row.index()].end(); ++bcol){
88 func(*
col, *bcol,
col.index(), bcol.index());
96 struct NonzeroPatternTraverser<2>
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,
103 if(
mat.M()!=matt.M())
104 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<
mat.N()<<
"!="<<matt.N());
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;
111 for(row_iterator mrow=
mat.begin(); mrow !=
mat.end(); ++mrow){
118 for(row_iterator_t mtcol=matt.begin(); mtcol != matt.end(); ++mtcol, func.nextCol()){
121 col_iterator_t mtrow=mtcol->
begin();
123 for(col_iterator mcol=mrow->begin(); mcol != mrow->end(); ++mcol){
126 for( ;mtrow != mtcol->end(); ++mtrow)
127 if(mtrow.index()>=mcol.index())
129 if(mtrow != mtcol->end() && mtrow.index()==mcol.index()){
130 func(*mcol, *mtrow, mtcol.index());
145 template<
class T,
class A,
int n,
int m>
146 class SparsityPatternInitializer
157 template<
class T1,
class T2>
158 void operator()(
const T1& t1,
const T2& t2, size_type j)
175 template<
int transpose,
class T,
class TA,
int n,
int m>
176 class MatrixInitializer
184 MatrixInitializer(
Matrix& A_, size_type rows)
187 template<
class T1,
class T2>
188 void operator()(
const T1& t1,
const T2& t2,
int j)
199 std::size_t nonzeros()
204 template<
class A1,
class A2,
int n2,
int m2,
int n3,
int m3>
205 void initPattern(
const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
206 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
208 SparsityPatternInitializer<T, TA, n, m> sparsity(
A.
createbegin());
209 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,sparsity);
217 template<
class T,
class TA,
int n,
int m>
218 class MatrixInitializer<1,T,TA,n,m>
221 enum{do_break=
false};
226 MatrixInitializer(
Matrix& A_, size_type rows)
227 :
A(A_), entries(rows)
230 template<
class T1,
class T2>
231 void operator()(
const T1& t1,
const T2& t2, size_type i, size_type j)
233 entries[i].insert(j);
242 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
243 for(Iter iter = entries.begin(); iter != entries.end(); ++iter)
247 template<
class A1,
class A2,
int n2,
int m2,
int n3,
int m3>
248 void initPattern(
const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
249 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
251 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
253 for(Iter iter = entries.begin(); iter != entries.end(); ++iter, ++citer){
254 typedef std::set<size_t>::const_iterator SetIter;
262 std::vector<std::set<size_t> > entries;
265 template<
class T,
class TA,
int n,
int m>
266 struct MatrixInitializer<0,T,TA,n,m>
267 :
public MatrixInitializer<1,T,TA,n,m>
271 :MatrixInitializer<1,T,TA,n,m>(A_,rows)
276 template<
class T,
class T1,
class T2,
int n,
int m,
int k>
277 void addMatMultTransposeMat(FieldMatrix<T,n,k>& res,
const FieldMatrix<T1,n,m>&
mat,
278 const FieldMatrix<T2,k,m>& matt)
280 typedef typename FieldMatrix<T,n,k>::size_type size_type;
284 for(size_type i=0; i < m; ++i)
289 template<
class T,
class T1,
class T2,
int n,
int m,
int k>
290 void addTransposeMatMultMat(FieldMatrix<T,n,k>& res,
const FieldMatrix<T1,m,n>& mat,
291 const FieldMatrix<T2,m,k>& matt)
293 typedef typename FieldMatrix<T,n,k>::size_type size_type;
294 for(size_type i=0; i<m; ++i)
301 template<
class T,
class T1,
class T2,
int n,
int m,
int k>
302 void addMatMultMat(FieldMatrix<T,n,m>& res,
const FieldMatrix<T1,n,k>& mat,
303 const FieldMatrix<T2,k,m>& matt)
305 typedef typename FieldMatrix<T,n,k>::size_type size_type;
308 for(size_type i=0; i < k; ++i)
314 template<
class T,
class A,
int n,
int m>
315 class EntryAccumulatorFather
318 enum{do_break=
false};
319 typedef BCRSMatrix<FieldMatrix<T,n,m>,
A>
Matrix;
323 EntryAccumulatorFather(
Matrix& mat_)
324 :mat(mat_),
row(mat.begin())
348 template<
class T,
class A,
int n,
int m,
int transpose>
349 class EntryAccumulator
350 :
public EntryAccumulatorFather<T,A,n,m>
353 typedef BCRSMatrix<FieldMatrix<T,n,m>,
A>
Matrix;
356 EntryAccumulator(
Matrix& mat_)
357 : EntryAccumulatorFather<T,
A,n,m>(mat_)
360 template<
class T1,
class T2>
361 void operator()(
const T1& t1,
const T2& t2, size_type i)
363 assert(this->
col.index()==i);
364 addMatMultMat(*(this->
col),t1,t2);
368 template<
class T,
class A,
int n,
int m>
369 class EntryAccumulator<T,
A,n,m,0>
370 :
public EntryAccumulatorFather<T,A,n,m>
373 typedef BCRSMatrix<FieldMatrix<T,n,m>,
A>
Matrix;
376 EntryAccumulator(
Matrix& mat_)
377 : EntryAccumulatorFather<T,
A,n,m>(mat_)
380 template<
class T1,
class T2>
381 void operator()(
const T1& t1,
const T2& t2, size_type i, size_type j)
383 addMatMultMat(this->mat[i][j], t1, t2);
387 template<
class T,
class A,
int n,
int m>
388 class EntryAccumulator<T,
A,n,m,1>
389 :
public EntryAccumulatorFather<T,A,n,m>
392 typedef BCRSMatrix<FieldMatrix<T,n,m>,
A>
Matrix;
395 EntryAccumulator(
Matrix& mat_)
396 : EntryAccumulatorFather<T,
A,n,m>(mat_)
399 template<
class T1,
class T2>
400 void operator()(
const T1& t1,
const T2& t2, size_type i, size_type j)
402 addTransposeMatMultMat(this->mat[i][j], t1, t2);
406 template<
class T,
class A,
int n,
int m>
407 class EntryAccumulator<T,
A,n,m,2>
408 :
public EntryAccumulatorFather<T,A,n,m>
411 typedef BCRSMatrix<FieldMatrix<T,n,m>,
A>
Matrix;
414 EntryAccumulator(
Matrix& mat_)
415 : EntryAccumulatorFather<T,
A,n,m>(mat_)
418 template<
class T1,
class T2>
419 void operator()(
const T1& t1,
const T2& t2, size_type i)
421 assert(this->
col.index()==i);
422 addMatMultTransposeMat(*this->
col,t1,t2);
427 template<
int transpose>
433 struct SizeSelector<0>
435 template<
class M1,
class M2>
436 static tuple<typename M1::size_type, typename M2::size_type>
437 size(
const M1& m1,
const M2& m2)
439 return make_tuple(m1.N(), m2.M());
444 struct SizeSelector<1>
446 template<
class M1,
class M2>
447 static tuple<typename M1::size_type, typename M2::size_type>
448 size(
const M1& m1,
const M2& m2)
450 return make_tuple(m1.M(), m2.M());
456 struct SizeSelector<2>
458 template<
class M1,
class M2>
459 static tuple<typename M1::size_type, typename M2::size_type>
460 size(
const M1& m1,
const M2& m2)
462 return make_tuple(m1.N(), m2.N());
466 template<
int transpose,
class T,
class A,
class A1,
class A2,
int n1,
int m1,
int n2,
int m2,
int n3,
int m3>
467 void matMultMat(BCRSMatrix<FieldMatrix<T,n1,m1>,
A>& res,
const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
468 const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
471 typename BCRSMatrix<FieldMatrix<T,n1,m1>,
A>::size_type rows, cols;
472 tie(rows,cols)=SizeSelector<transpose>::size(mat1, mat2);
473 MatrixInitializer<transpose,T,A,n1,m1> patternInit(res, rows);
475 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,patternInit);
476 res.setSize(rows, cols, patternInit.nonzeros());
477 res.setBuildMode(BCRSMatrix<FieldMatrix<T,n1,m1>,
A>::row_wise);
479 std::cout<<
"Counting nonzeros took "<<timer.elapsed()<<std::endl;
483 patternInit.initPattern(mat1, mat2);
485 std::cout<<
"Setting up sparsity pattern took "<<timer.elapsed()<<std::endl;
488 EntryAccumulator<T,A,n1,m1, transpose> entriesAccu(res);
489 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,entriesAccu);
490 std::cout<<
"Calculating entries took "<<timer.elapsed()<<std::endl;
502 template<
typename M1,
typename M2>
507 template<
typename T,
int n,
int k,
int m>
510 typedef FieldMatrix<T,n,m>
type;
513 template<
typename T,
typename A,
typename A1,
int n,
int k,
int m>
528 template<
class T,
class A,
class A1,
class A2,
int n,
int m,
int k>
530 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt,
bool tryHard=
false)
532 matMultMat<2>(res,
mat, matt);
543 template<
class T,
class A,
class A1,
class A2,
int n,
int m,
int k>
545 const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt,
bool tryHard=
false)
547 matMultMat<0>(res,
mat, matt);
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)
562 matMultMat<1>(res,
mat, matt);