3 #ifndef DUNE_ISTL_MATRIXMATRIX_HH
4 #define DUNE_ISTL_MATRIXMATRIX_HH
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/timer.hh>
36 struct NonzeroPatternTraverser
41 struct NonzeroPatternTraverser<0>
43 template<
class T,
class A1,
class A2,
class F,
int n,
int m,
int k>
49 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<A.M()<<
"!="<<B.N());
54 for(Row row= A.begin(); row != A.end(); ++row) {
56 for(Col
col = row->begin();
col != row->end(); ++
col) {
59 for(BCol bcol = B[
col.index()].begin(); bcol != B[
col.index()].end(); ++bcol) {
60 func(*
col, *bcol, row.index(), bcol.index());
69 struct NonzeroPatternTraverser<1>
71 template<
class T,
class A1,
class A2,
class F,
int n,
int m,
int k>
78 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<A.N()<<
"!="<<B.N());
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());
95 struct NonzeroPatternTraverser<2>
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,
102 if(
mat.M()!=matt.M())
103 DUNE_THROW(ISTLError,
"The sizes of the matrices do not match: "<<
mat.M()<<
"!="<<matt.M());
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;
110 for(row_iterator mrow=
mat.begin(); mrow !=
mat.end(); ++mrow) {
117 for(row_iterator_t mtcol=matt.begin(); mtcol != matt.end(); ++mtcol) {
120 col_iterator_t mtrow=mtcol->begin();
121 bool funcCalled =
false;
122 for(col_iterator mcol=mrow->begin(); mcol != mrow->end(); ++mcol) {
125 for( ; mtrow != mtcol->end(); ++mtrow)
126 if(mtrow.index()>=mcol.index())
128 if(mtrow != mtcol->end() && mtrow.index()==mcol.index()) {
129 func(*mcol, *mtrow, mtcol.index());
149 template<
class T,
class A,
int n,
int m>
150 class SparsityPatternInitializer
153 enum {do_break=
true};
161 template<
class T1,
class T2>
162 void operator()(
const T1&,
const T2&, size_type j)
175 CreateIterator rowiter;
179 template<
int transpose,
class T,
class TA,
int n,
int m>
180 class MatrixInitializer
183 enum {do_break=
true};
191 template<
class T1,
class T2>
192 void operator()(
const T1&,
const T2&,
int)
203 std::size_t nonzeros()
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)
212 SparsityPatternInitializer<T, TA, n, m> sparsity(A.
createbegin());
213 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,sparsity);
221 template<
class T,
class TA,
int n,
int m>
222 class MatrixInitializer<1,T,TA,n,m>
225 enum {do_break=
false};
231 : A(A_), entries(rows)
234 template<
class T1,
class T2>
235 void operator()(
const T1&,
const T2&, size_type i, size_type j)
237 entries[i].insert(j);
246 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
247 for(Iter iter = entries.begin(); iter != entries.end(); ++iter)
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>&)
255 typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
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)
266 std::vector<std::set<size_t> > entries;
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>
275 : MatrixInitializer<1,T,TA,n,m>(A_,rows)
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)
284 typedef typename FieldMatrix<T,n,k>::size_type size_type;
286 for(size_type row=0; row<n; ++row)
288 for(size_type i=0; i < m; ++i)
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)
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) {
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)
309 typedef typename FieldMatrix<T,n,k>::size_type size_type;
310 for(size_type row=0; row<n; ++row)
312 for(size_type i=0; i < k; ++i)
318 template<
class T,
class A,
int n,
int m>
319 class EntryAccumulatorFather
322 enum {do_break=
false};
327 EntryAccumulatorFather(
Matrix& mat_)
328 :
mat(mat_), row(
mat.begin())
352 template<
class T,
class A,
int n,
int m,
int transpose>
353 class EntryAccumulator
354 :
public EntryAccumulatorFather<T,A,n,m>
360 EntryAccumulator(
Matrix& mat_)
361 : EntryAccumulatorFather<T,A,n,m>(mat_)
364 template<
class T1,
class T2>
365 void operator()(
const T1& t1,
const T2& t2, size_type i)
367 assert(this->
col.index()==i);
368 addMatMultMat(*(this->
col),t1,t2);
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>
380 EntryAccumulator(
Matrix& mat_)
381 : EntryAccumulatorFather<T,A,n,m>(mat_)
384 template<
class T1,
class T2>
385 void operator()(
const T1& t1,
const T2& t2, size_type i, size_type j)
387 addMatMultMat(this->mat[i][j], t1, t2);
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>
399 EntryAccumulator(
Matrix& mat_)
400 : EntryAccumulatorFather<T,A,n,m>(mat_)
403 template<
class T1,
class T2>
404 void operator()(
const T1& t1,
const T2& t2, size_type i, size_type j)
406 addTransposeMatMultMat(this->mat[i][j], t1, t2);
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>
418 EntryAccumulator(
Matrix& mat_)
419 : EntryAccumulatorFather<T,A,n,m>(mat_)
422 template<
class T1,
class T2>
423 void operator()(
const T1& t1,
const T2& t2, [[maybe_unused]] size_type i)
425 assert(this->
col.index()==i);
426 addMatMultTransposeMat(*this->
col,t1,t2);
431 template<
int transpose>
436 struct SizeSelector<0>
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)
442 return std::make_tuple(m1.N(), m2.M());
447 struct SizeSelector<1>
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)
453 return std::make_tuple(m1.M(), m2.M());
459 struct SizeSelector<2>
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)
465 return std::make_tuple(m1.N(), m2.N());
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)
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);
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);
486 patternInit.initPattern(mat1, mat2);
492 EntryAccumulator<T,A,n1,m1, transpose> entriesAccu(res);
493 NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,entriesAccu);
506 template<
typename M1,
typename M2>
510 template<
typename T,
int n,
int k,
int m>
516 template<
typename T,
typename A,
typename A1,
int n,
int k,
int m>
531 template<
typename M1,
typename M2>
535 template<
typename T,
int n,
int k,
int m>
541 template<
typename T,
typename A,
typename A1,
int n,
int k,
int m>
557 template<
class T,
class A,
class A1,
class A2,
int n,
int m,
int k>
561 matMultMat<2>(res,
mat, matt);
572 template<
class T,
class A,
class A1,
class A2,
int n,
int m,
int k>
576 matMultMat<0>(res,
mat, matt);
587 template<
class T,
class A,
class A1,
class A2,
int n,
int m,
int k>
591 matMultMat<1>(res,
mat, matt);
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
FieldMatrix< T, n, m > type
Definition: matrixmatrix.hh:538
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