diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index fc364b370..68058f6a9 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -52,6 +52,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/linalg/FlexibleSolver5.cpp
   opm/simulators/linalg/FlexibleSolver6.cpp
   opm/simulators/linalg/MILU.cpp
+  opm/simulators/linalg/ParallelIstlInformation.cpp
   opm/simulators/linalg/PropertyTree.cpp
   opm/simulators/linalg/setupPropertyTree.cpp
   opm/simulators/utils/PartiallySupportedFlowKeywords.cpp
diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp
index e0c6643cd..673143127 100644
--- a/opm/simulators/flow/BlackoilModelEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelEbos.hpp
@@ -39,7 +39,6 @@
 
 #include <opm/grid/UnstructuredGrid.h>
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
-#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>
diff --git a/opm/simulators/flow/countGlobalCells.cpp b/opm/simulators/flow/countGlobalCells.cpp
index d72c550c3..240ddc85a 100644
--- a/opm/simulators/flow/countGlobalCells.cpp
+++ b/opm/simulators/flow/countGlobalCells.cpp
@@ -21,9 +21,12 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <config.h>
+
 #include "countGlobalCells.hpp"
 
 #include <cassert>
+#include <numeric>
 
 namespace Opm {
 namespace detail {
diff --git a/opm/simulators/flow/countGlobalCells.hpp b/opm/simulators/flow/countGlobalCells.hpp
index de5fe2459..050a668a2 100644
--- a/opm/simulators/flow/countGlobalCells.hpp
+++ b/opm/simulators/flow/countGlobalCells.hpp
@@ -24,7 +24,6 @@
 #ifndef OPM_COUNTGLOBALCELLS_HEADER_INCLUDED
 #define OPM_COUNTGLOBALCELLS_HEADER_INCLUDED
 
-#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 
 #include <dune/grid/common/gridview.hh>
diff --git a/opm/simulators/linalg/ParallelIstlInformation.cpp b/opm/simulators/linalg/ParallelIstlInformation.cpp
new file mode 100644
index 000000000..7cd233d26
--- /dev/null
+++ b/opm/simulators/linalg/ParallelIstlInformation.cpp
@@ -0,0 +1,355 @@
+/*
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2014, 2015 Statoil ASA
+  Copyright 2015 NTNU
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_MPI && HAVE_DUNE_ISTL
+
+#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
+
+#include <dune/common/enumset.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <cstddef>
+#include <exception>
+#include <mpi.h>
+#include <numeric>
+
+namespace
+{
+
+template<class T>
+class IndexSetInserter
+{
+public:
+    using ParallelIndexSet = T;
+    using LocalIndex = typename ParallelIndexSet::LocalIndex;
+    using GlobalIndex = typename ParallelIndexSet::GlobalIndex;
+
+    IndexSetInserter(ParallelIndexSet& indexSet, const GlobalIndex& component_size,
+                     std::size_t local_component_size, std::size_t num_components)
+        : indexSet_(&indexSet), component_size_(component_size),
+          local_component_size_(local_component_size),
+          num_components_(num_components)
+    {}
+
+    void operator()(const typename ParallelIndexSet::IndexPair& pair)
+    {
+        for(std::size_t i = 0; i < num_components_; i++)
+            indexSet_->add(i * component_size_ + pair.global(),
+                           LocalIndex(i * local_component_size_  + pair.local(),
+                                      pair.local().attribute()));
+    }
+private:
+    ParallelIndexSet* indexSet_;
+    /// \brief The global number of unknowns per component/equation.
+    GlobalIndex component_size_;
+    /// \brief The local number of unknowns per component/equation.
+    std::size_t local_component_size_;
+    /// \brief The number of components/equations.
+    std::size_t num_components_;
+};
+
+/** \brief gather/scatter callback for communcation */
+template<typename T>
+struct CopyGatherScatter
+{
+    using V = typename Dune::CommPolicy<T>::IndexedType;
+
+    static V gather(const T& a, std::size_t i)
+    {
+      return a[i];
+    }
+
+    static void scatter(T& a, V v, std::size_t i)
+    {
+      a[i] = v;
+    }
+};
+
+template<int I=0, typename... BinaryOperators, typename... ReturnValues>
+typename std::enable_if<I == sizeof...(BinaryOperators), void>::type
+computeGlobalReduction(const std::tuple<ReturnValues...>&,
+                       std::tuple<BinaryOperators...>&,
+                       std::tuple<ReturnValues...>&)
+{}
+
+template<int I=0, typename... BinaryOperators, typename... ReturnValues>
+typename std::enable_if<I !=sizeof...(BinaryOperators), void>::type
+computeGlobalReduction(const std::tuple<ReturnValues...>& receivedValues,
+                       std::tuple<BinaryOperators...>& operators,
+                       std::tuple<ReturnValues...>& values)
+{
+    auto& val = std::get<I>(values);
+    val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
+    computeGlobalReduction<I+1>(receivedValues, operators, values);
+}
+
+template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
+typename std::enable_if<I==sizeof...(Containers), void>::type
+computeLocalReduction(const std::tuple<Containers...>&,
+                      std::tuple<BinaryOperators...>&,
+                      std::tuple<ReturnValues...>&,
+                      const std::vector<double>&)
+{}
+
+template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
+typename std::enable_if<I!=sizeof...(Containers), void>::type
+computeLocalReduction(const std::tuple<Containers...>& containers,
+                      std::tuple<BinaryOperators...>& operators,
+                      std::tuple<ReturnValues...>& values,
+                      const std::vector<double>& ownerMask)
+{
+    const auto& container = std::get<I>(containers);
+    if (container.size())
+    {
+        auto& reduceOperator = std::get<I>(operators);
+        // Eigen:Block does not support STL iterators!!!!
+        // Therefore we need to rely on the harder random-access
+        // property of the containers. But this should be save, too.
+        // Just commenting out code in the hope that Eigen might improve
+        // in this regard in the future.
+        //auto newVal = container.begin();
+        auto mask   = ownerMask.begin();
+        auto& value = std::get<I>(values);
+        value = reduceOperator.getInitialValue();
+
+        for (auto endVal = ownerMask.end(); mask != endVal; /*++newVal,*/ ++mask )
+        {
+            value = reduceOperator(value, container[mask-ownerMask.begin()], *mask);
+        }
+    }
+    computeLocalReduction<I+1>(containers, operators, values, ownerMask);
+}
+
+}
+
+namespace Opm
+{
+namespace
+{
+
+    template<class T>
+    struct is_tuple
+        : std::integral_constant<bool, false>
+    {};
+    template<typename... T>
+    struct is_tuple<std::tuple<T...> >
+        : std::integral_constant<bool, true>
+    {};
+}
+
+ParallelISTLInformation::ParallelISTLInformation()
+    : indexSet_(new ParallelIndexSet),
+      remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
+      communicator_(MPI_COMM_WORLD)
+{}
+
+
+ParallelISTLInformation::ParallelISTLInformation(MPI_Comm communicator)
+    : indexSet_(new ParallelIndexSet),
+      remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
+      communicator_(communicator)
+{}
+
+
+ParallelISTLInformation::
+ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
+                        const std::shared_ptr<RemoteIndices>& remoteIndices,
+                        MPI_Comm communicator)
+    : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
+{}
+
+ParallelISTLInformation::ParallelISTLInformation(const ParallelISTLInformation& other)
+    : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
+      communicator_(other.communicator_)
+{}
+
+void ParallelISTLInformation::copyValuesTo(ParallelIndexSet& indexSet,
+                                           RemoteIndices& remoteIndices,
+                                           std::size_t local_component_size,
+                                           std::size_t num_components) const
+{
+    ParallelIndexSet::GlobalIndex global_component_size  = local_component_size;
+    if ( num_components > 1 )
+    {
+        ParallelIndexSet::GlobalIndex max_gi = 0;
+        // component the max global index
+        for( auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
+        {
+            max_gi = std::max(max_gi, i->global());
+        }
+        global_component_size = max_gi+1;
+        global_component_size = communicator_.max(global_component_size);
+    }
+    indexSet.beginResize();
+    IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
+                                                local_component_size, num_components);
+    std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
+    indexSet.endResize();
+    remoteIndices.rebuild<false>();
+}
+
+template<class T>
+void ParallelISTLInformation::copyOwnerToAll(const T& source, T& dest) const
+{
+    using AS = Dune::OwnerOverlapCopyAttributeSet;
+    using CopySet = Dune::EnumItem<AS, AS::copy>;
+    using OwnerSet = Dune::EnumItem<AS, AS::owner>;
+    using OverlapSet = Dune::EnumItem<AS, AS::overlap>;
+    using OwnerOverlapSet = Dune::Combine<OwnerSet, OverlapSet, AS::AttributeSet>;
+    using AllSet = Dune::Combine<OwnerOverlapSet, CopySet, AS::AttributeSet>;
+    OwnerSet sourceFlags;
+    AllSet destFlags;
+    Dune::Interface interface(communicator_);
+    if( !remoteIndices_->isSynced() )
+    {
+        remoteIndices_->rebuild<false>();
+    }
+    interface.build(*remoteIndices_,sourceFlags,destFlags);
+    Dune::BufferedCommunicator communicator;
+    communicator.template build<T>(interface);
+    communicator.template forward<CopyGatherScatter<T>>(source,dest);
+    communicator.free();
+}
+
+template<class T>
+const std::vector<double>&
+ParallelISTLInformation::updateOwnerMask(const T& container) const
+{
+    if (!indexSet_)
+    {
+        OPM_THROW(std::runtime_error, "Trying to update owner mask without parallel information!");
+    }
+    if (static_cast<std::size_t>(container.size()) != ownerMask_.size())
+    {
+        ownerMask_.resize(container.size(), 1.);
+        for (const auto& i : *indexSet_)
+        {
+            if (i.local().attribute() != Dune::OwnerOverlapCopyAttributeSet::owner)
+            {
+                ownerMask_[i.local().local()] = 0.;
+            }
+        }
+    }
+    return ownerMask_;
+}
+
+template<typename Container, typename BinaryOperator, typename T>
+void ParallelISTLInformation::computeReduction(const Container& container,
+                                               BinaryOperator binaryOperator,
+                                               T& value) const
+{
+    if constexpr (is_tuple<Container>())
+        computeTupleReduction(container, binaryOperator, value);
+    else
+    {
+        std::tuple<const Container&> containers = std::tuple<const Container&>(container);
+        auto values = std::make_tuple(value);
+        auto operators = std::make_tuple(binaryOperator);
+        computeTupleReduction(containers, operators, values);
+        value = std::get<0>(values);
+    }
+}
+
+template<typename... Containers, typename... BinaryOperators, typename... ReturnValues>
+void ParallelISTLInformation::computeTupleReduction(const std::tuple<Containers...>& containers,
+                                                    std::tuple<BinaryOperators...>& operators,
+                                                    std::tuple<ReturnValues...>& values) const
+{
+    static_assert(std::tuple_size<std::tuple<Containers...> >::value ==
+                  std::tuple_size<std::tuple<BinaryOperators...> >::value,
+                  "We need the same number of containers and binary operators");
+    static_assert(std::tuple_size<std::tuple<Containers...> >::value ==
+                  std::tuple_size<std::tuple<ReturnValues...> >::value,
+                  "We need the same number of containers and return values");
+    if (std::tuple_size<std::tuple<Containers...> >::value == 0)
+    {
+        return;
+    }
+
+    // Copy the initial values.
+    std::tuple<ReturnValues...> init = values;
+    updateOwnerMask(std::get<0>(containers));
+    computeLocalReduction(containers, operators, values, ownerMask_);
+    std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
+    communicator_.allgather(&values, 1, &(receivedValues[0]));
+    values = init;
+    for (auto& rval : receivedValues)
+    {
+        computeGlobalReduction(rval, operators, values);
+    }
+}
+
+template<class T1>
+auto
+accumulateMaskedValues(const T1& container, const std::vector<double>* maskContainer)
+    -> decltype(container[0]*(*maskContainer)[0])
+{
+    decltype(container[0]*(*maskContainer)[0]) initial = 0;
+
+    if (maskContainer)
+    {
+        return std::inner_product(container.begin(), container.end(),
+                                  maskContainer->begin(), initial);
+    }
+    else
+    {
+        return std::accumulate(container.begin(), container.end(), initial);
+    }
+}
+
+template<class T> using C1 = std::vector<T>;
+template<class T> using Ops1 = Reduction::MaskIDOperator<std::plus<T>>;
+
+template<class T>
+using C2 = std::tuple<std::vector<T>,
+                      std::vector<T>,
+                      std::vector<T>,
+                      std::vector<T>,
+                      std::vector<T>>;
+template<class T>
+using Ops2 = std::tuple<decltype(Reduction::makeGlobalSumFunctor<T>()),
+                        decltype(Reduction::makeGlobalMaxFunctor<T>()),
+                        decltype(Reduction::makeGlobalMinFunctor<T>()),
+                        decltype(Reduction::makeInnerProductFunctor<T>()),
+                        decltype(Reduction::makeLInfinityNormFunctor<T>())>;
+template<class T>
+using Vals2 = std::tuple<T,T,T,T,T>;
+
+#define INSTANCE1(T) \
+    template void ParallelISTLInformation::computeReduction<C1<T>,Ops1<T>,T>(const C1<T>&,Ops1<T>,T&) const;
+
+#define INSTANCE2(T) \
+    template void ParallelISTLInformation::computeReduction<C2<T>,Ops2<T>,Vals2<T>>(const C2<T>&,Ops2<T>,Vals2<T>&) const;
+
+#define INSTANCE(T) \
+    INSTANCE1(T) \
+    INSTANCE2(T)
+
+INSTANCE(int)
+INSTANCE(float)
+INSTANCE(std::size_t)
+
+#endif
+
+}
diff --git a/opm/simulators/linalg/ParallelIstlInformation.hpp b/opm/simulators/linalg/ParallelIstlInformation.hpp
index 624013e1b..29e623c25 100644
--- a/opm/simulators/linalg/ParallelIstlInformation.hpp
+++ b/opm/simulators/linalg/ParallelIstlInformation.hpp
@@ -18,48 +18,25 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
-#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
+#ifndef OPM_PARALLELISTLINFORMATION_HEADER_INCLUDED
+#define OPM_PARALLELISTLINFORMATION_HEADER_INCLUDED
 
-
-#include <opm/grid/UnstructuredGrid.h>
-#include <opm/common/ErrorMacros.hpp>
-#include <any>
-#include <exception>
-
-#include <algorithm>
-#include <functional>
-#include <limits>
-#include <numeric>
-#include <type_traits>
 #include <vector>
 
 #if HAVE_MPI && HAVE_DUNE_ISTL
 
-#include <opm/common/utility/platform_dependent/disable_warnings.h>
-#include <mpi.h>
+#include <algorithm>
+#include <limits>
+#include <memory>
+#include <tuple>
+#include <type_traits>
+
 #include <dune/istl/owneroverlapcopy.hh>
-#include <dune/common/parallel/interface.hh>
-#include <dune/common/parallel/communicator.hh>
-#include <dune/common/enumset.hh>
-#include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 namespace Opm
 {
-namespace
-{
-
-    template<class T>
-    struct is_tuple
-        : std::integral_constant<bool, false>
-    {};
-    template<typename... T>
-    struct is_tuple<std::tuple<T...> >
-        : std::integral_constant<bool, true>
-    {};
-}
 
 /// \brief Class that encapsulates the parallelization information needed by the
 /// ISTL solvers.
@@ -67,49 +44,42 @@ class ParallelISTLInformation
 {
 public:
     /// \brief The type of the parallel index set used.
-    typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
+    using ParallelIndexSet = Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet;
     /// \brief The type of the remote indices information used.
-    typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
+    using RemoteIndices = Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices;
 
     /// \brief Constructs an empty parallel information object using MPI_COMM_WORLD
-    ParallelISTLInformation()
-        : indexSet_(new ParallelIndexSet),
-          remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
-          communicator_(MPI_COMM_WORLD)
-    {}
+    ParallelISTLInformation();
+
     /// \brief Constructs an empty parallel information object using a communicator.
     /// \param communicator The communicator to use.
-    ParallelISTLInformation(MPI_Comm communicator)
-        : indexSet_(new ParallelIndexSet),
-          remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
-          communicator_(communicator)
-    {}
+    ParallelISTLInformation(MPI_Comm communicator);
+
     /// \brief Constructs a parallel information object from the specified information.
     /// \param indexSet The parallel index set to use.
     /// \param remoteIndices The remote indices information to use.
     /// \param communicator The communicator to use.
     ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
                             const std::shared_ptr<RemoteIndices>& remoteIndices,
-                            MPI_Comm communicator)
-        : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
-    {}
+                            MPI_Comm communicator);
+
     /// \brief Copy constructor.
     ///
     /// The information will be shared by the the two objects.
-    ParallelISTLInformation(const ParallelISTLInformation& other)
-    : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
-      communicator_(other.communicator_)
-    {}
+    ParallelISTLInformation(const ParallelISTLInformation& other);
+
     /// \brief Get a pointer to the underlying index set.
     std::shared_ptr<ParallelIndexSet> indexSet() const
     {
         return indexSet_;
     }
+
     /// \brief Get a pointer to the remote indices information.
     std::shared_ptr<RemoteIndices> remoteIndices() const
     {
         return remoteIndices_;
     }
+
     /// \brief Get the Collective MPI communicator that we use.
     Parallel::Communication communicator() const
     {
@@ -119,69 +89,17 @@ public:
     /// \param[out] indexSet The object to store the index set in.
     /// \param[out] remoteIndices The object to store the remote indices information in.
     void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices,
-                      std::size_t local_component_size = 0, std::size_t num_components = 1) const
-    {
-        ParallelIndexSet::GlobalIndex global_component_size  = local_component_size;
-        if ( num_components > 1 )
-        {
-            ParallelIndexSet::GlobalIndex max_gi = 0;
-            // component the max global index
-            for( auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
-            {
-                max_gi = std::max(max_gi, i->global());
-            }
-            global_component_size = max_gi+1;
-            global_component_size = communicator_.max(global_component_size);
-        }
-        indexSet.beginResize();
-        IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
-                                                    local_component_size, num_components);
-        std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
-        indexSet.endResize();
-        remoteIndices.rebuild<false>();
-    }
+                      std::size_t local_component_size = 0,
+                      std::size_t num_components = 1) const;
+
     /// \brief Communcate the dofs owned by us to the other process.
     ///
     /// Afterwards all associated dofs will contain the same data.
     template<class T>
-    void copyOwnerToAll (const T& source, T& dest) const
-    {
-        typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
-        typedef Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner> OwnerSet;
-        typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
-      OwnerSet sourceFlags;
-      AllSet destFlags;
-      Dune::Interface interface(communicator_);
-      if( !remoteIndices_->isSynced() )
-      {
-          remoteIndices_->rebuild<false>();
-      }
-      interface.build(*remoteIndices_,sourceFlags,destFlags);
-      Dune::BufferedCommunicator communicator;
-      communicator.template build<T>(interface);
-      communicator.template forward<CopyGatherScatter<T> >(source,dest);
-      communicator.free();
-    }
+    void copyOwnerToAll (const T& source, T& dest) const;
+
     template<class T>
-    const std::vector<double>& updateOwnerMask(const T& container) const
-    {
-        if( ! indexSet_ )
-        {
-            OPM_THROW(std::runtime_error, "Trying to update owner mask without parallel information!");
-        }
-        if( static_cast<std::size_t>(container.size())!= ownerMask_.size() )
-        {
-            ownerMask_.resize(container.size(), 1.);
-            for( auto i=indexSet_->begin(), end=indexSet_->end(); i!=end; ++i )
-            {
-                if (i->local().attribute()!=Dune::OwnerOverlapCopyAttributeSet::owner)
-                {
-                    ownerMask_[i->local().local()] = 0.;
-                }
-            }
-        }
-        return ownerMask_;
-    }
+    const std::vector<double>& updateOwnerMask(const T& container) const;
 
     /// \brief Get the owner Mask.
     ///
@@ -214,166 +132,14 @@ public:
     /// \param value The initial value or a tuple of them.
     template<typename Container, typename BinaryOperator, typename T>
     void computeReduction(const Container& container, BinaryOperator binaryOperator,
-                          T& value) const
-    {
-        computeReduction(container, binaryOperator, value, is_tuple<Container>());
-    }
+                          T& value) const;
+
 private:
-    /// \brief compute the reductions for tuples.
-    ///
-    /// This is a helper function to prepare for calling computeTupleReduction.
-    template<typename Container, typename BinaryOperator, typename T>
-    void computeReduction(const Container& container, BinaryOperator binaryOperator,
-                          T& value, std::integral_constant<bool,true>) const
-    {
-        computeTupleReduction(container, binaryOperator, value);
-    }
-    /// \brief compute the reductions for non-tuples.
-    ///
-    /// This is a helper function to prepare for calling computeTupleReduction.
-    template<typename Container, typename BinaryOperator, typename T>
-    void computeReduction(const Container& container, BinaryOperator binaryOperator,
-                          T& value, std::integral_constant<bool,false>) const
-    {
-        std::tuple<const Container&> containers=std::tuple<const Container&>(container);
-        auto values=std::make_tuple(value);
-        auto operators=std::make_tuple(binaryOperator);
-        computeTupleReduction(containers, operators, values);
-        value=std::get<0>(values);
-    }
-    /// \brief Compute the reductions for tuples.
     template<typename... Containers, typename... BinaryOperators, typename... ReturnValues>
     void computeTupleReduction(const std::tuple<Containers...>& containers,
                                std::tuple<BinaryOperators...>& operators,
-                               std::tuple<ReturnValues...>& values) const
-    {
-        static_assert(std::tuple_size<std::tuple<Containers...> >::value==
-                      std::tuple_size<std::tuple<BinaryOperators...> >::value,
-                      "We need the same number of containers and binary operators");
-        static_assert(std::tuple_size<std::tuple<Containers...> >::value==
-                      std::tuple_size<std::tuple<ReturnValues...> >::value,
-                      "We need the same number of containers and return values");
-        if( std::tuple_size<std::tuple<Containers...> >::value==0 )
-        {
-            return;
-        }
-        // Copy the initial values.
-        std::tuple<ReturnValues...> init=values;
-        updateOwnerMask(std::get<0>(containers));
-        computeLocalReduction(containers, operators, values);
-        std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
-        communicator_.allgather(&values, 1, &(receivedValues[0]));
-        values=init;
-        for( auto rvals=receivedValues.begin(), endvals=receivedValues.end(); rvals!=endvals;
-             ++rvals )
-        {
-            computeGlobalReduction(*rvals, operators, values);
-        }
-    }
-    /// \brief TMP for computing the the global reduction after receiving the local ones.
-    ///
-    /// End of recursion.
-    template<int I=0, typename... BinaryOperators, typename... ReturnValues>
-    typename std::enable_if<I == sizeof...(BinaryOperators), void>::type
-    computeGlobalReduction(const std::tuple<ReturnValues...>&,
-                                std::tuple<BinaryOperators...>&,
-                                std::tuple<ReturnValues...>&) const
-    {}
-    /// \brief TMP for computing the the global reduction after receiving the local ones.
-    template<int I=0, typename... BinaryOperators, typename... ReturnValues>
-    typename std::enable_if<I !=sizeof...(BinaryOperators), void>::type
-    computeGlobalReduction(const std::tuple<ReturnValues...>& receivedValues,
-                           std::tuple<BinaryOperators...>& operators,
-                           std::tuple<ReturnValues...>& values) const
-    {
-        auto& val=std::get<I>(values);
-        val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
-        computeGlobalReduction<I+1>(receivedValues, operators, values);
-    }
-    /// \brief TMP for computing the the local reduction on the DOF that the process owns.
-    ///
-    /// End of recursion.
-    template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
-    typename std::enable_if<I==sizeof...(Containers), void>::type
-    computeLocalReduction(const std::tuple<Containers...>&,
-                          std::tuple<BinaryOperators...>&,
-                          std::tuple<ReturnValues...>&) const
-    {}
-    /// \brief TMP for computing the the local reduction on the DOF that the process owns.
-    template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
-    typename std::enable_if<I!=sizeof...(Containers), void>::type
-    computeLocalReduction(const std::tuple<Containers...>& containers,
-                          std::tuple<BinaryOperators...>& operators,
-                          std::tuple<ReturnValues...>& values) const
-    {
-        const auto& container = std::get<I>(containers);
-        if( container.size() )
-        {
-            auto& reduceOperator  = std::get<I>(operators);
-            // Eigen:Block does not support STL iterators!!!!
-            // Therefore we need to rely on the harder random-access
-            // property of the containers. But this should be save, too.
-            // Just commenting out code in the hope that Eigen might improve
-            // in this regard in the future.
-            //auto newVal = container.begin();
-            auto mask   = ownerMask_.begin();
-            auto& value = std::get<I>(values);
-            value =  reduceOperator.getInitialValue();
+                               std::tuple<ReturnValues...>& values) const;
 
-            for( auto endVal=ownerMask_.end(); mask!=endVal;
-                 /*++newVal,*/ ++mask )
-            {
-                value = reduceOperator(value, container[mask-ownerMask_.begin()], *mask);
-            }
-        }
-        computeLocalReduction<I+1>(containers, operators, values);
-    }
-    /** \brief gather/scatter callback for communcation */
-    template<typename T>
-    struct CopyGatherScatter
-    {
-        typedef typename Dune::CommPolicy<T>::IndexedType V;
-
-      static V gather(const T& a, std::size_t i)
-      {
-        return a[i];
-      }
-
-      static void scatter(T& a, V v, std::size_t i)
-      {
-        a[i] = v;
-      }
-    };
-    template<class T>
-    class IndexSetInserter
-    {
-    public:
-        typedef T ParallelIndexSet;
-        typedef typename ParallelIndexSet::LocalIndex LocalIndex;
-        typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
-
-        IndexSetInserter(ParallelIndexSet& indexSet, const GlobalIndex& component_size,
-                         std::size_t local_component_size, std::size_t num_components)
-            : indexSet_(&indexSet), component_size_(component_size),
-              local_component_size_(local_component_size),
-              num_components_(num_components)
-        {}
-        void operator()(const typename ParallelIndexSet::IndexPair& pair)
-        {
-            for(std::size_t i = 0; i < num_components_; i++)
-                indexSet_->add(i * component_size_ + pair.global(),
-                               LocalIndex(i * local_component_size_  + pair.local(),
-                                          pair.local().attribute()));
-        }
-    private:
-        ParallelIndexSet* indexSet_;
-        /// \brief The global number of unknowns per component/equation.
-        GlobalIndex component_size_;
-        /// \brief The local number of unknowns per component/equation.
-        std::size_t local_component_size_;
-        /// \brief The number of components/equations.
-        std::size_t num_components_;
-    };
     std::shared_ptr<ParallelIndexSet> indexSet_;
     std::shared_ptr<RemoteIndices> remoteIndices_;
     Dune::CollectiveCommunication<MPI_Comm> communicator_;
@@ -663,21 +429,6 @@ private:
 
 namespace Opm
 {
-/// \brief Extracts the information about the data decomposition from the grid for dune-istl
-/// 
-/// In the case that grid is a parallel grid this method will query it to get the information
-/// about the data decompoisition and convert it to the format expected by the linear algebra
-/// of dune-istl.
-/// \warn for UnstructuredGrid this function doesn't do anything.
-/// \param anyComm The handle to store the information in. If grid is a parallel grid
-/// then this will ecapsulate an instance of ParallelISTLInformation.
-/// \param grid The grid to inspect.
-
-inline void extractParallelGridInformationToISTL(std::any& anyComm, const UnstructuredGrid& grid)
-{
-    (void)anyComm; (void)grid;
-}
-
 /// \brief Accumulates entries masked with 1.
 /// \param container The container whose values to accumulate.
 /// \param maskContainer null pointer or a pointer to a container
@@ -687,19 +438,8 @@ inline void extractParallelGridInformationToISTL(std::any& anyComm, const Unstru
 template<class T1>
 auto
 accumulateMaskedValues(const T1& container, const std::vector<double>* maskContainer)
-    -> decltype(container[0]*(*maskContainer)[0])
-{
-    decltype(container[0]*(*maskContainer)[0]) initial = 0;
+    -> decltype(container[0]*(*maskContainer)[0]);
 
-    if( maskContainer )
-    {
-        return std::inner_product(container.begin(), container.end(), maskContainer->begin(),
-                                  initial);
-    }else
-    {
-        return std::accumulate(container.begin(), container.end(), initial);
-    }
-}   
 } // end namespace Opm
 
 #endif
diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp
index 2f73f4bde..62afcf60a 100644
--- a/opm/simulators/wells/RateConverter.hpp
+++ b/opm/simulators/wells/RateConverter.hpp
@@ -24,7 +24,6 @@
 
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/grid/utility/RegionMapping.hpp>
-#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 #include <opm/simulators/wells/RegionAttributeHelpers.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <dune/grid/common/gridenums.hh>
diff --git a/opm/simulators/wells/RegionAttributeHelpers.hpp b/opm/simulators/wells/RegionAttributeHelpers.hpp
index b7d2d67a5..57003f9c0 100644
--- a/opm/simulators/wells/RegionAttributeHelpers.hpp
+++ b/opm/simulators/wells/RegionAttributeHelpers.hpp
@@ -25,7 +25,6 @@
 
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/grid/utility/RegionMapping.hpp>
-#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 
 #include <dune/grid/common/gridenums.hh>
 #include <algorithm>
diff --git a/opm/simulators/wells/RegionAverageCalculator.hpp b/opm/simulators/wells/RegionAverageCalculator.hpp
index 3140e9a58..0c0bea16f 100644
--- a/opm/simulators/wells/RegionAverageCalculator.hpp
+++ b/opm/simulators/wells/RegionAverageCalculator.hpp
@@ -22,7 +22,6 @@
 
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/simulators/wells/RegionAttributeHelpers.hpp>
-#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
 #include <dune/grid/common/gridenums.hh>