.gitlab-ci.yml

 ---
+include:
+  - "https://gitlab.dune-project.org/core/ci-config/raw/master/config/common/master.yml"
+  - "https://gitlab.dune-project.org/core/ci-config/raw/master/jobs/common/master.yml"
+
 before_script:
   - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
 
@@ -10,31 +14,3 @@ variables:
 # number of cores. We solve the issue by disallowing OpenMP to allocate more
 # than one thread.
   OMP_NUM_THREADS: 1
-
-debian:10  gcc:c++17:
-  image: duneci/base:10
-  script: duneci-standard-test
-  variables: {DUNECI_OPTS: /duneci/opts.gcc.c++17}
-  tags: [duneci]
-
-debian:9--gcc:
-  image: duneci/base:9
-  script: duneci-standard-test
-  tags: [duneci]
-
-debian:9--clang:
-  image: duneci/base:9
-  script: duneci-standard-test
-  variables: {DUNECI_OPTS: /duneci/opts.clang}
-  tags: [duneci]
-
-ubuntu:16.04--gcc:
-  image: duneci/base:16.04
-  script: duneci-standard-test
-  tags: [duneci]
-
-ubuntu:16.04--clang:
-  image: duneci/base:16.04
-  script: duneci-standard-test
-  variables: {DUNECI_OPTS: /duneci/opts.clang}
-  tags: [duneci]

CHANGELOG.md

 # Master (will become release 2.7)
 
+- `BDMatrix` objects now have the method `solve`, which implements that
+  canonical way to solve block-diagonal linear systems.
+
 - Deprecated the preconditioner implementations `SeqILU0` and `SeqILUn`.
   Use `SeqILU` instead, which implements incomplete LU decomposition
   of any order.
 
+- The class `VariableBlockVector::CreateIterator` is a true STL output iterator now.
+  This means that you can use STL algorithms like `std::fill` or `std::copy`
+  to set the block sizes.
+
+- Support for SuiteSparse's CHOLMOD providing a sparse Cholesky
+  factorization.
+
+- `MultiTypeBlockVector<Args...>` now inherits the constructors from its
+  parent type (`std::tuple<Args...>`). This means you can now also construct
+  `MultiTypeBlockVector`s from values or references of BlockVectors.
+
+- `MultiTypeBlockVector::count()` is now `const`
+
+- All matrix and vector classes can now be instantiated with number types
+  directly (A number type is any type for which `Dune::IsNumber<T>::value`
+  is true).  For example, you can now use `BlockVector<double>` instead of
+  the more cumbersome `BlockVector<FieldVector<double,1> >`.  Similarly, you can use
+  `BCRSMatrix<double>` instead of `BCRSMatrix<FieldMatrix<double,1,1>>`.
+  The old forms still work, and `FieldVector` and `FieldMatrix` types with
+  a single entry can still be cast to their `field_type`.  Therefore, the
+  change is completely backward-compatible.
+
+- Added a right-preconditioned flexible restarted GMRes solver
+
+- The UMFPack binding use the long int functions to compute larger systems.
+  With the \*_dl_\* versions instead of the \*_di_\* versions UMFPACK will not
+  have a memory limit of just 2 GiB.
+
+- Deprecated support for SuperLU 4.x. It will be removed after Dune 2.7.
+
 # Release 2.6
 
 - `BDMatrix` objects can now be constructed and assigned from `std::initializer_list`.

cmake/modules/DuneIstlMacros.cmake

@@ -10,7 +10,7 @@ find_package(SuperLU)
 include(AddSuperLUFlags)
 find_package(ARPACKPP)
 include(AddARPACKPPFlags)
-find_package(SuiteSparse OPTIONAL_COMPONENTS LDL SPQR UMFPACK)
+find_package(SuiteSparse OPTIONAL_COMPONENTS CHOLMOD LDL SPQR UMFPACK)
 include(AddSuiteSparseFlags)
 
 # enable / disable backwards compatibility w.r.t. category

cmake/modules/FindSuperLU.cmake

@@ -153,6 +153,10 @@ else()
     set(SUPERLU_WITH_VERSION "SuperLU <= 4.2 and >= 4.0" CACHE STRING
       "Human readable string containing SuperLU version information.")
   endif()
+
+  if(NOT SUPERLU_MIN_VERSION_5)
+    message(AUTHOR_WARNING "Support for SuperLU 4.x is deprecated an will be removed after Dune 2.7")
+  endif()
 endif()
 
 # behave like a CMake module is supposed to behave

doc/CMakeLists.txt

 add_subdirectory(doxygen)
 
-dune_add_latex_document(istl.tex
-  FORCE_DVI
-  BIBFILES istl.bib
-  IMAGES blockstructure.eps)
-create_doc_install(istl.pdf
-  ${CMAKE_INSTALL_DOCDIR} istl_safepdf)
+dune_add_latex_document(
+  SOURCE istl.tex
+  FATHER_TARGET doc
+  INSTALL ${CMAKE_INSTALL_DOCDIR})

doc/istl.tex

 \documentclass[11pt]{article}
 \usepackage{multicol}
-\usepackage{ifthen}
-%\usepackage{multitoc}
-%\usepackage{german}
-%\usepackage{bibgerm}
 \usepackage{amsmath}
 \usepackage{amsfonts}
 \usepackage{color}
 \usepackage{hyperref}
 \usepackage[dvips]{epsfig}
 \usepackage{subfigure}
-\usepackage[dvips]{graphicx}
 \usepackage[a4paper,body={148mm,240mm,nohead}]{geometry}
 \usepackage[ansinew]{inputenc}
 \usepackage{listings}
 \lstset{language=C++, basicstyle=\ttfamily,
   stringstyle=\ttfamily, commentstyle=\it, extendedchars=true}
-
-\newif\ifpdf
-\ifnum\ifx\pdfoutput\undefined0\else\pdfoutput\fi<1
-\pdffalse % we are not running PDFLaTeX
-\else
-\pdftrue % we are running PDFLaTeX
-\fi
-
-\ifpdf
-\usepackage[pdftex]{graphicx}
-\else
 \usepackage{graphicx}
-\fi
 
-\ifpdf
-\DeclareGraphicsExtensions{.pdf, .jpg, .tif}
-\else
-\DeclareGraphicsExtensions{.eps, .jpg}
-\fi
 
 \newcommand{\C}{\mathbb{C}}
 \newcommand{\R}{\mathbb{R}}
@@ -831,7 +809,7 @@ now lack of computational efficiency due to the generic implementation.
 
 % bibtex bibliography
 \bibliographystyle{plain}
-\bibliography{istl.bib}
+\bibliography{istl}
 
 % some links
 % http://www.netlib.org/blas/blast-forum/

dune/istl/CMakeLists.txt

@@ -10,6 +10,7 @@ install(FILES
    bdmatrix.hh
    btdmatrix.hh
    bvector.hh
+   cholmod.hh
    colcompmatrix.hh
    gsetc.hh
    ildl.hh

dune/istl/bcrsmatrix.hh

@@ -20,6 +20,8 @@
 #include <dune/common/iteratorfacades.hh>
 #include <dune/common/typetraits.hh>
 #include <dune/common/ftraits.hh>
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
 
 /*! \file
  * \brief Implementation of the BCRSMatrix class
@@ -442,7 +444,7 @@ namespace Dune {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -460,10 +462,7 @@ namespace Dune {
     typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
 
     //! increment block level counter
-    enum {
-      //! The number of blocklevels the matrix contains.
-      blocklevel = B::blocklevel+1
-    };
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     //! we support two modes
     enum BuildMode {
@@ -1579,12 +1578,16 @@ namespace Dune {
                                  "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         y[i.index()]=0;
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umv(x[j.index()],y[i.index()]);
+        {
+          auto&& xj = Impl::asVector(x[j.index()]);
+          auto&& yi = Impl::asVector(y[i.index()]);
+          Impl::asMatrix(*j).umv(xj, yi);
+        }
       }
     }
 
@@ -1599,11 +1602,15 @@ namespace Dune {
       if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umv(x[j.index()],y[i.index()]);
+        {
+          auto&& xj = Impl::asVector(x[j.index()]);
+          auto&& yi = Impl::asVector(y[i.index()]);
+          Impl::asMatrix(*j).umv(xj,yi);
+        }
       }
     }
 
@@ -1618,11 +1625,15 @@ namespace Dune {
       if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmv(x[j.index()],y[i.index()]);
+        {
+          auto&& xj = Impl::asVector(x[j.index()]);
+          auto&& yi = Impl::asVector(y[i.index()]);
+          Impl::asMatrix(*j).mmv(xj,yi);
+        }
       }
     }
 
@@ -1637,11 +1648,15 @@ namespace Dune {
       if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).usmv(alpha,x[j.index()],y[i.index()]);
+        {
+          auto&& xj = Impl::asVector(x[j.index()]);
+          auto&& yi = Impl::asVector(y[i.index()]);
+          Impl::asMatrix(*j).usmv(alpha,xj,yi);
+        }
       }
     }
 
@@ -1671,11 +1686,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umtv(x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).umtv(xi,yj);
+        }
       }
     }
 
@@ -1688,11 +1707,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmtv(x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).mmtv(xi,yj);
+        }
       }
     }
 
@@ -1707,11 +1730,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).usmtv(alpha,x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).usmtv(alpha,xi,yj);
+        }
       }
     }
 
@@ -1726,11 +1753,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umhv(x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).umhv(xi,yj);
+        }
       }
     }
 
@@ -1745,11 +1776,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmhv(x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).mmhv(xi,yj);
+        }
       }
     }
 
@@ -1764,11 +1799,15 @@ namespace Dune {
       if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
 #endif
       ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
+      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
       {
         ConstColIterator endj = (*i).end();
         for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).usmhv(alpha,x[i.index()],y[j.index()]);
+        {
+          auto&& xi = Impl::asVector(x[i.index()]);
+          auto&& yj = Impl::asVector(y[j.index()]);
+          Impl::asMatrix(*j).usmhv(alpha,xi,yj);
+        }
       }
     }
 
@@ -1783,15 +1822,11 @@ namespace Dune {
         DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
 #endif
 
-      double sum=0;
+      typename FieldTraits<field_type>::real_type sum=0;
 
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          sum += (*j).frobenius_norm2();
-      }
+      for (auto&& row : (*this))
+        for (auto&& entry : row)
+          sum += Impl::asMatrix(entry).frobenius_norm2();
 
       return sum;
     }
@@ -1816,7 +1851,7 @@ namespace Dune {
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm();
+          sum += Impl::asMatrix(y).infinity_norm();
         norm = max(sum, norm);
       }
       return norm;
@@ -1836,7 +1871,7 @@ namespace Dune {
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm_real();
+          sum += Impl::asMatrix(y).infinity_norm_real();
         norm = max(sum, norm);
       }
       return norm;
@@ -1857,12 +1892,12 @@ namespace Dune {
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm();
+          sum += Impl::asMatrix(y).infinity_norm();
         norm = max(sum, norm);
         isNaN += sum;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+
+      return norm * (isNaN / isNaN);
     }
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
@@ -1877,15 +1912,16 @@ namespace Dune {
 
       real_type norm = 0;
       real_type isNaN = 1;
+
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm_real();
+          sum += Impl::asMatrix(y).infinity_norm_real();
         norm = max(sum, norm);
         isNaN += sum;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+
+      return norm * (isNaN / isNaN);
     }
 
     //===== sizes

dune/istl/bdmatrix.hh

@@ -6,6 +6,7 @@
 #include <memory>
 
 #include <dune/common/rangeutilities.hh>
+#include <dune/common/scalarmatrixview.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 
@@ -32,7 +33,7 @@ namespace Dune {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -47,7 +48,7 @@ namespace Dune {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {blocklevel = B::blocklevel+1};
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     /** \brief Default constructor */
     BDMatrix() : BCRSMatrix<B,A>() {}
@@ -106,10 +107,25 @@ namespace Dune {
       return *this;
     }
 
+    /** \brief Solve the system Ax=b in O(n) time
+     *
+     * \exception ISTLError if the matrix is singular
+     *
+     */
+    template <class V>
+    void solve (V& x, const V& rhs) const {
+      for (size_type i=0; i<this->N(); i++)
+      {
+        auto&& xv = Impl::asVector(x[i]);
+        auto&& rhsv = Impl::asVector(rhs[i]);
+        Impl::asMatrix((*this)[i][i]).solve(xv,rhsv);
+      }
+    }
+
     /** \brief Inverts the matrix */
     void invert() {
-      for (int i=0; i<this->N(); i++)
-        (*this)[i][i].invert();
+      for (size_type i=0; i<this->N(); i++)
+        Impl::asMatrix((*this)[i][i]).invert();
     }
 
   private:

dune/istl/btdmatrix.hh

@@ -4,6 +4,8 @@
 #define DUNE_ISTL_BTDMATRIX_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
 #include <dune/istl/bcrsmatrix.hh>
 
 /** \file
@@ -29,7 +31,7 @@ namespace Dune {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -44,7 +46,7 @@ namespace Dune {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {blocklevel = B::blocklevel+1};
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     /** \brief Default constructor */
     BTDMatrix() : BCRSMatrix<B,A>() {}
@@ -120,7 +122,9 @@ namespace Dune {
 
       // special handling for 1x1 matrices.  The generic algorithm doesn't work for them
       if (this->N()==1) {
-        (*this)[0][0].solve(x[0],rhs[0]);
+        auto&& x0 = Impl::asVector(x[0]);
+        auto&& rhs0 = Impl::asVector(rhs[0]);
+        Impl::asMatrix((*this)[0][0]).solve(x0, rhs0);
         return;
       }
 
@@ -132,37 +136,38 @@ namespace Dune {
 
       /* Modify the coefficients. */
       block_type a_00_inv = (*this)[0][0];
-      a_00_inv.invert();
+      Impl::asMatrix(a_00_inv).invert();
 
-      //c[0] /= (*this)[0][0];	/* Division by zero risk. */
-      block_type c_0_tmp = c[0];
-      FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]);
+      //c[0] /= (*this)[0][0]; /* Division by zero risk. */
+      block_type tmp = a_00_inv;
+      Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[0]));
+      c[0] = tmp;
 
       // d = a^{-1} d        /* Division by zero would imply a singular matrix. */
-      typename V::block_type d_0_tmp = d[0];
-      a_00_inv.mv(d_0_tmp,d[0]);
+      auto d_0_tmp = d[0];
+      auto&& d_0 = Impl::asVector(d[0]);
+      Impl::asMatrix(a_00_inv).mv(Impl::asVector(d_0_tmp),d_0);
 
       for (unsigned int i = 1; i < this->N(); i++) {
 
         // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
         block_type tmp;
-        FMatrixHelp::multMatrix((*this)[i][i-1],c[i-1], tmp);
+        tmp = (*this)[i][i-1];
+        Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[i-1]));
+
         block_type id = (*this)[i][i];
         id -= tmp;
-        id.invert();         /* Division by zero risk. */
+        Impl::asMatrix(id).invert(); /* Division by zero risk. */
 
         if (i<c.size()) {
-          // c[i] *= id
-          tmp = c[i];
-          FMatrixHelp::multMatrix(id,tmp, c[i]);                      /* Last value calculated is redundant. */
+          Impl::asMatrix(c[i]).leftmultiply(Impl::asMatrix(id));            /* Last value calculated is redundant. */
         }
 
         // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
-        (*this)[i][i-1].mmv(d[i-1], d[i]);
-        typename V::block_type tmpVec = d[i];
-        id.mv(tmpVec, d[i]);
-        //d[i] *= id;
-
+        auto&& d_i = Impl::asVector(d[i]);
+        Impl::asMatrix((*this)[i][i-1]).mmv(Impl::asVector(d[i-1]), d_i);
+        auto tmpVec = d[i];
+        Impl::asMatrix(id).mv(Impl::asVector(tmpVec), d_i);
       }
 
       /* Now back substitute. */
@@ -170,7 +175,8 @@ namespace Dune {
       for (int i = this->N() - 2; i >= 0; i--) {
         //x[i] = d[i] - c[i] * x[i + 1];
         x[i] = d[i];
-        c[i].mmv(x[i+1], x[i]);
+        auto&& x_i = Impl::asVector(x[i]);
+        Impl::asMatrix(c[i]).mmv(Impl::asVector(x[i+1]), x_i);
       }
 
     }

dune/istl/bvector.hh

@@ -16,9 +16,12 @@
 #include <dune/common/deprecated.hh>
 #include <dune/common/dotproduct.hh>
 #include <dune/common/ftraits.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
 #include <dune/common/promotiontraits.hh>
 #include <dune/common/typetraits.hh>
 #include <dune/common/unused.hh>
+#include <dune/common/scalarvectorview.hh>
 
 #include "basearray.hh"
 #include "istlexception.hh"
@@ -35,6 +38,43 @@ namespace Dune {
 /** \brief Everything in this namespace is internal to dune-istl, and may change without warning */
 namespace Imp {
 
+  /** \brief Define some derived types transparently for number types and dune-istl vector types
+   *
+   * This is the actual implementation.  Calling code should use BlockTraits instead.
+   * \tparam isNumber Whether B is a number type (true) or a dune-istl matrix or vector type (false)
+   */
+  template <class B, bool isNumber>
+  class BlockTraitsImp;
+
+  template <class B>
+  class BlockTraitsImp<B,true>
+  {
+  public:
+    using field_type = B;
+
+    static constexpr unsigned int blockLevel()
+    {
+      return 0;
+    }
+  };
+
+  template <class B>
+  class BlockTraitsImp<B,false>
+  {
+  public:
+    using field_type = typename B::field_type;
+
+    static constexpr unsigned int blockLevel()
+    {
+      return B::blocklevel;
+    }
+  };
+
+  /** \brief Define some derived types transparently for number types and dune-istl matrix/vector types
+   */
+  template <class B>
+  using BlockTraits = BlockTraitsImp<B,IsNumber<B>::value>;
+
   /**
       \brief An unmanaged vector of blocks.
 
@@ -54,9 +94,7 @@ namespace Imp {
   public:
 
     //===== type definitions and constants
-
-    //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -133,7 +171,9 @@ namespace Imp {
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
 #endif
-      for (size_type i=0; i<this->n; ++i) (*this)[i].axpy(a,y[i]);
+      for (size_type i=0; i<this->n; ++i)
+        Impl::asVector((*this)[i]).axpy(a,Impl::asVector(y[i]));
+
       return *this;
     }
 
@@ -146,9 +186,9 @@ namespace Imp {
      * @return
      */
     template<class OtherB, class OtherA>
-    typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
+    auto operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
     {
-      typedef typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType PromotedType;
+      typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
       PromotedType sum(0);
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
@@ -167,14 +207,17 @@ namespace Imp {
      * @return
      */
     template<class OtherB, class OtherA>
-    typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
+    auto dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
     {
-      typedef typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType PromotedType;
+      typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
       PromotedType sum(0);
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
 #endif
-      for (size_type i=0; i<this->n; ++i) sum += ((*this)[i]).dot(y[i]);
+
+      for (size_type i=0; i<this->n; ++i)
+        sum += Impl::asVector((*this)[i]).dot(Impl::asVector(y[i]));
+
       return sum;
     }
 
@@ -184,7 +227,8 @@ namespace Imp {
     typename FieldTraits<field_type>::real_type one_norm () const
     {
       typename FieldTraits<field_type>::real_type sum=0;
-      for (size_type i=0; i<this->n; ++i) sum += (*this)[i].one_norm();
+      for (size_type i=0; i<this->n; ++i)
+        sum += Impl::asVector((*this)[i]).one_norm();
       return sum;
     }
 
@@ -192,7 +236,8 @@ namespace Imp {
     typename FieldTraits<field_type>::real_type one_norm_real () const
     {
       typename FieldTraits<field_type>::real_type sum=0;
-      for (size_type i=0; i<this->n; ++i) sum += (*this)[i].one_norm_real();
+      for (size_type i=0; i<this->n; ++i)
+        sum += Impl::asVector((*this)[i]).one_norm_real();
       return sum;
     }
 
@@ -200,16 +245,15 @@ namespace Imp {
     typename FieldTraits<field_type>::real_type two_norm () const
     {
       using std::sqrt;
-      typename FieldTraits<field_type>::real_type sum=0;
-      for (size_type i=0; i<this->n; ++i) sum += (*this)[i].two_norm2();
-      return sqrt(sum);
+      return sqrt(two_norm2());
     }
 
     //! Square of the two-norm (the sum over the squared values of the entries)
     typename FieldTraits<field_type>::real_type two_norm2 () const
     {
       typename FieldTraits<field_type>::real_type sum=0;
-      for (size_type i=0; i<this->n; ++i) sum += (*this)[i].two_norm2();
+      for (size_type i=0; i<this->n; ++i)
+        sum += Impl::asVector((*this)[i]).two_norm2();
       return sum;
     }
 
@@ -221,8 +265,8 @@ namespace Imp {
       using std::max;
 
       real_type norm = 0;
-      for (auto const &x : *this) {
-        real_type const a = x.infinity_norm();
+      for (auto const &xi : *this) {
+        real_type const a = Impl::asVector(xi).infinity_norm();
         norm = max(a, norm);
       }
       return norm;
@@ -236,8 +280,8 @@ namespace Imp {
       using std::max;
 
       real_type norm = 0;
-      for (auto const &x : *this) {
-        real_type const a = x.infinity_norm_real();
+      for (auto const &xi : *this) {
+        real_type const a = Impl::asVector(xi).infinity_norm_real();
         norm = max(a, norm);
       }
       return norm;
@@ -249,16 +293,17 @@ namespace Imp {
     typename FieldTraits<ft>::real_type infinity_norm() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
+      using std::abs;
 
       real_type norm = 0;
       real_type isNaN = 1;
-      for (auto const &x : *this) {
-        real_type const a = x.infinity_norm();
+
+      for (auto const &xi : *this) {
+        real_type const a = Impl::asVector(xi).infinity_norm();
         norm = max(a, norm);
         isNaN += a;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+      return norm * (isNaN / isNaN);
     }
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
@@ -270,13 +315,14 @@ namespace Imp {
 
       real_type norm = 0;
       real_type isNaN = 1;
-      for (auto const &x : *this) {
-        real_type const a = x.infinity_norm_real();
+
+      for (auto const &xi : *this) {
+        real_type const a = Impl::asVector(xi).infinity_norm_real();
         norm = max(a, norm);
         isNaN += a;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+
+      return norm * (isNaN / isNaN);
     }
 
     //===== sizes
@@ -291,8 +337,10 @@ namespace Imp {
     size_type dim () const
     {
       size_type d=0;
+
       for (size_type i=0; i<this->n; i++)
-        d += (*this)[i].dim();
+        d += Impl::asVector((*this)[i]).dim();
+
       return d;
     }
 
@@ -358,7 +406,7 @@ namespace Imp {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -370,10 +418,7 @@ namespace Imp {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {
-      //! The number of blocklevel we contain.
-      blocklevel = B::blocklevel+1
-    };
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     //! make iterators available as types
     typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
@@ -418,7 +463,7 @@ namespace Imp {
     {
       static_assert(std::numeric_limits<S>::is_integer,
         "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
-      if(_capacity > _n)
+      if((size_type)_capacity > _n)
         storage_.reserve(_capacity);
       storage_.resize(_n);
       syncBaseArray();
@@ -645,7 +690,7 @@ namespace Imp {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -657,10 +702,7 @@ namespace Imp {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {
-      //! The number of blocklevels we contain
-      blocklevel = B::blocklevel+1
-    };
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     //! make iterators available as types
     typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
@@ -774,7 +816,7 @@ namespace Imp {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -832,7 +874,8 @@ namespace Imp {
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
 #endif
-      for (size_type i=0; i<y.n; ++i) (this->operator[](y.j[i])).axpy(a,y.p[i]);
+      for (size_type i=0; i<y.n; ++i)
+        Impl::asVector((*this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
       return *this;
     }
 
@@ -946,8 +989,7 @@ namespace Imp {
         norm = max(a, norm);
         isNaN += a;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+      return norm * (isNaN / isNaN);
     }
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
@@ -964,8 +1006,7 @@ namespace Imp {
         norm = max(a, norm);
         isNaN += a;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+      return norm * (isNaN / isNaN);
     }
 
     //===== sizes
@@ -1029,7 +1070,7 @@ namespace Imp {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -1041,10 +1082,7 @@ namespace Imp {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {
-      //! The number of block level this vector contains.
-      blocklevel = B::blocklevel+1
-    };
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     //! make iterators available as types
     typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;

dune/istl/cholmod.hh

+#pragma once
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include<dune/istl/solver.hh>
+
+
+#include <memory>
+
+#include <cholmod.h>
+
+namespace Dune {
+
+namespace Impl{
+
+  /** @brief Dummy class for empty ignore nodes
+   *
+   * This class implements "no" ignore nodes with a
+   * compatible interface for the "setMatrix" method.
+   *
+   * It should be optimized out by the compiler, so no
+   * overhead should be preduced
+   */
+  struct NoIgnore
+  {
+    const NoIgnore& operator[](std::size_t) const { return *this; }
+    explicit operator bool() const { return false; }
+    std::size_t count() const { return 0; }
+  };
+
+} //namespace Impl
+
+template<class T>
+class Cholmod;
+
+
+template<class T, class A, int k>
+class Cholmod<BCRSMatrix<FieldMatrix<T,k,k>, A>>
+  : public InverseOperator<
+      BlockVector<FieldVector<T,k>,
+        typename A::template rebind<FieldVector<T,k>>::other>,
+      BlockVector<FieldVector<T,k>,
+        typename A::template rebind<FieldVector<T,k>>::other>>
+{
+public:
+
+  // type of unknown
+  using X = BlockVector<FieldVector<T,k>,
+    typename A::template rebind<FieldVector<T,k>>::other>;
+  // type of rhs
+  using B = BlockVector<FieldVector<T,k>,
+    typename A::template rebind<FieldVector<T,k>>::other>;
+  // type of external matrix
+  using Matrix = BCRSMatrix<FieldMatrix<T,k,k>, A>;
+
+
+  /** @brief Default constructor
+   *
+   *  Calls the the cholmod runtime,
+   *  see CHOLMOD doc
+   */
+  Cholmod()
+  {
+    cholmod_start(&c_);
+  }
+
+  /** @brief Destructor
+   *
+   *  Free space and calls cholmod_finish,
+   *  see CHOLMOD doc
+   */
+  ~Cholmod()
+  {
+    if (L_)
+      cholmod_free_factor(&L_, &c_);
+    cholmod_finish(&c_);
+  }
+
+  // forbit this to avoid freeing memory twice
+  Cholmod(const Cholmod&) = delete;
+  Cholmod& operator=(const Cholmod&) = delete;
+
+
+  /**
+    *  \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
+    */
+  void apply (X& x, B& b, double reduction, InverseOperatorResult& res)
+  {
+    DUNE_UNUSED_PARAMETER(reduction);
+    apply(x,b,res);
+  }
+
+  /**
+  *  \copydoc InverseOperator::apply(X&, Y&, InverseOperatorResult&)
+  */
+  void apply(X& x, B& b, InverseOperatorResult& res)
+  {
+    // cast to double array
+    auto b2 = std::make_unique<double[]>(L_->n);
+    auto x2 = std::make_unique<double[]>(L_->n);
+
+    // copy to cholmod
+    auto bp = b2.get();
+    for (const auto& bi : b)
+      for (const auto& bii : bi)
+        *bp++ = bii;
+
+    // create a cholmod dense object
+    auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
+    // cast because void-ptr
+    auto b4 = static_cast<double*>(b3->x);
+    std::copy(b2.get(), b2.get() + L_->n, b4);
+
+    // solve for a cholmod x object
+    auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
+    // cast because void-ptr
+    auto xp = static_cast<double*>(x3->x);
+
+    // copy into x
+    x.resize(L_->n);
+    for (auto& xi : x)
+      for (auto& xii : xi)
+        xii = *xp++;
+
+    // statistics for a direct solver
+    res.iterations = 1;
+    res.converged = true;
+  }
+
+
+  /** @brief Set matrix without ignore nodes
+   *
+   * This method forwards a nullptr to the setMatrix method
+   * with ignore nodes
+   */
+  void setMatrix(const Matrix& matrix)
+  {
+    const Impl::NoIgnore* noIgnore = nullptr;
+    setMatrix(matrix, noIgnore);
+  }
+
+  /** @brief Set matrix and ignore nodes
+   *
+   * The input matrix is copied to CHOLMOD compatible form.
+   * The ignore argument consists of a compatible BitVector,
+   * indicating the dof's which has to be deleted from the matrix
+   * Decomposing the matrix at the end takes a lot of time
+   * \param [in] matrix Matrix to be decomposed. In BCRS compatible form
+   * \param [in] ignore Pointer to a compatible BitVector
+   */
+  template<class Ignore>
+  void setMatrix(const Matrix& matrix, const Ignore* ignore)
+  {
+
+    // Total number of rows
+    int N = k * matrix.N();
+    if ( ignore )
+      N -= ignore->count();
+
+    // number of nonzeroes
+    const int nnz = k * k * matrix.nonzeroes();
+    // number of diagonal entries
+    const int nDiag = k * k * matrix.N();
+    // number of nonzeroes in the dialgonal
+    const int nnzDiag = (k * (k+1)) / 2 * matrix.N();
+    // number of nonzeroes in triangular submatrix (including diagonal)
+    const int nnzTri = (nnz - nDiag) / 2 + nnzDiag;
+
+    /*
+    * CHOLMOD uses compressed-column sparse matrices, but for symmetric
+    * matrices this is the same as the compressed-row sparse matrix used
+    * by DUNE.  So we can just store Mᵀ instead of M (as M = Mᵀ).
+    */
+    const auto deleter = [c = &this->c_](auto* p) {
+      cholmod_free_sparse(&p, c);
+    };
+    auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
+      cholmod_allocate_sparse(N,             // # rows
+                              N,             // # cols
+                              nnzTri,        // # of nonzeroes
+                              1,             // indices are sorted ( 1 = true)
+                              1,             // matrix is "packed" ( 1 = true)
+                              -1,            // stype of matrix ( -1 = cosider the lower part only )
+                              CHOLMOD_REAL,  // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
+                              &c_            // cholmod_common ptr
+                             ), deleter);
+
+    // copy the data of BCRS matrix to Cholmod Sparse matrix
+    int* Ap = static_cast<int*>(M->p);
+    int* Ai = static_cast<int*>(M->i);
+    double* Ax = static_cast<double*>(M->x);
+
+    std::size_t blocksize = k;
+
+    // create a vector that maps each remaining matrix index to it's number in the condensed matrix
+    std::vector<std::size_t> subIndices;
+
+    if ( ignore )
+    {
+      subIndices.resize(matrix.M()*blocksize);
+
+      std::size_t j=0;
+      for( std::size_t block=0; block<matrix.N(); block++ )
+      {
+        for( std::size_t i=0; i<blocksize; i++ )
+        {
+          if( not (*ignore)[block][i] )
+          {
+            subIndices[ block*blocksize + i ] = j++;
+          }
+        }
+      }
+    }
+
+    // Copy the data to the CHOLMOD array
+    int n = 0;
+    for (auto rowIt = matrix.begin(); rowIt != matrix.end(); rowIt++)
+    {
+      const auto row = rowIt.index();
+      for (std::size_t i=0; i<blocksize; i++)
+      {
+        if( ignore and (*ignore)[row][i] )
+          continue;
+
+        // col start
+        *Ap++ = n;
+
+        for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+        {
+          const auto col = colIt.index();
+
+          // are we already in the lower part?
+          if (col < row)
+            continue;
+
+          for (auto j = (row == col ? i : 0); j<blocksize; j++)
+          {
+            if( ignore and (*ignore)[col][j] )
+              continue;
+
+            const auto jj = ignore ? subIndices[  blocksize*col + j ] : blocksize*col + j;
+
+            // set the current index and entry
+            *Ai++ = jj;
+            *Ax++ = (*colIt)[i][j];
+            // increase number of set values
+            n++;
+          }
+        }
+      }
+    }
+    // set last col start
+    *Ap = n;
+
+    // Now analyse the pattern and optimal row order
+    L_ = cholmod_analyze(M.get(), &c_);
+
+    // Do the factorization (this may take some time)
+    cholmod_factorize(M.get(), L_, &c_);
+  }
+
+  virtual SolverCategory::Category category() const
+  {
+    return SolverCategory::Category::sequential;
+  }
+
+private:
+
+  // create a destrucable unique_ptr
+  auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
+  {
+    const auto deleter = [c](auto* p) {
+      cholmod_free_dense(&p, c);
+    };
+    return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
+  }
+
+  cholmod_common c_;
+  cholmod_factor* L_ = nullptr;
+};
+
+
+} /* namespace Dune */

dune/istl/colcompmatrix.hh

@@ -8,6 +8,7 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/typetraits.hh>
 #include <dune/common/unused.hh>
+#include <dune/common/scalarmatrixview.hh>
 #include <limits>
 
 namespace Dune
@@ -136,22 +137,14 @@ namespace Dune
     const RowIndexSet& s_;
   };
 
-  /**
-   * @brief Utility class for converting an ISTL Matrix into a column-compressed matrix
-   * @tparam M the matrix type
-   */
-  template<class M>
-  struct ColCompMatrix
-  {};
-
   /**
    * @brief Inititializer for the ColCompMatrix
    * as needed by OverlappingSchwarz
    * @tparam M the matrix type
+   * @tparam I the internal index type
    */
-  template<class M>
-  struct ColCompMatrixInitializer
-  {};
+  template<class M, class I = int>
+  class ColCompMatrixInitializer;
 
   template<class M, class X, class TM, class TD, class T1>
   class SeqOverlappingSchwarz;
@@ -160,21 +153,27 @@ namespace Dune
   struct SeqOverlappingSchwarzAssemblerHelper;
 
   /**
-   * @brief Converter for BCRSMatrix to column-compressed Matrix.
-   * specialization for BCRSMatrix
+   * @brief Utility class for converting an ISTL Matrix into a column-compressed matrix
+   * @tparam M The matrix type
+   * @tparam I the internal index type
    */
-  template<class B, class TA, int n, int m>
-  class ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
+  template<class Mat, class I = int>
+  class ColCompMatrix
   {
-    friend struct ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
+    friend class ColCompMatrixInitializer<Mat, I>;
+
+    using B = typename Mat::field_type;
 
   public:
     /** @brief The type of the matrix to convert. */
-    typedef BCRSMatrix<FieldMatrix<B,n,m>,TA> Matrix;
+    using Matrix = Mat;
+
     friend struct SeqOverlappingSchwarzAssemblerHelper<ColCompMatrix<Matrix>, true>;
 
     typedef typename Matrix::size_type size_type;
 
+    using Index = I;
+
     /**
      * @brief Constructor that initializes the data.
      * @param mat The matrix to convert.
@@ -197,6 +196,10 @@ namespace Dune
 
     size_type nnz() const
     {
+      // TODO: The following code assumes that the blocks are dense
+      // and that they all have the same dimensions.
+      const auto n = MatrixDimension<typename Matrix::block_type>::rowdim();
+      const auto m = MatrixDimension<typename Matrix::block_type>::coldim();
       return Nnz_/n/m;
     }
 
@@ -214,12 +217,12 @@ namespace Dune
       return values;
     }
 
-    int* getRowIndex() const
+    Index* getRowIndex() const
     {
       return rowindex;
     }
 
-    int* getColStart() const
+    Index* getColStart() const
     {
       return colstart;
     }
@@ -240,34 +243,50 @@ namespace Dune
     virtual void setMatrix(const Matrix& mat);
 
   public:
-    int N_, M_, Nnz_;
+    size_type N_, M_, Nnz_;
     B* values;
-    int* rowindex;
-    int* colstart;
+    Index* rowindex;
+    Index* colstart;
   };
 
-  template<class T, class A, int n, int m>
-  class ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >
+  template<class M, class I>
+  class ColCompMatrixInitializer
   {
-    template<class I, class S, class D>
+    template<class IList, class S, class D>
     friend class OverlappingSchwarzInitializer;
   public:
-    typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
-    typedef Dune::ColCompMatrix<Matrix> ColCompMatrix;
+    using Matrix = M;
+    using Index = I;
+    typedef Dune::ColCompMatrix<Matrix, Index> ColCompMatrix;
     typedef typename Matrix::row_type::const_iterator CIter;
     typedef typename Matrix::size_type size_type;
 
-    ColCompMatrixInitializer(ColCompMatrix& lum);
+    /** \brief Constructor for scalar-valued matrices
+     *
+     * \tparam Block Dummy parameter to make SFINAE work
+     */
+    template <class Block=typename M::block_type>
+    ColCompMatrixInitializer(ColCompMatrix& lum,
+                             typename std::enable_if_t<Dune::IsNumber<Block>::value>* sfinae = nullptr);
 
-    ColCompMatrixInitializer();
+    /** \brief Constructor for dense matrix-valued matrices
+     *
+     * \tparam Block Dummy parameter to make SFINAE work
+     */
+    template <class Block=typename M::block_type>
+    ColCompMatrixInitializer(ColCompMatrix& lum,
+                             typename std::enable_if_t<!Dune::IsNumber<Block>::value>* sfinae = nullptr);
 
-    virtual ~ColCompMatrixInitializer();
+    ColCompMatrixInitializer();
 
     template<typename Iter>
     void addRowNnz(const Iter& row) const;
 
-    template<typename Iter, typename Set>
-    void addRowNnz(const Iter& row, const Set& s) const;
+    template<typename Iter, typename FullMatrixIndex>
+    void addRowNnz(const Iter& row, const std::set<FullMatrixIndex>& indices) const;
+
+    template<typename Iter, typename SubMatrixIndex>
+    void addRowNnz(const Iter& row, const std::vector<SubMatrixIndex>& indices) const;
 
     void allocate();
 
@@ -293,42 +312,56 @@ namespace Dune
 
     ColCompMatrix* mat;
     size_type cols;
-    mutable size_type *marker;
+
+    // Number of rows/columns of the matrix entries
+    // (assumed to be scalars or dense matrices)
+    size_type n, m;
+
+    mutable std::vector<size_type> marker;
   };
 
-  template<class T, class A, int n, int m>
-  ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer(ColCompMatrix& mat_)
-    : mat(&mat_), cols(mat_.M()), marker(0)
+  template<class M, class I>
+  template <class Block>
+  ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer(ColCompMatrix& mat_,typename std::enable_if_t<Dune::IsNumber<Block>::value>* sfinae)
+    : mat(&mat_), cols(mat_.M())
   {
+    n = 1;
+    m = 1;
+
     mat->Nnz_=0;
   }
 
-  template<class T, class A, int n, int m>
-  ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer()
-    : mat(0), cols(0), marker(0)
-  {}
-
-  template<class T, class A, int n, int m>
-  ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::~ColCompMatrixInitializer()
+  template<class M, class I>
+  template <class Block>
+  ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer(ColCompMatrix& mat_,typename std::enable_if_t<!Dune::IsNumber<Block>::value>* sfinae)
+    : mat(&mat_), cols(mat_.M())
   {
-    if(marker)
-      delete[] marker;
+    // WARNING: This assumes that all blocks are dense and identical
+    n = M::block_type::rows;
+    m = M::block_type::cols;
+
+    mat->Nnz_=0;
   }
 
-  template<class T, class A, int n, int m>
+  template<class M, class I>
+  ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer()
+    : mat(0), cols(0), n(0), m(0)
+  {}
+
+  template<class M, class I>
   template<typename Iter>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row) const
+  void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row) const
   {
     mat->Nnz_+=row->getsize();
   }
 
-  template<class T, class A, int n, int m>
-  template<typename Iter, typename Map>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
-                                                                            const Map& indices) const
+  template<class M, class I>
+  template<typename Iter, typename FullMatrixIndex>
+  void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row,
+                                                                            const std::set<FullMatrixIndex>& indices) const
   {
     typedef typename  Iter::value_type::const_iterator RIter;
-    typedef typename Map::const_iterator MIter;
+    typedef typename std::set<FullMatrixIndex>::const_iterator MIter;
     MIter siter =indices.begin();
     for(RIter entry=row->begin(); entry!=row->end(); ++entry)
     {
@@ -341,42 +374,51 @@ namespace Dune
     }
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocate()
+  template<class M, class I>
+  template<typename Iter, typename SubMatrixIndex>
+  void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row,
+                                                                            const std::vector<SubMatrixIndex>& indices) const
+  {
+    using RIter = typename Iter::value_type::const_iterator;
+    for(RIter entry=row->begin(); entry!=row->end(); ++entry)
+      if (indices[entry.index()]!=std::numeric_limits<SubMatrixIndex>::max())
+          ++mat->Nnz_;
+  }
+
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::allocate()
   {
     allocateMatrixStorage();
     allocateMarker();
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage() const
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::allocateMatrixStorage() const
   {
     mat->Nnz_*=n*m;
     // initialize data
-    mat->values=new T[mat->Nnz_];
-    mat->rowindex=new int[mat->Nnz_];
-    mat->colstart=new int[cols+1];
+    mat->values=new typename M::field_type[mat->Nnz_];
+    mat->rowindex=new I[mat->Nnz_];
+    mat->colstart=new I[cols+1];
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMarker()
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::allocateMarker()
   {
-    marker = new typename Matrix::size_type[cols];
-
-    for(size_type i=0; i < cols; ++i)
-      marker[i]=0;
+    marker.resize(cols);
+    std::fill(marker.begin(), marker.end(), 0);
   }
 
-  template<class T, class A, int n, int m>
+  template<class M, class I>
   template<typename Iter>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(const Iter& row, const CIter& col) const
+  void ColCompMatrixInitializer<M, I>::countEntries(const Iter& row, const CIter& col) const
   {
     DUNE_UNUSED_PARAMETER(row);
     countEntries(col.index());
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(size_type colindex) const
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::countEntries(size_type colindex) const
   {
     for(size_type i=0; i < m; ++i)
     {
@@ -385,8 +427,8 @@ namespace Dune
     }
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::calcColstart() const
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::calcColstart() const
   {
     mat->colstart[0]=0;
     for(size_type i=0; i < cols; ++i) {
@@ -396,32 +438,31 @@ namespace Dune
     }
   }
 
-  template<class T, class A, int n, int m>
+  template<class M, class I>
   template<typename Iter>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col) const
+  void ColCompMatrixInitializer<M, I>::copyValue(const Iter& row, const CIter& col) const
   {
     copyValue(col, row.index(), col.index());
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const CIter& col, size_type rowindex, size_type colindex) const
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::copyValue(const CIter& col, size_type rowindex, size_type colindex) const
   {
     for(size_type i=0; i<n; i++) {
       for(size_type j=0; j<m; j++) {
-        assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
-        assert((int)marker[colindex*m+j]<mat->Nnz_);
+        assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
+        assert((size_type)marker[colindex*m+j]<mat->Nnz_);
         mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
-        mat->values[marker[colindex*m+j]]=(*col)[i][j];
+        mat->values[marker[colindex*m+j]]=Impl::asMatrix(*col)[i][j];
         ++marker[colindex*m+j]; // index for next entry in column
       }
     }
   }
 
-  template<class T, class A, int n, int m>
-  void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::createMatrix() const
+  template<class M, class I>
+  void ColCompMatrixInitializer<M, I>::createMatrix() const
   {
-    delete[] marker;
-    marker=0;
+    marker.clear();
   }
 
   template<class F, class MRS>
@@ -461,12 +502,6 @@ namespace Dune
     typedef typename std::iterator_traits<Iter>::value_type row_type;
     typedef typename row_type::const_iterator CIter;
 
-    // Calculate upper Bound for nonzeros
-    for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
-      initializer.addRowNnz(row, mrs.rowIndexSet());
-
-    initializer.allocate();
-
     typedef typename MRS::Matrix::size_type size_type;
 
     // A vector containing the corresponding indices in
@@ -479,6 +514,12 @@ namespace Dune
     for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
       subMatrixIndex[*index]=s++;
 
+    // Calculate upper Bound for nonzeros
+    for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
+      initializer.addRowNnz(row, subMatrixIndex);
+
+    initializer.allocate();
+
     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
       for(CIter col=row->begin(); col != row->end(); ++col) {
         if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
@@ -499,20 +540,26 @@ namespace Dune
 
 #ifndef DOXYGEN
 
-  template<class B, class TA, int n, int m>
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::ColCompMatrix()
+  template<class Mat, class I>
+  ColCompMatrix<Mat, I>::ColCompMatrix()
     : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
   {}
 
-  template<class B, class TA, int n, int m>
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
+  template<class Mat, class I>
+  ColCompMatrix<Mat, I>
   ::ColCompMatrix(const Matrix& mat)
-    : N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.nonzeroes())
-  {}
+  {
+    // WARNING: This assumes that all blocks are dense and identical
+    const size_type n = MatrixDimension<typename Mat::block_type>::rowdim();
+    const size_type m = MatrixDimension<typename Mat::block_type>::coldim();
+    N_ = n*mat.N();
+    M_ = m*mat.N();
+    Nnz_ = n*m*mat.nonzeroes();
+  }
 
-  template<class B, class TA, int n, int m>
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
+  template<class Mat, class I>
+  ColCompMatrix<Mat, I>&
+  ColCompMatrix<Mat, I>::operator=(const Matrix& mat)
   {
     if(N_+M_+Nnz_!=0)
       free();
@@ -520,9 +567,9 @@ namespace Dune
     return *this;
   }
 
-  template<class B, class TA, int n, int m>
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const ColCompMatrix& mat)
+  template<class Mat, class I>
+  ColCompMatrix<Mat, I>&
+  ColCompMatrix<Mat, I>::operator=(const ColCompMatrix& mat)
   {
     if(N_+M_+Nnz_!=0)
       free();
@@ -530,57 +577,58 @@ namespace Dune
     M_=mat.M_;
     Nnz_= mat.Nnz_;
     if(M_>0) {
-      colstart=new int[M_+1];
-      for(int i=0; i<=M_; ++i)
+      colstart=new size_type[M_+1];
+      for(size_type i=0; i<=M_; ++i)
         colstart[i]=mat.colstart[i];
     }
 
     if(Nnz_>0) {
       values = new B[Nnz_];
-      rowindex = new int[Nnz_];
+      rowindex = new size_type[Nnz_];
 
-      for(int i=0; i<Nnz_; ++i)
+      for(size_type i=0; i<Nnz_; ++i)
         values[i]=mat.values[i];
 
-      for(int i=0; i<Nnz_; ++i)
+      for(size_type i=0; i<Nnz_; ++i)
         rowindex[i]=mat.rowindex[i];
     }
     return *this;
   }
 
-  template<class B, class TA, int n, int m>
-  void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
+  template<class Mat, class I>
+  void ColCompMatrix<Mat, I>
   ::setMatrix(const Matrix& mat)
   {
-    N_=n*mat.N();
-    M_=m*mat.M();
-    ColCompMatrixInitializer<Matrix> initializer(*this);
+    N_=MatrixDimension<Mat>::rowdim(mat);
+    M_=MatrixDimension<Mat>::coldim(mat);
+    ColCompMatrixInitializer<Mat, I> initializer(*this);
 
     copyToColCompMatrix(initializer, MatrixRowSet<Matrix>(mat));
   }
 
-  template<class B, class TA, int n, int m>
-  void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
+  template<class Mat, class I>
+  void ColCompMatrix<Mat, I>
   ::setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs)
   {
     if(N_+M_+Nnz_!=0)
       free();
-    N_=mrs.size()*n;
-    M_=mrs.size()*m;
-    ColCompMatrixInitializer<Matrix> initializer(*this);
 
-    copyToColCompMatrix(initializer, MatrixRowSubset<Matrix,std::set<std::size_t> >(mat,mrs));
+    N_=mrs.size()*MatrixDimension<Mat>::rowdim(mat) / mat.N();
+    M_=mrs.size()*MatrixDimension<Mat>::coldim(mat) / mat.M();
+    ColCompMatrixInitializer<Mat, I> initializer(*this);
+
+    copyToColCompMatrix(initializer, MatrixRowSubset<Mat,std::set<std::size_t> >(mat,mrs));
   }
 
-  template<class B, class TA, int n, int m>
-  ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~ColCompMatrix()
+  template<class Mat, class I>
+  ColCompMatrix<Mat, I>::~ColCompMatrix()
   {
     if(N_+M_+Nnz_!=0)
       free();
   }
 
-  template<class B, class TA, int n, int m>
-  void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
+  template<class Mat, class I>
+  void ColCompMatrix<Mat, I>::free()
   {
     delete[] values;
     delete[] rowindex;

dune/istl/gsetc.hh

@@ -8,6 +8,9 @@
 #include <iostream>
 #include <iomanip>
 #include <string>
+
+#include <dune/common/hybridutilities.hh>
+
 #include "multitypeblockvector.hh"
 #include "multitypeblockmatrix.hh"
 
@@ -385,12 +388,23 @@ namespace Dune {
         rhs = b[i.index()];           // rhs = b_i
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j++;           // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
+            for (; j != endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            x[i.index()] = rhs / id(*diag);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)           // iterate over a_ij with j < i
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
         coliterator diag=j++;           // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
         for (; j != endj; ++j)           // iterate over a_ij with j > i
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j>i} a_ij * xold_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j>i} a_ij * xold_j
         algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w);           // if I==1: xnew_i = rhs/a_ii
+        });
       }
       // next two lines: xnew_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j) + (1-w)*xold;
       x *= w;
@@ -417,13 +431,25 @@ namespace Dune {
         rhs = b[i.index()];           // rhs = b_i
         coliterator endj=(*i).end();           // iterate over a_ij with j < i
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];               //  rhs -= sum_{j<i} a_ij * xnew_j
+            coliterator diag=j;           // *diag = a_ii
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];               // rhs -= sum_{j<i} a_ij * xnew_j
+            v = rhs / id(*diag);
+            id(x[i.index()]) += w*id(v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);               //  rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               //  rhs -= sum_{j<i} a_ij * xnew_j
         coliterator diag=j;           // *diag = a_ii
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
         algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w);           // if blocksize I==1: v = rhs/a_ii
-        x[i.index()].axpy(w,v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+        id(x[i.index()]).axpy(w,v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+          });
       }
     }
 
@@ -447,13 +473,25 @@ namespace Dune {
         rhs = b[i.index()];
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j;
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            v = rhs / id(*diag);
+            x[i.index()] += w*id(v);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         coliterator diag=j;
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         algmeta_itsteps<I-1,typename M::block_type>::bsorb(*diag,v,rhs,w);
-        x[i.index()].axpy(w,v);
+        id(x[i.index()]).axpy(w,v);
+          });
       }
     }
 
@@ -473,12 +511,23 @@ namespace Dune {
         rhs = b[i.index()];
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j;
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            v[i.index()] = rhs / id(*diag);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         coliterator diag=j;
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
+          });
       }
       x.axpy(w,v);
     }

dune/istl/ildl.hh

 #ifndef DUNE_ISTL_ILDL_HH
 #define DUNE_ISTL_ILDL_HH
 
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
 #include "ilu.hh"
 
 /**
@@ -29,8 +31,16 @@ namespace Dune
     }
   }
 
+  template< class K >
+  inline static void bildl_subtractBCT ( const K &B, const K &CT, K &A,
+                                         typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr )
+  {
+    A -= B * CT;
+  }
+
   template< class Matrix >
-  inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matrix &A )
+  inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matrix &A,
+                                         typename std::enable_if_t<!Dune::IsNumber<Matrix>::value>* sfinae = nullptr )
   {
     for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
     {
@@ -114,12 +124,12 @@ namespace Dune
         const auto &A_k = A[ ik.index() ];
 
         auto B = A_ik;
-        A_ik.rightmultiply( *A_k.find( ik.index() ) );
+        Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
         bildl_subtractBCT( B, A_ik, A_ii );
       }
       try
       {
-        A_ii.invert();
+        Impl::asMatrix(A_ii).invert();
       }
       catch( const Dune::FMatrixError &e )
       {
@@ -142,7 +152,10 @@ namespace Dune
       const auto &A_i = *i;
       v[ i.index() ] = d[ i.index() ];
       for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
-        (*ij).mmv( v[ ij.index() ], v[ i.index() ] );
+      {
+        auto&& vi = Impl::asVector( v[ i.index() ] );
+        Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
+      }
     }
 
     // solve D w = v, note: diagonal stores Dii^{-1}
@@ -155,8 +168,16 @@ namespace Dune
         const auto &A_i = *i;
         const auto ii = A_i.beforeEnd();
         assert( ii.index() == i.index() );
-        auto rhs = v[ i.index() ];
-        ii->mv( rhs, v[ i.index() ] );
+        // We need to be careful here: Directly using
+        // auto rhs = Impl::asVector(v[ i.index() ]);
+        // is not OK in case this is a proxy. Hence
+        // we first have to copy the value. Notice that
+        // this is still not OK, if the vector type itself returns
+        // proxy references.
+        auto rhsValue = v[ i.index() ];
+        auto&& rhs = Impl::asVector(rhsValue);
+        auto&& vi = Impl::asVector( v[ i.index() ] );
+        Impl::asMatrix(*ii).mv(rhs, vi);
       }
     }
     else
@@ -168,8 +189,16 @@ namespace Dune
         const auto &A_i = *i;
         const auto ii = A_i.find( i.index() );
         assert( ii.index() == i.index() );
-        auto rhs = v[ i.index() ];
-        ii->mv( rhs, v[ i.index() ] );
+        // We need to be careful here: Directly using
+        // auto rhs = Impl::asVector(v[ i.index() ]);
+        // is not OK in case this is a proxy. Hence
+        // we first have to copy the value. Notice that
+        // this is still not OK, if the vector type itself returns
+        // proxy references.
+        auto rhsValue = v[ i.index() ];
+        auto&& rhs = Impl::asVector(rhsValue);
+        auto&& vi = Impl::asVector( v[ i.index() ] );
+        Impl::asMatrix(*ii).mv(rhs, vi);
       }
     }
 
@@ -179,7 +208,10 @@ namespace Dune
     {
       const auto &A_i = *i;
       for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
-        (*ij).mmtv( v[ i.index() ], v[ ij.index() ] );
+      {
+        auto&& vij = Impl::asVector( v[ ij.index() ] );
+        Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
+      }
     }
   }
 

dune/istl/ilu.hh

@@ -5,14 +5,12 @@
 
 #include <cmath>
 #include <complex>
-#include <iostream>
-#include <iomanip>
-#include <string>
-#include <set>
 #include <map>
 #include <vector>
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
 #include "istlexception.hh"
 
 /** \file
@@ -55,7 +53,7 @@ namespace Dune {
         coliterator jj = A[ij.index()].find(ij.index());
 
         // compute L_ij = A_jj^-1 * A_ij
-        (*ij).rightmultiply(*jj);
+        Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
 
         // modify row
         coliterator endjk=A[ij.index()].end();                 // end of row j
@@ -65,7 +63,7 @@ namespace Dune {
           if (ik.index()==jk.index())
           {
             block B(*jk);
-            B.leftmultiply(*ij);
+            Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
             *ik -= B;
             ++ik; ++jk;
           }
@@ -82,7 +80,7 @@ namespace Dune {
       if (ij.index()!=i.index())
         DUNE_THROW(ISTLError,"diagonal entry missing");
       try {
-        (*ij).invert();   // compute inverse of diagonal block
+        Impl::asMatrix(*ij).invert();   // compute inverse of diagonal block
       }
       catch (Dune::FMatrixError & e) {
         DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
@@ -106,22 +104,36 @@ namespace Dune {
     rowiterator endi=A.end();
     for (rowiterator i=A.begin(); i!=endi; ++i)
     {
-      dblock rhs(d[i.index()]);
+      // We need to be careful here: Directly using
+      // auto rhs = Impl::asVector(d[ i.index() ]);
+      // is not OK in case this is a proxy. Hence
+      // we first have to copy the value. Notice that
+      // this is still not OK, if the vector type itself returns
+      // proxy references.
+      dblock rhsValue(d[i.index()]);
+      auto&& rhs = Impl::asVector(rhsValue);
       for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
-        (*j).mmv(v[j.index()],rhs);
-      v[i.index()] = rhs;           // Lii = I
+        Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
+      Impl::asVector(v[i.index()]) = rhs;           // Lii = I
     }
 
     // upper triangular solve
     rowiterator rendi=A.beforeBegin();
     for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
     {
-      vblock rhs(v[i.index()]);
+      // We need to be careful here: Directly using
+      // auto rhs = Impl::asVector(v[ i.index() ]);
+      // is not OK in case this is a proxy. Hence
+      // we first have to copy the value. Notice that
+      // this is still not OK, if the vector type itself returns
+      // proxy references.
+      vblock rhsValue(v[i.index()]);
+      auto&& rhs = Impl::asVector(rhsValue);
       coliterator j;
       for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
-        (*j).mmv(v[j.index()],rhs);
-      v[i.index()] = 0;
-      (*j).umv(rhs,v[i.index()]);           // diagonal stores inverse!
+        Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
+      auto&& vi = Impl::asVector(v[i.index()]);
+      Impl::asMatrix(*j).mv(rhs,vi);           // diagonal stores inverse!
     }
   }
 
@@ -129,19 +141,19 @@ namespace Dune {
 
   // recursive function template to access first entry of a matrix
   template<class M>
-  typename M::field_type& firstmatrixelement (M& A)
+  typename M::field_type& firstmatrixelement (M& A, typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
   {
     return firstmatrixelement(*(A.begin()->begin()));
   }
 
-  template<class K, int n, int m>
-  K& firstmatrixelement (FieldMatrix<K,n,m>& A)
+  template<class K>
+  K& firstmatrixelement (K& A, typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
   {
-    return A[0][0];
+    return A;
   }
 
-  template<class K>
-  K& firstmatrixelement (FieldMatrix<K,1,1>& A)
+  template<class K, int n, int m>
+  K& firstmatrixelement (FieldMatrix<K,n,m>& A)
   {
     return A[0][0];
   }
@@ -170,7 +182,6 @@ namespace Dune {
     createiterator ci=ILU.createbegin();
     for (crowiterator i=A.begin(); i!=endi; ++i)
     {
-      // std::cout << "in row " << i.index() << std::endl;
       map rowpattern; // maps column index to generation
 
       // initialize pattern with row of A
@@ -182,9 +193,6 @@ namespace Dune {
       {
         if ((*ik).second<n)
         {
-          //                            std::cout << "  eliminating " << i.index() << "," << (*ik).first
-          //                                              << " level " << (*ik).second << std::endl;
-
           coliterator endk = ILU[(*ik).first].end();                       // end of row k
           coliterator kj = ILU[(*ik).first].find((*ik).first);                       // diagonal in k
           for (++kj; kj!=endk; ++kj)                       // row k eliminates in row i
@@ -198,9 +206,6 @@ namespace Dune {
               mapiterator ij = rowpattern.find(kj.index());
               if (ij==rowpattern.end())
               {
-                //std::cout << "    new entry " << i.index() << "," << kj.index()
-                //                                                << " level " << (*ik).second+1 << std::endl;
-
                 rowpattern[kj.index()] = generation+1;
               }
             }
@@ -219,8 +224,6 @@ namespace Dune {
         firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
     }
 
-    //  printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
-
     // copy entries of A
     for (crowiterator i=A.begin(); i!=endi; ++i)
     {
@@ -391,33 +394,31 @@ namespace Dune {
       // lower triangular solve
       for( size_type i=0; i<iEnd; ++ i )
       {
-        dblock rhs( d[ i ] );
+        dblock rhsValue( d[ i ] );
+        auto&& rhs = Impl::asVector(rhsValue);
         const size_type rowI     = lower.rows_[ i ];
         const size_type rowINext = lower.rows_[ i+1 ];
 
         for( size_type col = rowI; col < rowINext; ++ col )
-        {
-          lower.values_[ col ].mmv( v[ lower.cols_[ col ] ], rhs );
-        }
+          Impl::asMatrix(lower.values_[ col ]).mmv( Impl::asVector(v[ lower.cols_[ col ] ] ), rhs );
 
-        v[ i ] = rhs;  // Lii = I
+        Impl::asVector(v[ i ]) = rhs;  // Lii = I
       }
 
       // upper triangular solve
       for( size_type i=0; i<iEnd; ++ i )
       {
-        vblock& vBlock = v[ lastRow - i ];
-        vblock rhs ( vBlock );
+        auto&& vBlock = Impl::asVector(v[ lastRow - i ]);
+        vblock rhsValue ( v[ lastRow - i ] );
+        auto&& rhs = Impl::asVector(rhsValue);
         const size_type rowI     = upper.rows_[ i ];
         const size_type rowINext = upper.rows_[ i+1 ];
 
         for( size_type col = rowI; col < rowINext; ++ col )
-        {
-          upper.values_[ col ].mmv( v[ upper.cols_[ col ] ], rhs );
-        }
+          Impl::asMatrix(upper.values_[ col ]).mmv( Impl::asVector(v[ upper.cols_[ col ] ]), rhs );
 
         // apply inverse and store result
-        inv[ i ].mv( rhs, vBlock);
+        Impl::asMatrix(inv[ i ]).mv(rhs, vBlock);
       }
     }
 

dune/istl/io.hh

@@ -15,13 +15,10 @@
 #include "istlexception.hh"
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/common/dynmatrix.hh>
-#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/hybridutilities.hh>
 #include <dune/common/unused.hh>
 
-#include <dune/istl/matrix.hh>
-#include <dune/istl/scaledidmatrix.hh>
-#include "bcrsmatrix.hh"
+#include <dune/istl/bcrsmatrix.hh>
 
 
 namespace Dune {
@@ -43,7 +40,7 @@ namespace Dune {
   //
 
   /**
-   * \brief Recursively print all the blocks
+   * \brief Recursively print a vector
    *
    * \code
    * #include <dune/istl/io.hh>
@@ -51,43 +48,30 @@ namespace Dune {
    */
   template<class V>
   void recursive_printvector (std::ostream& s, const V& v, std::string rowtext,
-                              int& counter, int columns, int width,
-                              int precision)
+                              int& counter, int columns, int width)
   {
-    for (typename V::ConstIterator i=v.begin(); i!=v.end(); ++i)
-      recursive_printvector(s,*i,rowtext,counter,columns,width,precision);
-  }
-
-  /**
-   * \brief Recursively print all the blocks -- specialization for FieldVector
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   */
-  template<class K, int n>
-  void recursive_printvector (std::ostream& s, const FieldVector<K,n>& v,
-                              std::string rowtext, int& counter, int columns,
-                              int width, int precision)
-  {
-    DUNE_UNUSED_PARAMETER(precision);
-    // we now can print n numbers
-    for (int i=0; i<n; i++)
-    {
-      if (counter%columns==0)
-      {
-        s << rowtext; // start a new row
-        s << " ";     // space in front of each entry
-        s.width(4);   // set width for counter
-        s << counter; // number of first entry in a line
-      }
-      s << " ";         // space in front of each entry
-      s.width(width);   // set width for each entry anew
-      s << v[i];        // yeah, the number !
-      counter++;        // increment the counter
-      if (counter%columns==0)
-        s << std::endl; // start a new line
-    }
+    Hybrid::ifElse(IsNumber<V>(),
+      [&](auto id) {
+        // Print one number
+        if (counter%columns==0)
+        {
+          s << rowtext; // start a new row
+          s << " ";     // space in front of each entry
+          s.width(4);   // set width for counter
+          s << counter; // number of first entry in a line
+        }
+        s << " ";         // space in front of each entry
+        s.width(width);   // set width for each entry anew
+        s << v;        // yeah, the number !
+        counter++;        // increment the counter
+        if (counter%columns==0)
+          s << std::endl; // start a new line
+      },
+      [&](auto id) {
+        // Recursively print a vector
+        for (const auto& entry : id(v))
+          recursive_printvector(s,id(entry),rowtext,counter,columns,width);
+      });
   }
 
 
@@ -119,7 +103,7 @@ namespace Dune {
       << std::endl;
 
     // print data from all blocks
-    recursive_printvector(s,v,rowtext,counter,columns,width,precision);
+    recursive_printvector(s,v,rowtext,counter,columns,width);
 
     // check if new line is required
     if (counter%columns!=0)
@@ -154,6 +138,31 @@ namespace Dune {
     }
   }
 
+  /**
+   * \brief Print one row of a matrix, specialization for number types
+   *
+   * \code
+   * #include <dune/istl/io.hh>
+   * \endcode
+   */
+  template<class K>
+  void print_row (std::ostream& s, const K& value,
+                  typename FieldMatrix<K,1,1>::size_type I,
+                  typename FieldMatrix<K,1,1>::size_type J,
+                  typename FieldMatrix<K,1,1>::size_type therow,
+                  int width, int precision,
+                  typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
+  {
+    DUNE_UNUSED_PARAMETER(I);
+    DUNE_UNUSED_PARAMETER(J);
+    DUNE_UNUSED_PARAMETER(therow);
+    DUNE_UNUSED_PARAMETER(precision);
+
+    s << " ";         // space in front of each entry
+    s.width(width);   // set width for each entry anew
+    s << value;
+  }
+
   /**
    * \brief Print one row of a matrix
    *
@@ -164,7 +173,8 @@ namespace Dune {
   template<class M>
   void print_row (std::ostream& s, const M& A, typename M::size_type I,
                   typename M::size_type J, typename M::size_type therow,
-                  int width, int precision)
+                  int width, int precision,
+                  typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
   {
     typename M::size_type i0=I;
     for (typename M::size_type i=0; i<A.N(); i++)
@@ -193,60 +203,6 @@ namespace Dune {
     }
   }
 
-  /**
-   * \brief Print one row of a matrix, specialization for FieldMatrix
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   */
-  template<class K, int n, int m>
-  void print_row (std::ostream& s, const FieldMatrix<K,n,m>& A,
-                  typename FieldMatrix<K,n,m>::size_type I,
-                  typename FieldMatrix<K,n,m>::size_type J,
-                  typename FieldMatrix<K,n,m>::size_type therow, int width,
-                  int precision)
-  {
-    DUNE_UNUSED_PARAMETER(J);
-    DUNE_UNUSED_PARAMETER(precision);
-
-    typedef typename FieldMatrix<K,n,m>::size_type size_type;
-
-    for (size_type i=0; i<n; i++)
-      if (I+i==therow)
-        for (int j=0; j<m; j++)
-        {
-          s << " ";         // space in front of each entry
-          s.width(width);   // set width for each entry anew
-          s << A[i][j];     // yeah, the number !
-        }
-  }
-
-  /**
-   * \brief Print one row of a matrix, specialization for FieldMatrix<K,1,1>
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   */
-  template<class K>
-  void print_row (std::ostream& s, const FieldMatrix<K,1,1>& A,
-                  typename FieldMatrix<K,1,1>::size_type I,
-                  typename FieldMatrix<K,1,1>::size_type J,
-                  typename FieldMatrix<K,1,1>::size_type therow,
-                  int width, int precision)
-  {
-    DUNE_UNUSED_PARAMETER(J);
-    DUNE_UNUSED_PARAMETER(precision);
-
-    if (I==therow)
-    {
-      s << " ";         // space in front of each entry
-      s.width(width);   // set width for each entry anew
-      s << static_cast<K>(A);         // yeah, the number !
-    }
-  }
-
   /**
    * \brief Print a generic block matrix
    *
@@ -419,80 +375,17 @@ namespace Dune {
    * #include <dune/istl/io.hh>
    * \endcode
    *
-   * This specialization for DiagonalMatrices ends the recursion
-   */
-  template <class FieldType, int dim>
-  void writeMatrixToMatlabHelper(const ScaledIdentityMatrix<FieldType,dim>& matrix, int rowOffset, int colOffset, std::ostream& s)
-  {
-    for (int i=0; i<dim; i++)
-    {
-      //+1 for Matlab numbering
-      s << rowOffset + i + 1 << " " << colOffset + i + 1 << " ";
-      MatlabPODWriter<FieldType>::write(matrix.scalar(), s)<< std::endl;
-    }
-  }
-
-  /**
-   * \brief Helper method for the writeMatrixToMatlab routine.
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   *
-   * This specialization for DiagonalMatrices ends the recursion
-   */
-  template <class FieldType, int dim>
-  void writeMatrixToMatlabHelper(const DiagonalMatrix<FieldType,dim>& matrix, int rowOffset, int colOffset, std::ostream& s)
-  {
-    for (int i=0; i<dim; i++)
-    {
-      //+1 for Matlab numbering
-      s << rowOffset + i + 1 << " " << colOffset + i + 1 << " ";
-      MatlabPODWriter<FieldType>::write(matrix.diagonal(i), s)<< std::endl;
-    }
-  }
-
-  /**
-   * \brief Helper method for the writeMatrixToMatlab routine.
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   *
-   * This specialization for FieldMatrices ends the recursion
-   */
-  template <class FieldType, int rows, int cols>
-  void writeMatrixToMatlabHelper
-    (   const FieldMatrix<FieldType,rows,cols>& matrix, int rowOffset,
-    int colOffset, std::ostream& s)
-  {
-    for (int i=0; i<rows; i++)
-      for (int j=0; j<cols; j++) {
-        //+1 for Matlab numbering
-        s << rowOffset + i + 1 << " " << colOffset + j + 1 << " ";
-        MatlabPODWriter<FieldType>::write(matrix[i][j], s)<< std::endl;
-      }
-  }
-
-  /**
-   * \brief Helper method for the writeMatrixToMatlab routine.
-   *
-   * \code
-   * #include <dune/istl/io.hh>
-   * \endcode
-   *
-   * This specialization for DynamicMatrices ends the recursion
+   * This specialization for numbers ends the recursion
    */
   template <class FieldType>
-  void writeMatrixToMatlabHelper(const DynamicMatrix<FieldType>& matrix, int rowOffset,
-                                 int colOffset, std::ostream& s)
+  void writeMatrixToMatlabHelper(const FieldType& value,
+                                 int rowOffset, int colOffset,
+                                 std::ostream& s,
+                                 typename std::enable_if_t<Dune::IsNumber<FieldType>::value>* sfinae = nullptr)
   {
-    for (int i=0; i<matrix.N(); i++)
-      for (int j=0; j<matrix.M(); j++) {
-        //+1 for Matlab numbering
-        s << rowOffset + i + 1 << " " << colOffset + j + 1 << " ";
-        MatlabPODWriter<FieldType>::write(matrix[i][j], s)<< std::endl;
-      }
+    //+1 for Matlab numbering
+    s << rowOffset + 1 << " " << colOffset + 1 << " ";
+    MatlabPODWriter<FieldType>::write(value, s)<< std::endl;
   }
 
   /**
@@ -505,7 +398,8 @@ namespace Dune {
   template <class MatrixType>
   void writeMatrixToMatlabHelper(const MatrixType& matrix,
                                  int externalRowOffset, int externalColOffset,
-                                 std::ostream& s)
+                                 std::ostream& s,
+                                 typename std::enable_if_t<!Dune::IsNumber<MatrixType>::value>* sfinae = nullptr)
   {
     // Precompute the accumulated sizes of the columns
     std::vector<typename MatrixType::size_type> colOffset(matrix.M());
@@ -521,11 +415,8 @@ namespace Dune {
     // Loop over all matrix rows
     for (typename MatrixType::size_type rowIdx=0; rowIdx<matrix.N(); rowIdx++)
     {
-
-      const typename MatrixType::row_type& row = matrix[rowIdx];
-
-      typename MatrixType::row_type::ConstIterator cIt   = row.begin();
-      typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+      auto cIt    = matrix[rowIdx].begin();
+      auto cEndIt = matrix[rowIdx].end();
 
       // Loop over all columns in this row
       for (; cIt!=cEndIt; ++cIt)
@@ -570,42 +461,18 @@ namespace Dune {
     outStream.precision(oldPrecision);
   }
 
-  // Recursively print all the blocks
+  // Recursively write vector entries to a stream
   template<class V>
   void writeVectorToMatlabHelper (const V& v, std::ostream& stream)
   {
-    for (const auto& entry : v)
-      writeVectorToMatlabHelper(entry, stream);
-  }
-
-  // Recursively print all the blocks -- specialization for FieldVector
-  template<class K, int n>
-  void writeVectorToMatlabHelper (const FieldVector<K,n>& v, std::ostream& s)
-  {
-    for (const auto& entry : v)
-    {
-      s << entry << std::endl;
-    }
-  }
-
-  // Recursively print all the blocks -- specialization for std::vector
-  template<class K>
-  void writeVectorToMatlabHelper (const std::vector<K>& v, std::ostream& s)
-  {
-    for (const auto& entry : v)
-    {
-      s << entry << std::endl;
-    }
-  }
-
-  // Recursively print all the blocks -- specialization for std::array
-  template<class K, std::size_t n>
-  void writeVectorToMatlabHelper (const std::array<K,n>& v, std::ostream& s)
-  {
-    for (const auto& entry : v)
-    {
-      s << entry << std::endl;
-    }
+    Hybrid::ifElse(IsNumber<V>(),
+      [&](auto id) {
+        stream << id(v) << std::endl;
+      },
+      [&](auto id) {
+        for (const auto& entry : id(v))
+          writeVectorToMatlabHelper(entry, stream);
+      });
   }
 
   /**

dune/istl/matrix.hh

@@ -12,6 +12,8 @@
 
 #include <dune/common/ftraits.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
 
 #include <dune/istl/bvector.hh>
 #include <dune/istl/istlexception.hh>
@@ -42,7 +44,7 @@ namespace MatrixImp
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the allocator type
     typedef A allocator_type;
@@ -557,7 +559,7 @@ namespace MatrixImp
   public:
 
     /** \brief Export the type representing the underlying field */
-    typedef typename T::field_type field_type;
+    using field_type = typename Imp::BlockTraits<T>::field_type;
 
     /** \brief Export the type representing the components */
     typedef T block_type;
@@ -583,10 +585,8 @@ namespace MatrixImp
     /** \brief Const iterator for the entries of each row */
     typedef typename row_type::const_iterator ConstColIterator;
 
-    enum {
-      //! The number of nesting levels the matrix contains.
-      blocklevel = T::blocklevel+1
-    };
+    //! The number of nesting levels the matrix contains.
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<T>::blockLevel()+1;
 
     /** \brief Create empty matrix */
     Matrix() : data_(0,0), cols_(0)
@@ -787,14 +787,15 @@ namespace MatrixImp
       if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-
       for (size_type i=0; i<data_.N(); i++) {
         y[i]=0;
         for (size_type j=0; j<cols_; j++)
-          (*this)[i][j].umv(x[j], y[i]);
-
+        {
+          auto&& xj = Impl::asVector(x[j]);
+          auto&& yi = Impl::asVector(y[i]);
+          Impl::asMatrix((*this)[i][j]).umv(xj, yi);
+        }
       }
-
     }
 
     //! y = A^T x
@@ -805,7 +806,6 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
 #endif
-
       for(size_type i=0; i<y.N(); ++i)
         y[i]=0;
       umtv(x,y);
@@ -819,14 +819,13 @@ namespace MatrixImp
       if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-
-      for (size_type i=0; i<data_.N(); i++) {
-
+      for (size_type i=0; i<data_.N(); i++)
         for (size_type j=0; j<cols_; j++)
-          (*this)[i][j].umv(x[j], y[i]);
-
-      }
-
+        {
+          auto&& xj = Impl::asVector(x[j]);
+          auto&& yi = Impl::asVector(y[i]);
+          Impl::asMatrix((*this)[i][j]).umv(xj, yi);
+        }
     }
 
     //! y -= A x
@@ -837,13 +836,13 @@ namespace MatrixImp
       if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmv(x[j.index()],y[i.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xj = Impl::asVector(x[j]);
+          auto&& yi = Impl::asVector(y[i]);
+          Impl::asMatrix((*this)[i][j]).mmv(xj, yi);
+        }
     }
 
     /** \brief \f$ y += \alpha A x \f$ */
@@ -854,14 +853,13 @@ namespace MatrixImp
       if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-
-      for (size_type i=0; i<data_.N(); i++) {
-
+      for (size_type i=0; i<data_.N(); i++)
         for (size_type j=0; j<cols_; j++)
-          (*this)[i][j].usmv(alpha, x[j], y[i]);
-
-      }
-
+        {
+          auto&& xj = Impl::asVector(x[j]);
+          auto&& yi = Impl::asVector(y[i]);
+          Impl::asMatrix((*this)[i][j]).usmv(alpha, xj, yi);
+        }
     }
 
     //! y += A^T x
@@ -872,13 +870,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umtv(x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).umtv(xi, yj);
+        }
     }
 
     //! y -= A^T x
@@ -889,13 +887,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmtv(x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).mmtv(xi, yj);
+        }
     }
 
     //! y += alpha A^T x
@@ -906,13 +904,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).usmtv(alpha,x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).usmtv(alpha, xi, yj);
+        }
     }
 
     //! y += A^H x
@@ -923,13 +921,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).umhv(x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).umhv(xi,yj);
+        }
     }
 
     //! y -= A^H x
@@ -940,13 +938,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).mmhv(x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).mmhv(xi,yj);
+        }
     }
 
     //! y += alpha A^H x
@@ -957,13 +955,13 @@ namespace MatrixImp
       if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
       if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
 #endif
-      ConstRowIterator endi=end();
-      for (ConstRowIterator i=begin(); i!=endi; ++i)
-      {
-        ConstColIterator endj = (*i).end();
-        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
-          (*j).usmhv(alpha,x[i.index()],y[j.index()]);
-      }
+      for (size_type i=0; i<data_.N(); i++)
+        for (size_type j=0; j<cols_; j++)
+        {
+          auto&& xi = Impl::asVector(x[i]);
+          auto&& yj = Impl::asVector(y[j]);
+          Impl::asMatrix((*this)[i][j]).usmhv(alpha,xi,yj);
+        }
     }
 
     //===== norms
@@ -977,10 +975,10 @@ namespace MatrixImp
     //! square of frobenius norm, need for block recursion
     typename FieldTraits<field_type>::real_type frobenius_norm2 () const
     {
-      double sum=0;
-      for (size_type i=0; i<N(); ++i)
-        for (size_type j=0; j<M(); ++j)
-          sum += data_[i][j].frobenius_norm2();
+      typename FieldTraits<field_type>::real_type sum=0;
+      for (size_type i=0; i<this->N(); i++)
+        for (size_type j=0; j<this->M(); j++)
+          sum += Impl::asMatrix(data_[i][j]).frobenius_norm2();
       return sum;
     }
 
@@ -995,9 +993,11 @@ namespace MatrixImp
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm();
+          sum += Impl::asMatrix(y).infinity_norm();
         norm = max(sum, norm);
+        isNaN += sum;
       }
+
       return norm;
     }
 
@@ -1012,7 +1012,7 @@ namespace MatrixImp
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm_real();
+          sum += Impl::asMatrix(y).infinity_norm_real();
         norm = max(sum, norm);
       }
       return norm;
@@ -1030,12 +1030,12 @@ namespace MatrixImp
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm();
+          sum += Impl::asMatrix(y).infinity_norm();
         norm = max(sum, norm);
         isNaN += sum;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+
+      return norm * (isNaN / isNaN);
     }
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
@@ -1050,12 +1050,12 @@ namespace MatrixImp
       for (auto const &x : *this) {
         real_type sum = 0;
         for (auto const &y : x)
-          sum += y.infinity_norm_real();
+          sum += Impl::asMatrix(y).infinity_norm_real();
         norm = max(sum, norm);
         isNaN += sum;
       }
-      isNaN /= isNaN;
-      return norm * isNaN;
+
+      return norm * (isNaN / isNaN);
     }
 
     //===== query

dune/istl/matrixmarket.hh

@@ -23,6 +23,7 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/unused.hh>
+#include <dune/common/hybridutilities.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
@@ -169,13 +170,13 @@ namespace Dune
     template<class M>
     struct mm_header_printer;
 
-    template<typename T, typename A, int i, int j>
-    struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
+    template<typename T, typename A>
+    struct mm_header_printer<BCRSMatrix<T,A> >
     {
       static void print(std::ostream& os)
       {
         os<<"%%MatrixMarket matrix coordinate ";
-        os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
+        os<<mm_numeric_type<typename Imp::BlockTraits<T>::field_type>::str()<<" general"<<std::endl;
       }
     };
 
@@ -185,7 +186,7 @@ namespace Dune
       static void print(std::ostream& os)
       {
         os<<"%%MatrixMarket matrix array ";
-        os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
+        os<<mm_numeric_type<typename Imp::BlockTraits<B>::field_type>::str()<<" general"<<std::endl;
       }
     };
 
@@ -220,6 +221,19 @@ namespace Dune
     template<class M>
     struct mm_block_structure_header;
 
+    template<typename T, typename A>
+    struct mm_block_structure_header<BlockVector<T,A> >
+    {
+      typedef BlockVector<T,A> M;
+      static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
+
+      static void print(std::ostream& os, const M&)
+      {
+        os<<"% ISTL_STRUCT blocked ";
+        os<<"1 1"<<std::endl;
+      }
+    };
+
     template<typename T, typename A, int i>
     struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
     {
@@ -232,6 +246,19 @@ namespace Dune
       }
     };
 
+    template<typename T, typename A>
+    struct mm_block_structure_header<BCRSMatrix<T,A> >
+    {
+      typedef BCRSMatrix<T,A> M;
+      static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
+
+      static void print(std::ostream& os, const M&)
+      {
+        os<<"% ISTL_STRUCT blocked ";
+        os<<"1 1"<<std::endl;
+      }
+    };
+
     template<typename T, typename A, int i, int j>
     struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
     {
@@ -647,23 +674,39 @@ namespace Dune
     struct MatrixValuesSetter
     {
       /**
-       * @brief Sets the matrixvalues.
+       * @brief Sets the matrix values.
        * @param row The row data as read from file.
        * @param matrix The matrix whose data we set.
        */
-      template<typename M>
+      template<typename T>
       void operator()(const std::vector<std::set<IndexData<D> > >& rows,
-                      M& matrix)
+                      BCRSMatrix<T>& matrix)
+      {
+        static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
+        for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
+        {
+          auto brow=iter.index();
+          for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
+            (*iter)[siter->index] = siter->number;
+        }
+      }
+
+      /**
+       * @brief Sets the matrix values.
+       * @param row The row data as read from file.
+       * @param matrix The matrix whose data we set.
+       */
+      template<typename T>
+      void operator()(const std::vector<std::set<IndexData<D> > >& rows,
+                      BCRSMatrix<FieldMatrix<T,brows,bcols> >& matrix)
       {
-        for(typename M::RowIterator iter=matrix.begin();
-            iter!= matrix.end(); ++iter)
+        for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
         {
-          for(typename M::size_type brow=iter.index()*brows,
+          for (auto brow=iter.index()*brows,
               browend=iter.index()*brows+brows;
               brow<browend; ++brow)
           {
-            typedef typename std::set<IndexData<D> >::const_iterator Siter;
-            for(Siter siter=rows[brow].begin(), send=rows[brow].end();
+            for (auto siter=rows[brow].begin(), send=rows[brow].end();
                 siter != send; ++siter)
               (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
           }
@@ -694,14 +737,41 @@ namespace Dune
       return std::conj(r);
     }
 
-    template<typename T, typename A, int brows, int bcols, typename D>
-    void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
+    template<typename M>
+    struct mm_multipliers
+    {};
+
+    template<typename B, typename A>
+    struct mm_multipliers<BCRSMatrix<B,A> >
+    {
+      enum {
+        rows = 1,
+        cols = 1
+      };
+    };
+
+    template<typename B, int i, int j, typename A>
+    struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
+    {
+      enum {
+        rows = i,
+        cols = j
+      };
+    };
+
+    template<typename T, typename A, typename D>
+    void readSparseEntries(Dune::BCRSMatrix<T,A>& matrix,
                            std::istream& file, std::size_t entries,
                            const MMHeader& mmHeader, const D&)
     {
-      typedef Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A> Matrix;
+      typedef Dune::BCRSMatrix<T,A> Matrix;
+
+      // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
+      constexpr int brows = mm_multipliers<Matrix>::rows;
+      constexpr int bcols = mm_multipliers<Matrix>::cols;
+
       // First path
-      // store entries together with column index in a speparate
+      // store entries together with column index in a separate
       // data structure
       std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
 
@@ -810,6 +880,15 @@ namespace Dune
     istr >> cols;
   }
 
+  template<typename T, typename A>
+  void mm_read_vector_entries(Dune::BlockVector<T,A>& vector,
+                              std::size_t size,
+                              std::istream& istr)
+  {
+    for (int i=0; size>0; ++i, --size)
+      istr>>vector[i];
+  }
+
   template<typename T, typename A, int entries>
   void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
                               std::size_t size,
@@ -829,8 +908,8 @@ namespace Dune
    * @param istr The input stream to read the data from.
    * @warning Not all formats are supported!
    */
-  template<typename T, typename A, int entries>
-  void readMatrixMarket(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
+  template<typename T, typename A>
+  void readMatrixMarket(Dune::BlockVector<T,A>& vector,
                         std::istream& istr)
   {
     using namespace MatrixMarketImpl;
@@ -844,29 +923,37 @@ namespace Dune
     if(header.type!=array_type)
       DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
 
-    std::size_t size=rows/entries;
-    if(size*entries!=rows)
-      DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
 
-    vector.resize(size);
+    Dune::Hybrid::ifElse(Dune::IsNumber<T>(),
+      [&](auto id){
+        vector.resize(rows);
+      },
+      [&](auto id){  // Assume that T is a FieldVector
+        T dummy;
+        auto blocksize = id(dummy).size();
+        std::size_t size=rows/blocksize;
+        if(size*blocksize!=rows)
+          DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
 
-    istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+        vector.resize(size);
+      });
 
+    istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
     mm_read_vector_entries(vector, rows, istr);
   }
 
-
   /**
    * @brief Reads a sparse matrix from a matrix market file.
    * @param matrix The matrix to store the data in.
    * @param istr The input stream to read the data from.
    * @warning Not all formats are supported!
    */
-  template<typename T, typename A, int brows, int bcols>
-  void readMatrixMarket(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
+  template<typename T, typename A>
+  void readMatrixMarket(Dune::BCRSMatrix<T,A>& matrix,
                         std::istream& istr)
   {
     using namespace MatrixMarketImpl;
+    using Matrix = Dune::BCRSMatrix<T,A>;
 
     MMHeader header;
     if(!readMatrixMarketBanner(istr, header)) {
@@ -896,46 +983,42 @@ namespace Dune
 
     std::size_t nnz, blockrows, blockcols;
 
+    // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
+    constexpr int brows = mm_multipliers<Matrix>::rows;
+    constexpr int bcols = mm_multipliers<Matrix>::cols;
+
     std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
 
     istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
 
 
     matrix.setSize(blockrows, blockcols);
-    matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
+    matrix.setBuildMode(Dune::BCRSMatrix<T,A>::row_wise);
 
     if(header.type==array_type)
       DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
 
-    readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
+    readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
   }
 
-  template<typename M>
-  struct mm_multipliers
-  {};
-
-  template<typename B, int i, int j, typename A>
-  struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
-  {
-    enum {
-      rows = i,
-      cols = j
-    };
-  };
-
-  template<typename B, int i, int j>
-  void mm_print_entry(const FieldMatrix<B,i,j>& entry,
-                      typename FieldMatrix<B,i,j>::size_type rowidx,
-                      typename FieldMatrix<B,i,j>::size_type colidx,
+  // Print a scalar entry
+  template<typename B>
+  void mm_print_entry(const B& entry,
+                      std::size_t rowidx,
+                      std::size_t colidx,
                       std::ostream& ostr)
   {
-    typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
-    typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
-    for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
-      int coli=colidx;
-      for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
-        ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
-    }
+    Hybrid::ifElse(IsNumber<B>(),
+      [&](auto id) {
+        ostr << rowidx << " " << colidx << " " << entry << std::endl;
+      },
+      [&](auto id) {
+        for (auto row=id(entry).begin(); row != id(entry).end(); ++row, ++rowidx) {
+          int coli=colidx;
+          for (auto col = row->begin(); col != row->end(); ++col, ++coli)
+            ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
+        }
+      });
   }
 
   // Write a vector entry
@@ -963,6 +1046,12 @@ namespace Dune
                             std::integral_constant<int,isnumeric>());
   }
 
+  template<typename T, typename A>
+  std::size_t countEntries(const BlockVector<T,A>& vector)
+  {
+    return vector.size();
+  }
+
   template<typename T, typename A, int i>
   std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
   {
@@ -987,8 +1076,8 @@ namespace Dune
                          std::ostream& ostr,
                          const std::integral_constant<int,1>&)
   {
-    ostr<<matrix.N()*mm_multipliers<M>::rows<<" "
-        <<matrix.M()*mm_multipliers<M>::cols<<" "
+    ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
+        <<matrix.M()*MatrixMarketImpl::mm_multipliers<M>::cols<<" "
         <<countNonZeros(matrix)<<std::endl;
 
     typedef typename M::const_iterator riterator;
@@ -996,8 +1085,8 @@ namespace Dune
     for(riterator row=matrix.begin(); row != matrix.end(); ++row)
       for(citerator col = row->begin(); col != row->end(); ++col)
         // Matrix Market indexing start with 1!
-        mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
-                       col.index()*mm_multipliers<M>::cols+1, ostr);
+        mm_print_entry(*col, row.index()*MatrixMarketImpl::mm_multipliers<M>::rows+1,
+                       col.index()*MatrixMarketImpl::mm_multipliers<M>::cols+1, ostr);
   }
 
 

dune/istl/matrixredistribute.hh

@@ -446,10 +446,12 @@ namespace Dune
           {
             *c=0;
             if(c.index()==i->local()) {
-              typedef typename M::block_type::RowIterator RIter;
-              for(RIter r=c->begin(), rend=c->end();
-                  r != rend; ++r)
-                (*r)[r.index()]=1;
+              auto setDiagonal = [](auto&& scalarOrMatrix, const auto& value) {
+                auto&& matrix = Dune::Impl::asMatrix(scalarOrMatrix);
+                for (auto rowIt = matrix.begin(); rowIt != matrix.end(); ++rowIt)
+                  (*rowIt)[rowIt.index()] = value;
+              };
+              setDiagonal(*c, 1);
             }
           }
         }