This file is indexed.

/usr/include/dune/istl/matrixmatrix.hh is in libdune-istl-dev 2.4.1-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_MATRIXMATRIX_HH
#define DUNE_ISTL_MATRIXMATRIX_HH
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/tuples.hh>
#include <dune/common/timer.hh>
namespace Dune
{

  /**
   * @addtogroup ISTL_SPMV
   *
   * @{
   */
  /** @file
   * @author Markus Blatt
   * @brief provides functions for sparse matrix matrix multiplication.
   */

  namespace
  {

    /**
     * @brief Traverses over the nonzero pattern of the matrix-matrix product.
     *
     * Template parameter b is used to select the matrix product:
     * <dt>0</dt><dd>\f$A\cdot B\f$</dd>
     * <dt>1</dt><dd>\f$A^T\cdot B\f$</dd>
     * <dt>2</dt><dd>\f$A\cdot B^T\f$</dd>
     */
    template<int b>
    struct NonzeroPatternTraverser
    {};


    template<>
    struct NonzeroPatternTraverser<0>
    {
      template<class T,class A1, class A2, class F, int n, int m, int k>
      static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>& A,
                           const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>& B,
                           F& func)
      {
        if(A.M()!=B.N())
          DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<A.M()<<"!="<<B.N());

        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>::ConstRowIterator Row;
        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,k>,A1>::ConstColIterator Col;
        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;
        for(Row row= A.begin(); row != A.end(); ++row) {
          // Loop over all column entries
          for(Col col = row->begin(); col != row->end(); ++col) {
            // entry at i,k
            // search for all nonzeros in row k
            for(BCol bcol = B[col.index()].begin(); bcol != B[col.index()].end(); ++bcol) {
              func(*col, *bcol, row.index(), bcol.index());
            }
          }
        }
      }

    };

    template<>
    struct NonzeroPatternTraverser<1>
    {
      template<class T, class A1, class A2, class F, int n, int m, int k>
      static void traverse(const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>& A,
                           const Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>& B,
                           F& func)
      {

        if(A.N()!=B.N())
          DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<A.N()<<"!="<<B.N());

        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstRowIterator Row;
        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,n>,A1>::ConstColIterator Col;
        typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<T,k,m>,A2>::ConstColIterator BCol;

        for(Row row=A.begin(); row!=A.end(); ++row) {
          for(Col col=row->begin(); col!=row->end(); ++col) {
            for(BCol bcol  = B[row.index()].begin(); bcol !=  B[row.index()].end(); ++bcol) {
              func(*col, *bcol, col.index(), bcol.index());
            }
          }
        }
      }
    };

    template<>
    struct NonzeroPatternTraverser<2>
    {
      template<class T, class A1, class A2, class F, int n, int m, int k>
      static void traverse(const BCRSMatrix<FieldMatrix<T,n,m>,A1>& mat,
                           const BCRSMatrix<FieldMatrix<T,k,m>,A2>& matt,
                           F& func)
      {
        if(mat.M()!=matt.M())
          DUNE_THROW(ISTLError, "The sizes of the matrices do not match: "<<mat.M()<<"!="<<matt.M());

        typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstRowIterator row_iterator;
        typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A1>::ConstColIterator col_iterator;
        typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstRowIterator row_iterator_t;
        typedef typename BCRSMatrix<FieldMatrix<T,k,m>,A2>::ConstColIterator col_iterator_t;

        for(row_iterator mrow=mat.begin(); mrow != mat.end(); ++mrow) {
          //iterate over the column entries
          // mt is a transposed matrix crs therefore it is treated as a ccs matrix
          // and the row_iterator iterates over the columns of the transposed matrix.
          // search the row of the transposed matrix for an entry with the same index
          // as the mcol iterator

          for(row_iterator_t mtcol=matt.begin(); mtcol != matt.end(); ++mtcol) {
            //Search for col entries in mat that have a corrsponding row index in matt
            // (i.e. corresponding col index in the as this is the transposed matrix
            col_iterator_t mtrow=mtcol->begin();
            bool funcCalled = false;
            for(col_iterator mcol=mrow->begin(); mcol != mrow->end(); ++mcol) {
              // search
              // TODO: This should probably be substituted by a binary search
              for( ; mtrow != mtcol->end(); ++mtrow)
                if(mtrow.index()>=mcol.index())
                  break;
              if(mtrow != mtcol->end() && mtrow.index()==mcol.index()) {
                func(*mcol, *mtrow, mtcol.index());
                funcCalled = true;
                // In some cases we only search for one pair, then we break here
                // and continue with the next column.
                if(F::do_break)
                  break;
              }
            }
            // move on with func only if func was called, otherwise they might
            // get out of sync
            if (funcCalled)
              func.nextCol();
          }
          func.nextRow();
        }
      }
    };



    template<class T, class A, int n, int m>
    class SparsityPatternInitializer
    {
    public:
      enum {do_break=true};
      typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::CreateIterator CreateIterator;
      typedef typename BCRSMatrix<FieldMatrix<T,n,m>,A>::size_type size_type;

      SparsityPatternInitializer(CreateIterator iter)
        : rowiter(iter)
      {}

      template<class T1, class T2>
      void operator()(const T1&, const T2&, size_type j)
      {
        rowiter.insert(j);
      }

      void nextRow()
      {
        ++rowiter;
      }
      void nextCol()
      {}

    private:
      CreateIterator rowiter;
    };


    template<int transpose, class T, class TA, int n, int m>
    class MatrixInitializer
    {
    public:
      enum {do_break=true};
      typedef typename Dune::BCRSMatrix<FieldMatrix<T,n,m>,TA> Matrix;
      typedef typename Matrix::CreateIterator CreateIterator;
      typedef typename Matrix::size_type size_type;

      MatrixInitializer(Matrix& A_, size_type)
        : count(0), A(A_)
      {}
      template<class T1, class T2>
      void operator()(const T1&, const T2&, int)
      {
        ++count;
      }

      void nextCol()
      {}

      void nextRow()
      {}

      std::size_t nonzeros()
      {
        return count;
      }

      template<class A1, class A2, int n2, int m2, int n3, int m3>
      void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
                       const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
      {
        SparsityPatternInitializer<T, TA, n, m> sparsity(A.createbegin());
        NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,sparsity);
      }

    private:
      std::size_t count;
      Matrix& A;
    };

    template<class T, class TA, int n, int m>
    class MatrixInitializer<1,T,TA,n,m>
    {
    public:
      enum {do_break=false};
      typedef Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA> Matrix;
      typedef typename Matrix::CreateIterator CreateIterator;
      typedef typename Matrix::size_type size_type;

      MatrixInitializer(Matrix& A_, size_type rows)
        :  A(A_), entries(rows)
      {}

      template<class T1, class T2>
      void operator()(const T1&, const T2&, size_type i, size_type j)
      {
        entries[i].insert(j);
      }

      void nextCol()
      {}

      size_type nonzeros()
      {
        size_type nnz=0;
        typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
        for(Iter iter = entries.begin(); iter != entries.end(); ++iter)
          nnz+=(*iter).size();
        return nnz;
      }
      template<class A1, class A2, int n2, int m2, int n3, int m3>
      void initPattern(const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>&,
                       const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>&)
      {
        typedef typename std::vector<std::set<size_t> >::const_iterator Iter;
        CreateIterator citer = A.createbegin();
        for(Iter iter = entries.begin(); iter != entries.end(); ++iter, ++citer) {
          typedef std::set<size_t>::const_iterator SetIter;
          for(SetIter index=iter->begin(); index != iter->end(); ++index)
            citer.insert(*index);
        }
      }

    private:
      Matrix& A;
      std::vector<std::set<size_t> > entries;
    };

    template<class T, class TA, int n, int m>
    struct MatrixInitializer<0,T,TA,n,m>
      : public MatrixInitializer<1,T,TA,n,m>
    {
      MatrixInitializer(Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>& A_,
                        typename Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,TA>::size_type rows)
        : MatrixInitializer<1,T,TA,n,m>(A_,rows)
      {}
    };


    template<class T, class T1, class T2, int n, int m, int k>
    void addMatMultTransposeMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,n,m>& mat,
                                const FieldMatrix<T2,k,m>& matt)
    {
      typedef typename FieldMatrix<T,n,k>::size_type size_type;

      for(size_type row=0; row<n; ++row)
        for(size_type col=0; col<k; ++col) {
          for(size_type i=0; i < m; ++i)
            res[row][col]+=mat[row][i]*matt[col][i];
        }
    }

    template<class T, class T1, class T2, int n, int m, int k>
    void addTransposeMatMultMat(FieldMatrix<T,n,k>& res, const FieldMatrix<T1,m,n>& mat,
                                const FieldMatrix<T2,m,k>& matt)
    {
      typedef typename FieldMatrix<T,n,k>::size_type size_type;
      for(size_type i=0; i<m; ++i)
        for(size_type row=0; row<n; ++row) {
          for(size_type col=0; col < k; ++col)
            res[row][col]+=mat[i][row]*matt[i][col];
        }
    }

    template<class T, class T1, class T2, int n, int m, int k>
    void addMatMultMat(FieldMatrix<T,n,m>& res, const FieldMatrix<T1,n,k>& mat,
                       const FieldMatrix<T2,k,m>& matt)
    {
      typedef typename FieldMatrix<T,n,k>::size_type size_type;
      for(size_type row=0; row<n; ++row)
        for(size_type col=0; col<m; ++col) {
          for(size_type i=0; i < k; ++i)
            res[row][col]+=mat[row][i]*matt[i][col];
        }
    }


    template<class T, class A, int n, int m>
    class EntryAccumulatorFather
    {
    public:
      enum {do_break=false};
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
      typedef typename Matrix::RowIterator Row;
      typedef typename Matrix::ColIterator Col;

      EntryAccumulatorFather(Matrix& mat_)
        : mat(mat_), row(mat.begin())
      {
        mat=0;
        col=row->begin();
      }
      void nextRow()
      {
        ++row;
        if(row!=mat.end())
          col=row->begin();
      }

      void nextCol()
      {
        ++this->col;
      }
    protected:
      Matrix& mat;
    private:
      Row row;
    protected:
      Col col;
    };

    template<class T, class A, int n, int m, int transpose>
    class EntryAccumulator
      : public EntryAccumulatorFather<T,A,n,m>
    {
    public:
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
      typedef typename Matrix::size_type size_type;

      EntryAccumulator(Matrix& mat_)
        : EntryAccumulatorFather<T,A,n,m>(mat_)
      {}

      template<class T1, class T2>
      void operator()(const T1& t1, const T2& t2, size_type i)
      {
        assert(this->col.index()==i);
        addMatMultMat(*(this->col),t1,t2);
      }
    };

    template<class T, class A, int n, int m>
    class EntryAccumulator<T,A,n,m,0>
      : public EntryAccumulatorFather<T,A,n,m>
    {
    public:
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
      typedef typename Matrix::size_type size_type;

      EntryAccumulator(Matrix& mat_)
        : EntryAccumulatorFather<T,A,n,m>(mat_)
      {}

      template<class T1, class T2>
      void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
      {
        addMatMultMat(this->mat[i][j], t1, t2);
      }
    };

    template<class T, class A, int n, int m>
    class EntryAccumulator<T,A,n,m,1>
      : public EntryAccumulatorFather<T,A,n,m>
    {
    public:
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
      typedef typename Matrix::size_type size_type;

      EntryAccumulator(Matrix& mat_)
        : EntryAccumulatorFather<T,A,n,m>(mat_)
      {}

      template<class T1, class T2>
      void operator()(const T1& t1, const T2& t2, size_type i, size_type j)
      {
        addTransposeMatMultMat(this->mat[i][j], t1, t2);
      }
    };

    template<class T, class A, int n, int m>
    class EntryAccumulator<T,A,n,m,2>
      : public EntryAccumulatorFather<T,A,n,m>
    {
    public:
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
      typedef typename Matrix::size_type size_type;

      EntryAccumulator(Matrix& mat_)
        : EntryAccumulatorFather<T,A,n,m>(mat_)
      {}

      template<class T1, class T2>
      void operator()(const T1& t1, const T2& t2, size_type i)
      {
        DUNE_UNUSED_PARAMETER(i);
        assert(this->col.index()==i);
        addMatMultTransposeMat(*this->col,t1,t2);
      }
    };


    template<int transpose>
    struct SizeSelector
    {};

    template<>
    struct SizeSelector<0>
    {
      template<class M1, class M2>
      static tuple<typename M1::size_type, typename M2::size_type>
      size(const M1& m1, const M2& m2)
      {
        return make_tuple(m1.N(), m2.M());
      }
    };

    template<>
    struct SizeSelector<1>
    {
      template<class M1, class M2>
      static tuple<typename M1::size_type, typename M2::size_type>
      size(const M1& m1, const M2& m2)
      {
        return make_tuple(m1.M(), m2.M());
      }
    };


    template<>
    struct SizeSelector<2>
    {
      template<class M1, class M2>
      static tuple<typename M1::size_type, typename M2::size_type>
      size(const M1& m1, const M2& m2)
      {
        return make_tuple(m1.N(), m2.N());
      }
    };

    template<int transpose, class T, class A, class A1, class A2, int n1, int m1, int n2, int m2, int n3, int m3>
    void matMultMat(BCRSMatrix<FieldMatrix<T,n1,m1>,A>& res, const BCRSMatrix<FieldMatrix<T,n2,m2>,A1>& mat1,
                    const BCRSMatrix<FieldMatrix<T,n3,m3>,A2>& mat2)
    {
      // First step is to count the number of nonzeros
      typename BCRSMatrix<FieldMatrix<T,n1,m1>,A>::size_type rows, cols;
      tie(rows,cols)=SizeSelector<transpose>::size(mat1, mat2);
      MatrixInitializer<transpose,T,A,n1,m1> patternInit(res, rows);
      Timer timer;
      NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,patternInit);
      res.setSize(rows, cols, patternInit.nonzeros());
      res.setBuildMode(BCRSMatrix<FieldMatrix<T,n1,m1>,A>::row_wise);

      //std::cout<<"Counting nonzeros took "<<timer.elapsed()<<std::endl;
      timer.reset();

      // Second step is to allocate the storage for the result and initialize the nonzero pattern
      patternInit.initPattern(mat1, mat2);

      //std::cout<<"Setting up sparsity pattern took "<<timer.elapsed()<<std::endl;
      timer.reset();
      // As a last step calculate the entries
      res = 0.0;
      EntryAccumulator<T,A,n1,m1, transpose> entriesAccu(res);
      NonzeroPatternTraverser<transpose>::traverse(mat1,mat2,entriesAccu);
      //std::cout<<"Calculating entries took "<<timer.elapsed()<<std::endl;
    }

  }

  /**
   * @brief Helper TMP to get the result type of a sparse matrix matrix multiplication (\f$C=A*B\f$)
   *
   * The type of matrix C will be stored as the associated type MatMultMatResult::type.
   * @tparam M1 The type of matrix A.
   * @tparam M2 The type of matrix B.
   */
  template<typename M1, typename M2>
  struct MatMultMatResult
  {};

  template<typename T, int n, int k, int m>
  struct MatMultMatResult<FieldMatrix<T,n,k>,FieldMatrix<T,k,m> >
  {
    typedef FieldMatrix<T,n,m> type;
  };

  template<typename T, typename A, typename A1, int n, int k, int m>
  struct MatMultMatResult<BCRSMatrix<FieldMatrix<T,n,k>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
  {
    typedef 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;
  };


  /**
   * @brief Helper TMP to get the result type of a sparse matrix matrix multiplication (\f$C=A*B\f$)
   *
   * The type of matrix C will be stored as the associated type MatMultMatResult::type.
   * @tparam M1 The type of matrix A.
   * @tparam M2 The type of matrix B.
   */
  template<typename M1, typename M2>
  struct TransposedMatMultMatResult
  {};

  template<typename T, int n, int k, int m>
  struct TransposedMatMultMatResult<FieldMatrix<T,k,n>,FieldMatrix<T,k,m> >
  {
    typedef FieldMatrix<T,n,m> type;
  };

  template<typename T, typename A, typename A1, int n, int k, int m>
  struct TransposedMatMultMatResult<BCRSMatrix<FieldMatrix<T,k,n>,A >,BCRSMatrix<FieldMatrix<T,k,m>,A1 > >
  {
    typedef 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;
  };


  /**
   * @brief Calculate product of a sparse matrix with a transposed sparse matrices (\f$C=A*B^T\f$).
   *
   * @param res Matrix for the result of the computation.
   * @param mat Matrix A.
   * @param matt Matrix B, which will be transposed before the multiplication.
   * @param tryHard <i>ignored</i>
   */
  template<class T, class A, class A1, class A2, int n, int m, int k>
  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)
  {
    DUNE_UNUSED_PARAMETER(tryHard);
    matMultMat<2>(res,mat, matt);
  }

  /**
   * @brief Calculate product of two sparse matrices (\f$C=A*B\f$).
   *
   * @param res Matrix for the result of the computation.
   * @param mat Matrix A.
   * @param matt Matrix B.
   * @param tryHard <i>ignored</i>
   */
  template<class T, class A, class A1, class A2, int n, int m, int k>
  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)
  {
    matMultMat<0>(res,mat, matt);
  }

  /**
   * @brief Calculate product of a transposed sparse matrix with another sparse matrices (\f$C=A^T*B\f$).
   *
   * @param res Matrix for the result of the computation.
   * @param mat Matrix A, which will be transposed before the multiplication.
   * @param matt Matrix B.
   * @param tryHard <i>ignored</i>
   */
  template<class T, class A, class A1, class A2, int n, int m, int k>
  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)
  {
    DUNE_UNUSED_PARAMETER(tryHard);
    matMultMat<1>(res,mat, matt);
  }

}
#endif