diff --git a/dune/istl/matrixmatrix.hh b/dune/istl/matrixmatrix.hh
index d7d9416bcae1b4f0f4aa4abf0f9504e7735b357b..c255cdcc695bf623f3e5c056f8f2c46b70ece53a 100644
--- a/dune/istl/matrixmatrix.hh
+++ b/dune/istl/matrixmatrix.hh
@@ -117,11 +117,11 @@ namespace Dune
           // 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, func.nextCol()) {
+          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
@@ -130,12 +130,17 @@ namespace Dune
                   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();
         }
diff --git a/dune/istl/test/mmtest.cc b/dune/istl/test/mmtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fae7c693c61eff5c27240ab14fd6665941b64945
--- /dev/null
+++ b/dune/istl/test/mmtest.cc
@@ -0,0 +1,41 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+#include <iostream>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/matrixmatrix.hh>
+
+int main(int argc, char** argv)
+{
+  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1> > MatrixType;
+  MatrixType m1(2,2,MatrixType::random) ,
+  m2(2,2,MatrixType::random) ,
+  res(2,2,MatrixType::random);
+
+  // initialize first matrix [1,0;0,1]
+  m1.setrowsize(0,1);
+  m1.setrowsize(1,1);
+  m1.endrowsizes();
+  m1.addindex(0,0);
+  m1.addindex(1,1);
+  m1.endindices();
+  m1[0][0] = 1;
+  m1[1][1] = 1;
+  // initialize second matrix [0,1;1,0]
+  m2.setrowsize(0,1);
+  m2.setrowsize(1,1);
+  m2.endrowsizes();
+  m2.addindex(0,1);
+  m2.addindex(1,0);
+  m2.endindices();
+  m2[0][1] = 1;
+  m2[1][0] = 1;
+  Dune::printmatrix(std::cout, m1, "m1", "");
+  Dune::printmatrix(std::cout, m2, "m2", "");
+  Dune::matMultTransposeMat(res, m1, m2);
+  Dune::printmatrix(std::cout, res, "res", "");
+  return 0;
+}