diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index ceea210fa..8a9e6174e 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -142,6 +142,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_blackoil_amg.cpp
   tests/test_block.cpp
   tests/test_boprops_ad.cpp
+  tests/test_graphcoloring.cpp
   tests/test_rateconverter.cpp
   tests/test_span.cpp
   tests/test_syntax.cpp
@@ -150,6 +151,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_welldensitysegmented.cpp
   tests/test_vfpproperties.cpp
   tests/test_singlecellsolves.cpp
+  tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
   tests/test_multiphaseupwind.cpp
   tests/test_wellmodel.cpp
@@ -290,6 +292,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/FlowMainEbos.hpp
   opm/autodiff/FlowMainSequential.hpp
   opm/autodiff/GeoProps.hpp
+  opm/autodiff/GraphColoring.hpp
   opm/autodiff/GridHelpers.hpp
   opm/autodiff/GridInit.hpp
   opm/autodiff/ImpesTPFAAD.hpp
diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index 2a425f7b5..b565ad0d1 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -29,6 +29,7 @@
 
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
 
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -64,6 +65,9 @@ struct CPRParameter
     double cpr_relax_;
     double cpr_solver_tol_;
     int cpr_ilu_n_;
+    MILU_VARIANT cpr_ilu_milu_;
+    bool cpr_ilu_redblack_;
+    bool cpr_ilu_reorder_sphere_;
     int cpr_max_ell_iter_;
     bool cpr_use_amg_;
     bool cpr_use_bicgstab_;
@@ -77,25 +81,33 @@ struct CPRParameter
         // reset values to default
         reset();
 
-        cpr_relax_          = param.getDefault("cpr_relax", cpr_relax_);
-        cpr_solver_tol_     = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
-        cpr_ilu_n_          = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
-        cpr_max_ell_iter_   = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
-        cpr_use_amg_        = param.getDefault("cpr_use_amg", cpr_use_amg_);
-        cpr_use_bicgstab_   = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
-        cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
+        cpr_relax_                = param.getDefault("cpr_relax", cpr_relax_);
+        cpr_solver_tol_           = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
+        cpr_ilu_n_                = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
+        cpr_ilu_redblack_         = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
+        cpr_ilu_reorder_sphere_   = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
+        cpr_max_ell_iter_         = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
+        cpr_use_amg_              = param.getDefault("cpr_use_amg", cpr_use_amg_);
+        cpr_use_bicgstab_         = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
+        cpr_solver_verbose_       = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
         cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_);
+
+        std::string milu("ILU");
+        cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
     }
 
     void reset()
     {
-        cpr_relax_          = 1.0;
-        cpr_solver_tol_     = 1e-2;
-        cpr_ilu_n_          = 0;
-        cpr_max_ell_iter_   = 25;
-        cpr_use_amg_        = true;
-        cpr_use_bicgstab_   = true;
-        cpr_solver_verbose_ = false;
+        cpr_relax_                = 1.0;
+        cpr_solver_tol_           = 1e-2;
+        cpr_ilu_n_                = 0;
+        cpr_ilu_milu_             = MILU_VARIANT::ILU;
+        cpr_ilu_redblack_         = false;
+        cpr_ilu_reorder_sphere_   = true;
+        cpr_max_ell_iter_         = 25;
+        cpr_use_amg_              = true;
+        cpr_use_bicgstab_         = true;
+        cpr_solver_verbose_       = false;
         cpr_pressure_aggregation_ = false;
     }
 };
@@ -103,6 +115,31 @@ struct CPRParameter
 
 namespace ISTLUtility
 {
+
+template<class T>
+void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
+                      const CPRParameter& params)
+{
+    args.setN(params.cpr_ilu_n_);
+    args.setMilu(params.cpr_ilu_milu_);
+}
+
+template<class T>
+void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
+                      MILU_VARIANT milu, int n=0)
+{
+    args.setN(n);
+    args.setMilu(milu);
+}
+
+template<class S, class P>
+void setILUParameters(S&, const P&)
+{}
+
+template<class S, class P>
+void setILUParameters(S&, bool, int)
+{}
+
 ///
 /// \brief A traits class for selecting the types of the preconditioner.
 ///
@@ -167,9 +204,9 @@ struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
 //! \param relax The relaxation factor to use.
 template<class M, class X, class C>
 std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-createILU0Ptr(const M& A, const C& comm, double relax)
+createILU0Ptr(const M& A, const C& comm, double relax, MILU_VARIANT milu)
 {
-    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax);
+    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax, milu);
 }
 //! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
 //! \param A     The matrix of the linear system to solve.
@@ -177,25 +214,28 @@ createILU0Ptr(const M& A, const C& comm, double relax)
 //! \param relax The relaxation factor to use.
 template<class M, class X, class C>
 std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-createILUnPtr(const M& A, const C& comm, int ilu_n, double relax)
+createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, MILU_VARIANT milu)
 {
-    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax);
+    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax, milu );
 }
 /// \brief Creates the elliptic preconditioner (ILU0)
 /// \param Ae    The matrix of the elliptic system.
 /// \param relax The relaxation parameter for ILU0.
+/// \param milu  If true, the modified ilu approach is used. Dropped elements
+///              will get added to the diagonal of U to preserve the row sum
+///              for constant vectors (Ae = LUe).
 /// \param comm  The object describing the parallelization information and communication.
 template<class M, class X=typename M::range_type, class P>
 typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
 createEllipticPreconditionerPointer(const M& Ae, double relax,
-                                    const P& comm)
+                                    MILU_VARIANT milu, const P& comm)
 {
     typedef typename CPRSelector<M,X,X,P >
         ::EllipticPreconditioner ParallelPreconditioner;
 
     typedef typename CPRSelector<M,X,X,P>
         ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
-    return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
+    return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
 }
 
 template < class C, class Op, class P, class S, std::size_t index >
@@ -220,13 +260,14 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
     SmootherArgs  smootherArgs;
     smootherArgs.iterations = 1;
     smootherArgs.relaxationFactor = relax;
+    setILUParameters(smootherArgs, params);
 
     amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) );
 }
 
 template < class C, class Op, class P, class AMG >
 inline void
-createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr)
+createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr)
 {
     // TODO: revise choice of parameters
     int coarsenTarget=1200;
@@ -243,6 +284,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
     SmootherArgs  smootherArgs;
     smootherArgs.iterations = 1;
     smootherArgs.relaxationFactor = relax;
+    setILUParameters(smootherArgs, milu);
 
     amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
 }
@@ -254,7 +296,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
 //  \param amgPtr  The unique_ptr to be filled (return)
 template < int pressureIndex=0, class Op, class P, class AMG >
 inline void
-createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr )
+createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
 {
     // type of matrix
     typedef typename Op::matrix_type  M;
@@ -268,7 +310,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
     // The coarsening criterion used in the AMG
     typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
 
-    createAMGPreconditionerPointer<Criterion>(opA, relax, comm, amgPtr);
+    createAMGPreconditionerPointer<Criterion>(opA, relax, milu, comm, amgPtr);
 }
 
 } // end namespace ISTLUtility
@@ -370,10 +412,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
 
             // create the preconditioner for the whole system.
             if( param_.cpr_ilu_n_ == 0 ) {
-                pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_ );
+                pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_, param_.cpr_ilu_milu_ );
             }
             else {
-                pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_);
+                pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_);
             }
         }
 
@@ -531,11 +573,11 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
         {
             if( amg )
             {
-                ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ );
+                ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ );
             }
             else
             {
-                precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
+                precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm);
             }
        }
     };
diff --git a/opm/autodiff/GraphColoring.hpp b/opm/autodiff/GraphColoring.hpp
new file mode 100644
index 000000000..71164868e
--- /dev/null
+++ b/opm/autodiff/GraphColoring.hpp
@@ -0,0 +1,245 @@
+/*
+  Copyright 2018 Equinor
+
+  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/>.
+*/
+#ifndef OPM_GRAPHCOLORING_HEADER_INCLUDED
+#define OPM_GRAPHCOLORING_HEADER_INCLUDED
+
+#include <vector>
+#include <deque>
+#include <tuple>
+#include <algorithm>
+#include <numeric>
+#include <queue>
+#include <cstddef>
+
+namespace Opm
+{
+
+namespace Detail
+{
+template<class Graph>
+std::size_t colorGraphWelshPowell(const Graph& graph,
+                           std::deque<typename Graph::VertexDescriptor>& orderedVertices,
+                           std::vector<int>& colors,
+                           int color, int noVertices)
+{
+    std::vector<int> forbidden(noVertices, false);
+    std::size_t noColored = 0;
+
+    for(auto vertex = orderedVertices.begin(),
+            vertexEnd = orderedVertices.end();
+        vertex != vertexEnd; ++vertex)
+    {
+        // Skip forbidden vertices
+        while(vertex != vertexEnd && forbidden[*vertex])
+            ++vertex;
+        if ( vertex == vertexEnd )
+        {
+            break;
+        }
+
+        // Color Vertex
+        colors[*vertex] = color;
+        ++noColored;
+        // Forbid neighors
+        for(auto edge = graph.beginEdges(*vertex), endEdge = graph.endEdges(*vertex);
+            edge != endEdge; ++edge)
+        {
+            forbidden[edge.target()] = true;
+        }
+    }
+    // forbidden vertices will be colored next for coloring
+    using Vertex = typename Graph::VertexDescriptor;
+    auto newEnd = std::remove_if(orderedVertices.begin(), orderedVertices.end(),
+                                 [&forbidden](const Vertex& vertex)
+                                 {
+                                     return !forbidden[vertex];
+                                 });
+    orderedVertices.resize(newEnd-orderedVertices.begin());
+    return noColored;
+}
+template<class Graph, class Functor>
+std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescriptor root,
+                        Functor functor)
+{
+    std::vector<int> visited(graph.maxVertex() + 1, false);
+    using Vertex = typename Graph::VertexDescriptor;
+    std::queue<Vertex> nextVertices;
+    std::size_t noVisited = 0;
+    nextVertices.push(root);
+    visited[root] = true; // We do not visit root.
+
+    while( !nextVertices.empty() )
+    {
+        auto current = nextVertices.front();
+        for(auto edge = graph.beginEdges(current),
+                endEdge = graph.endEdges(current);
+            edge != endEdge; ++edge)
+        {
+            if ( ! visited[edge.target()] )
+            {
+                visited[edge.target()] = true;
+                nextVertices.push(edge.target());
+                functor(edge.target());
+                ++noVisited;
+            }
+        }
+        nextVertices.pop();
+    }
+    return noVisited;
+}
+} // end namespace Detail
+
+
+/// \brief Color the vertices of graph.
+///
+/// It uses the algorithm of Welsh and Powell for this.
+/// \param graph The graph to color. Must adhere to the graph interface of dune-istl.
+/// \return A pair of a vector with the colors of the vertices and the number of colors
+///         assigned
+template<class Graph>
+std::tuple<std::vector<int>, int, std::vector<std::size_t> >
+colorVerticesWelshPowell(const Graph& graph)
+{
+    using Vertex = typename Graph::VertexDescriptor;
+    std::deque<Vertex> orderedVertices;
+    auto noVertices = graph.maxVertex()+1;
+    std::vector<int> degrees(noVertices, 0);
+    int maxDegree = 0;
+    std::ptrdiff_t firstDegreeChange = 0;
+
+    // populate deque
+    for( auto vertex = graph.begin(), endVertex = graph.end();
+         vertex != endVertex; ++vertex)
+    {
+        auto currentVertex = *vertex;
+        auto& degree = degrees[currentVertex];
+
+        for(auto edge = graph.beginEdges(currentVertex),
+                endEdge = graph.endEdges(currentVertex);
+            edge != endEdge; ++edge)
+        {
+            ++degree;
+        }
+
+
+        if( degree >= maxDegree )
+        {
+            orderedVertices.emplace_front(currentVertex);
+            ++firstDegreeChange;
+            if(degree > maxDegree)
+            {
+                firstDegreeChange = 1;
+                maxDegree = degree;
+            }
+        }
+        else
+        {
+            orderedVertices.emplace_back(currentVertex);
+        }
+    }
+
+    // order deque by descending degree
+    std::stable_sort(orderedVertices.begin() + firstDegreeChange,
+                     orderedVertices.end(),
+              [&degrees](const Vertex& v1, const Vertex& v2)
+              {
+                  return degrees[v1] > degrees[v2];
+              });
+
+    // Overwrite degree with color
+    auto& colors = degrees;
+    std::fill(colors.begin(), colors.end(), -1);
+
+    int color = 0;
+    std::vector<std::size_t> verticesPerColor;
+    verticesPerColor.reserve(10);
+
+    while(!orderedVertices.empty())
+    {
+        verticesPerColor
+            .push_back(Detail::colorGraphWelshPowell(graph, orderedVertices, colors,
+                                                     color++, noVertices));
+    }
+    return std::make_tuple(colors, color, verticesPerColor);
+}
+
+/// \! Reorder colored graph preserving order of vertices with the same color.
+template<class Graph>
+std::vector<std::size_t>
+reorderVerticesPreserving(const std::vector<int>& colors, int noColors,
+                          const std::vector<std::size_t>& verticesPerColor,
+                          const Graph& graph)
+{
+    std::vector<std::size_t> colorIndex(noColors, 0);
+    std::vector<std::size_t> indices(graph.maxVertex() + 1);
+    std::partial_sum(verticesPerColor.begin(),
+                     verticesPerColor.begin()+verticesPerColor.size() - 1,
+                     colorIndex.begin() + 1);
+
+    for(const auto& vertex: graph)
+    {
+        indices[vertex] = colorIndex[colors[vertex]]++;
+    }
+    return indices;
+}
+
+/// \! Reorder Vetrices in spheres
+template<class Graph>
+std::vector<std::size_t>
+reorderVerticesSpheres(const std::vector<int>& colors, int noColors,
+                       const std::vector<std::size_t>& verticesPerColor,
+                       const Graph& graph,
+                       typename Graph::VertexDescriptor root)
+{
+    std::vector<std::size_t> colorIndex(noColors, 0);
+    const auto notVisitedTag = std::numeric_limits<std::size_t>::max();
+    std::vector<std::size_t> indices(graph.maxVertex() + 1, notVisitedTag);
+    using Vertex = typename Graph::VertexDescriptor;
+    std::partial_sum(verticesPerColor.begin(),
+                     verticesPerColor.begin()+verticesPerColor.size() - 1,
+                     colorIndex.begin() + 1);
+    std::size_t noVisited = 0;
+    auto numberer = [&colorIndex, &colors, &indices](Vertex vertex)
+        {
+            indices[vertex] = colorIndex[colors[vertex]]++;
+        };
+
+    while ( noVisited < graph.maxVertex() + 1 )
+    {
+        numberer(root);
+        ++noVisited; //root node already visited and not visited in BFS
+        noVisited += Detail::breadthFirstSearch(graph, root, numberer);
+        if ( noVisited < graph.maxVertex() + 1 )
+        {
+            // Graph is disconnected search for not yet visited node
+            for(auto vertex: graph)
+            {
+                if ( indices[vertex] ==  notVisitedTag )
+                {
+                    // \todo make sure that this is a peripheral node!
+                    root = vertex;
+                    break;
+                }
+            }
+        }
+    }
+    return indices;
+}
+} // end namespace Opm
+#endif
diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp
index dc31183fb..7692dbbdd 100644
--- a/opm/autodiff/ISTLSolver.hpp
+++ b/opm/autodiff/ISTLSolver.hpp
@@ -440,6 +440,7 @@ namespace Detail
                 }
 
                 const double relax = parameters_.ilu_relaxation_;
+                const MILU_VARIANT ilu_milu  = parameters_.ilu_milu_;
                 if (  parameters_.use_cpr_ )
                 {
                     using Matrix         = typename MatrixOperator::matrix_type;
@@ -452,7 +453,7 @@ namespace Detail
                     std::unique_ptr< AMG > amg;
                     // Construct preconditioner.
                     Criterion crit(15, 2000);
-                    constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax );
+                    constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
 
                     // Solve.
                     solve(linearOperator, x, istlb, *sp, *amg, result);
@@ -463,7 +464,7 @@ namespace Detail
                     std::unique_ptr< AMG > amg;
 
                     // Construct preconditioner.
-                    constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
+                    constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
 
                     // Solve.
                     solve(linearOperator, x, istlb, *sp, *amg, result);
@@ -495,7 +496,10 @@ namespace Detail
         {
             const double relax   = parameters_.ilu_relaxation_;
             const int ilu_fillin = parameters_.ilu_fillin_level_;
-            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax));
+            const MILU_VARIANT ilu_milu  = parameters_.ilu_milu_;
+            const bool ilu_redblack = parameters_.ilu_redblack_;
+            const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
+            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
             return precond;
         }
 
@@ -517,38 +521,46 @@ namespace Detail
         {
             typedef std::unique_ptr<ParPreconditioner> Pointer;
             const double relax  = parameters_.ilu_relaxation_;
-            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
+            const MILU_VARIANT ilu_milu  = parameters_.ilu_milu_;
+            const bool ilu_redblack = parameters_.ilu_redblack_;
+            const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
+            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
         }
 #endif
 
         template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
+        constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
+            ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, comm, amg );
         }
 
 
         template <class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
+        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax,
+                            const MILU_VARIANT milu) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
+            ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax,
+                                                                                 milu, comm, amg );
         }
 
         template <class C, class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
+        constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax,
+                            const MILU_VARIANT milu ) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg, parameters_ );
+            ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax,
+                                                                     comm, amg, parameters_ );
         }
 
 
         template <class C, class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
+        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const MILU_VARIANT milu ) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg, parameters_ );
+            ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, milu,
+                                                                     comm, amg, parameters_ );
         }
         /// \brief Solve the system using the given preconditioner and scalar product.
         template <class Operator, class ScalarProd, class Precond>
diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
index 377aecd1d..683993bfc 100644
--- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
@@ -28,6 +28,7 @@
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
 
 #include <array>
 #include <memory>
@@ -44,6 +45,9 @@ namespace Opm
         int    linear_solver_restart_;
         int    linear_solver_verbosity_;
         int    ilu_fillin_level_;
+        Opm::MILU_VARIANT   ilu_milu_;
+        bool   ilu_redblack_;
+        bool   ilu_reorder_sphere_;
         bool   newton_use_gmres_;
         bool   require_full_sparsity_pattern_;
         bool   ignoreConvergenceFailure_;
@@ -69,6 +73,10 @@ namespace Opm
             linear_solver_use_amg_    = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
             ilu_relaxation_           = param.getDefault("ilu_relaxation", ilu_relaxation_ );
             ilu_fillin_level_         = param.getDefault("ilu_fillin_level",  ilu_fillin_level_ );
+            ilu_redblack_             = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
+            ilu_reorder_sphere_       = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
+            std::string milu("ILU");
+            ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
 
             // Check whether to use cpr approach
             const std::string cprSolver = "cpr";
@@ -89,6 +97,9 @@ namespace Opm
             linear_solver_use_amg_    = false;
             ilu_fillin_level_         = 0;
             ilu_relaxation_           = 0.9;
+            ilu_milu_                 = MILU_VARIANT::ILU;
+            ilu_redblack_             = false;
+            ilu_reorder_sphere_       = true;
         }
     };
 
diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp
index e458fb0dd..4726ea12e 100644
--- a/opm/autodiff/ParallelOverlappingILU0.hpp
+++ b/opm/autodiff/ParallelOverlappingILU0.hpp
@@ -20,14 +20,20 @@
 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
 
+#include <opm/autodiff/GraphColoring.hpp>
 #include <opm/common/Exceptions.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <dune/common/version.hh>
 #include <dune/istl/preconditioner.hh>
 #include <dune/istl/paamg/smoother.hh>
+#include <dune/istl/paamg/graph.hh>
 #include <dune/istl/paamg/pinfo.hh>
 
 #include <type_traits>
+#include <numeric>
+#include <limits>
+#include <cstddef>
+#include <string>
 
 namespace Opm
 {
@@ -37,6 +43,64 @@ namespace Opm
 template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
 class ParallelOverlappingILU0;
 
+enum class MILU_VARIANT{
+    /// \brief Do not perform modified ILU
+    ILU = 0,
+    /// \brief \f$U_{ii} = U_{ii} +\f$  sum(dropped entries)
+    MILU_1 = 1,
+    /// \brief \f$U_{ii} = U_{ii} + sign(U_{ii}) * \f$ sum(dropped entries)
+    MILU_2 = 2,
+    /// \brief \f$U_{ii} = U_{ii} sign(U_{ii}) * \f$ sum(|dropped entries|)
+    MILU_3 = 3,
+    /// \brief \f$U_{ii} = U_{ii} + (U_{ii}>0?1:0) * \f$ sum(dropped entries)
+    MILU_4 = 4
+};
+
+inline MILU_VARIANT convertString2Milu(std::string milu)
+{
+    if( 0 == milu.compare("MILU_1") )
+    {
+        return MILU_VARIANT::MILU_1;
+    }
+    if ( 0 == milu.compare("MILU_2") )
+    {
+        return MILU_VARIANT::MILU_2;
+    }
+    if ( 0 == milu.compare("MILU_3") )
+    {
+        return MILU_VARIANT::MILU_3;
+    }
+    return MILU_VARIANT::ILU;
+}
+
+template<class F>
+class ParallelOverlappingILU0Args
+    : public Dune::Amg::DefaultSmootherArgs<F>
+{
+ public:
+    ParallelOverlappingILU0Args(MILU_VARIANT milu = MILU_VARIANT::ILU )
+        : milu_(milu)
+    {}
+    void setMilu(MILU_VARIANT milu)
+    {
+        milu_ = milu;
+    }
+    MILU_VARIANT getMilu() const
+    {
+        return milu_;
+    }
+    void setN(int n)
+    {
+        n_ = n;
+    }
+    int getN() const
+    {
+        return n_;
+    }
+ private:
+    MILU_VARIANT milu_;
+    int n_;
+};
 } // end namespace Opm
 
 namespace Dune
@@ -45,6 +109,13 @@ namespace Dune
 namespace Amg
 {
 
+
+template<class M, class X, class Y, class C>
+struct SmootherTraits<Opm::ParallelOverlappingILU0<M,X,Y,C> >
+{
+    using Arguments = Opm::ParallelOverlappingILU0Args<typename M::field_type>;
+};
+
 /// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother
 /// \tparam Matrix The type of the Matrix.
 /// \tparam Domain The type of the Vector representing the domain.
@@ -54,17 +125,18 @@ namespace Amg
 template<class Matrix, class Domain, class Range, class ParallelInfo>
 struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
 {
-    typedef Dune::SeqILU0<Matrix,Domain,Range> T;
+    typedef Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> T;
     typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
-    typedef ConstructionTraits<T> SeqConstructionTraits;
     static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
     {
-        return new Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>(args.getMatrix(),
-                                                         args.getComm(),
-                                                         args.getArgs().relaxationFactor);
+        return new T(args.getMatrix(),
+                     args.getComm(),
+                     args.getArgs().getN(),
+                     args.getArgs().relaxationFactor,
+                     args.getArgs().getMilu());
     }
 
-    static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
+    static inline void deconstruct(T* bp)
     {
         delete bp;
     }
@@ -81,6 +153,278 @@ namespace Opm
 {
     namespace detail
     {
+    struct Reorderer
+    {
+        virtual std::size_t operator[](std::size_t i) const = 0;
+    };
+
+    struct NoReorderer : public Reorderer
+    {
+        virtual std::size_t operator[](std::size_t i) const
+        {
+            return i;
+        }
+    };
+
+    struct RealReorderer : public Reorderer
+    {
+        RealReorderer(const std::vector<std::size_t>& ordering)
+            : ordering_(&ordering)
+        {}
+        virtual std::size_t operator[](std::size_t i) const
+        {
+            return (*ordering_)[i];
+        }
+        const std::vector<std::size_t>* ordering_;
+    };
+
+    struct IdentityFunctor
+    {
+        template<class T>
+        T operator()(const T& t)
+        {
+            return t;
+        }
+    };
+
+    struct OneFunctor
+    {
+        template<class T>
+        T operator()(const T&)
+        {
+            return 1.0;
+        }
+    };
+    struct SignFunctor
+    {
+        template<class T>
+        double operator()(const T& t)
+        {
+            if ( t < 0.0 )
+            {
+                return -1;
+            }
+            else
+            {
+                return 1.0;
+            }
+        }
+    };
+
+    struct IsPositiveFunctor
+    {
+        template<class T>
+        double operator()(const T& t)
+        {
+            if ( t < 0.0 )
+            {
+                return 0;
+            }
+            else
+            {
+                return 1;
+            }
+        }
+    };
+    struct AbsFunctor
+    {
+        template<class T>
+        T operator()(const T& t)
+        {
+            using std::abs;
+            return abs(t);
+        }
+    };
+
+    template<class M, class F1=detail::IdentityFunctor, class F2=detail::OneFunctor >
+    void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(),
+                             std::vector<typename M::block_type>* diagonal = nullptr)
+    {
+        if( diagonal )
+        {
+            diagonal->reserve(A.N());
+        }
+
+        for ( auto irow = A.begin(), iend = A.end(); irow != iend; ++irow)
+        {
+            auto a_i_end = irow->end();
+            auto a_ik    = irow->begin();
+
+            std::array<typename M::field_type, M::block_type::rows> sum_dropped{};
+
+            // Eliminate entries in lower triangular matrix
+            // and store factors for L
+            for ( ; a_ik.index() < irow.index(); ++a_ik )
+            {
+                auto k = a_ik.index();
+                auto a_kk = A[k].find(k);
+                // L_ik = A_kk^-1 * A_ik
+                a_ik->rightmultiply(*a_kk);
+
+                // modify the rest of the row, everything right of a_ik
+                // a_i* -=a_ik * a_k*
+                auto a_k_end = A[k].end();
+                auto a_kj = a_kk, a_ij = a_ik;
+                ++a_kj; ++a_ij;
+
+                while ( a_kj != a_k_end)
+                {
+                    auto modifier = *a_kj;
+                    modifier.leftmultiply(*a_ik);
+
+                    while( a_ij != a_i_end && a_ij.index() < a_kj.index())
+                    {
+                        ++a_ij;
+                    }
+
+                    if ( a_ij != a_i_end && a_ij.index() == a_kj.index() )
+                    {
+                        // Value is not dropped
+                        *a_ij -= modifier;
+                        ++a_ij; ++a_kj;
+                    }
+                    else
+                    {
+                        auto entry = sum_dropped.begin();
+                        for( const auto& row: modifier )
+                        {
+                            for( const auto& colEntry: row )
+                            {
+                                *entry += absFunctor(-colEntry);
+                            }
+                            ++entry;
+                        }
+                        ++a_kj;
+                    }
+                }
+            }
+
+            if ( a_ik.index() != irow.index() )
+                OPM_THROW(std::logic_error, "Matrix is missing diagonal for row " << irow.index());
+
+            int index = 0;
+            for(const auto& entry: sum_dropped)
+            {
+                auto& bdiag = (*a_ik)[index][index];
+                if(entry<0)
+                    std::cout << entry << std::endl;
+                bdiag += signFunctor(bdiag) * entry;
+                ++index;
+            }
+
+            if ( diagonal )
+            {
+                diagonal->push_back(*a_ik);
+            }
+            a_ik->invert();   // compute inverse of diagonal block
+        }
+    }
+
+    template<class M>
+    void milu0_decomposition(M& A,
+                             std::vector<typename M::block_type>* diagonal)
+    {
+        milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(),
+                            diagonal);
+    }
+
+    template<class M>
+    void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
+                             Reorderer& ordering, Reorderer& inverseOrdering)
+    {
+        using Map = std::map<std::size_t, int>;
+
+        auto iluRow = ILU.createbegin();
+
+        for(std::size_t i = 0, iend = A.N(); i < iend; ++i)
+        {
+            auto& orow = A[inverseOrdering[i]];
+
+            Map rowPattern;
+            for ( auto col = orow.begin(), cend = orow.end(); col != cend; ++col)
+            {
+                rowPattern[ordering[col.index()]] = 0;
+            }
+
+            for(auto ik = rowPattern.begin(); ik->first < i; ++ik)
+            {
+                if ( ik->second < n )
+                {
+                    auto& rowk = ILU[ik->first];
+
+                    for ( auto kj = rowk.find(ik->first), endk = rowk.end();
+                          kj != endk; ++kj)
+                    {
+                        // Assume double and block_type FieldMatrix
+                        // first element is misused to store generation number
+                        int generation = (*kj)[0][0];
+                        if(generation < n)
+                        {
+                            auto ij = rowPattern.find(kj.index());
+                            if ( ij == rowPattern.end() )
+                            {
+                                rowPattern[ordering[kj.index()]] = generation + 1;
+                            }
+                        }
+                    }
+                }
+            }
+            // create the row
+            for(const auto entry: rowPattern)
+            {
+                iluRow.insert(entry.first);
+            }
+            ++iluRow;
+
+            // write generation to newly created row.
+            auto generationPair = rowPattern.begin();
+            for ( auto col = ILU[i].begin(), cend = ILU[i].end(); col != cend;
+                  ++col, ++generationPair)
+            {
+                assert(col.index() == generationPair->first);
+                (*col)[0][0] = generationPair->second;
+            }
+        }
+
+        // copy Entries from A
+        for(auto iter=A.begin(), iend = A.end(); iter != iend; ++iter)
+        {
+            auto& newRow = ILU[ordering[iter.index()]];
+            // reset stored generation
+            for ( auto& col: newRow)
+            {
+                col = 0;
+            }
+            // copy row.
+            for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
+            {
+                newRow[ordering[col.index()]] = *col;
+            }
+        }
+        // call decomposition on pattern        
+        switch ( milu )
+        {
+        case MILU_VARIANT::MILU_1:
+            detail::milu0_decomposition ( ILU);
+            break;
+        case MILU_VARIANT::MILU_2:
+            detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
+                                          detail::SignFunctor() );
+            break;
+        case MILU_VARIANT::MILU_3:
+            detail::milu0_decomposition ( ILU, detail::AbsFunctor(),
+                                          detail::SignFunctor() );
+            break;
+        case MILU_VARIANT::MILU_4:
+            detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
+                                          detail::IsPositiveFunctor() );
+            break;
+        default:
+            bilu0_decomposition( ILU );
+            break;
+        }
+    }
+
       //! compute ILU decomposition of A. A is overwritten by its decomposition
       template<class M, class CRS, class InvVector>
       void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
@@ -154,6 +498,7 @@ namespace Opm
       }
     } // end namespace detail
 
+
 /// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as
 ///
 /// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with
@@ -257,10 +602,18 @@ public:
       \param A The matrix to operate on.
       \param n ILU fill in level (for testing). This does not work in parallel.
       \param w The relaxation factor.
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \param redblack Whether to use a red-black ordering.
+      \param reorder_sphere If true, we start the reordering at a root node.
+                            The vertices on each layer aound it (same distance) are
+                            ordered consecutivly. If false, we preserver the order of
+                            the vertices with the same color.
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
-                             const int n, const field_type w )
+                             const int n, const field_type w,
+                             MILU_VARIANT milu, bool redblack=false,
+                             bool reorder_sphere=true)
         : lower_(),
           upper_(),
           inv_(),
@@ -269,7 +622,8 @@ public:
     {
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
-        init( reinterpret_cast<const Matrix&>(A), n );
+        init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
+              reorder_sphere  );
     }
 
     /*! \brief Constructor gets all parameters to operate the prec.
@@ -277,10 +631,18 @@ public:
       \param comm   communication object, e.g. Dune::OwnerOverlapCopyCommunication
       \param n ILU fill in level (for testing). This does not work in parallel.
       \param w The relaxation factor.
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \param redblack Whether to use a red-black ordering.
+      \param reorder_sphere If true, we start the reordering at a root node.
+                            The vertices on each layer aound it (same distance) are
+                            ordered consecutivly. If false, we preserver the order of
+                            the vertices with the same color.
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
-                             const ParallelInfo& comm, const int n, const field_type w )
+                             const ParallelInfo& comm, const int n, const field_type w,
+                             MILU_VARIANT milu, bool redblack=false,
+                             bool reorder_sphere=true)
         : lower_(),
           upper_(),
           inv_(),
@@ -289,7 +651,8 @@ public:
     {
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
-        init( reinterpret_cast<const Matrix&>(A), n );
+        init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
+              reorder_sphere );
     }
 
     /*! \brief Constructor.
@@ -297,11 +660,18 @@ public:
       Constructor gets all parameters to operate the prec.
       \param A The matrix to operate on.
       \param w The relaxation factor.
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \param redblack Whether to use a red-black ordering.
+      \param reorder_sphere If true, we start the reordering at a root node.
+                  The vertices on each layer aound it (same distance) are
+                  ordered consecutivly. If false, we preserver the order of
+                  the vertices with the same color.
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
-                             const field_type w)
-        : ParallelOverlappingILU0( A, 0, w )
+                             const field_type w, MILU_VARIANT milu, bool redblack=false,
+                             bool reorder_sphere=true)
+        : ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
     {
     }
 
@@ -311,10 +681,18 @@ public:
       \param A      The matrix to operate on.
       \param comm   communication object, e.g. Dune::OwnerOverlapCopyCommunication
       \param w      The relaxation factor.
+      \param milu   The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \param redblack Whether to use a red-black ordering.
+      \param reorder_sphere If true, we start the reordering at a root node.
+                            The vertices on each layer aound it (same distance) are
+                            ordered consecutivly. If false, we preserver the order of
+                            the vertices with the same color.
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
-                             const ParallelInfo& comm, const field_type w)
+                             const ParallelInfo& comm, const field_type w,
+                             MILU_VARIANT milu, bool redblack=false,
+                             bool reorder_sphere=true)
         : lower_(),
           upper_(),
           inv_(),
@@ -323,7 +701,8 @@ public:
     {
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
-        init( reinterpret_cast<const Matrix&>(A), 0 );
+        init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
+              reorder_sphere );
     }
 
     /*!
@@ -344,7 +723,8 @@ public:
     */
     virtual void apply (Domain& v, const Range& d)
     {
-        Range& md = const_cast<Range&>(d);
+        Range& md = reorderD(d);
+        Domain& mv = reorderV(v);
         copyOwnerToAll( md );
 
         // iterator types
@@ -361,41 +741,42 @@ public:
         // lower triangular solve
         for( size_type i=0; i<iEnd; ++ i )
         {
-          dblock rhs( d[ i ] );
+          dblock rhs( md[ i ] );
           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 );
+            lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
           }
 
-          v[ i ] = rhs;  // Lii = I
+          mv[ i ] = rhs;  // Lii = I
         }
 
-        copyOwnerToAll( v );
+        copyOwnerToAll( mv );
 
         for( size_type i=0; i<iEnd; ++ i )
         {
-            vblock& vBlock = v[ lastRow - i ];
+            vblock& vBlock = mv[ lastRow - i ];
             vblock rhs ( vBlock );
             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 );
+                upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
             }
 
             // apply inverse and store result
             inv_[ i ].mv( rhs, vBlock);
         }
 
-        copyOwnerToAll( v );
+        copyOwnerToAll( mv );
 
         if( relaxation_ ) {
-            v *= w_;
+            mv *= w_;
         }
+        reorderBack(mv, v);
     }
 
     template <class V>
@@ -417,7 +798,7 @@ public:
     }
 
 protected:
-    void init( const Matrix& A, const int iluIteration )
+    void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
     {
         // (For older DUNE versions the communicator might be
         // invalid if redistribution in AMG happened on the coarset level.
@@ -441,17 +822,103 @@ protected:
 
         std::unique_ptr< Matrix > ILU;
 
+        if ( redBlack )
+        {
+            using Graph = Dune::Amg::MatrixGraph<const Matrix>;
+            Graph graph(A);
+            auto colorsTuple = colorVerticesWelshPowell(graph);
+            const auto& colors = std::get<0>(colorsTuple);
+            const auto& verticesPerColor = std::get<2>(colorsTuple);
+            auto noColors = std::get<1>(colorsTuple);
+            if ( reorderSpheres )
+            {
+                ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
+                                                   graph, 0);
+            }
+            else
+            {
+                ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
+                                                      graph);
+            }
+        }
+
+        std::vector<std::size_t> inverseOrdering(ordering_.size());
+        std::size_t index = 0;
+        for( auto newIndex: ordering_)
+        {
+            inverseOrdering[newIndex] = index++;
+        }
+
         try
         {
             if( iluIteration == 0 ) {
                 // create ILU-0 decomposition
-                ILU.reset( new Matrix( A ) );
-                bilu0_decomposition( *ILU );
+                if ( ordering_.empty() )
+                {
+                    ILU.reset( new Matrix( A ) );
+                }
+                else
+                {
+                    ILU.reset( new Matrix(A.N(), A.M(), A.nonzeroes(), Matrix::row_wise));
+                    auto& newA = *ILU;
+                    // Create sparsity pattern
+                    for(auto iter=newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
+                    {
+                        const auto& row = A[inverseOrdering[iter.index()]];
+                        for(auto col = row.begin(), cend = row.end(); col != cend; ++col)
+                        {
+                            iter.insert(ordering_[col.index()]);
+                        }
+                    }
+                    // Copy values
+                    for(auto iter = A.begin(), iend = A.end(); iter != iend; ++iter)
+                    {
+                        auto& newRow = newA[ordering_[iter.index()]];
+                        for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
+                        {
+                            newRow[ordering_[col.index()]] = *col;
+                        }
+                    }
+                }
+
+                switch ( milu )
+                {
+                case MILU_VARIANT::MILU_1:
+                    detail::milu0_decomposition ( *ILU);
+                    break;
+                case MILU_VARIANT::MILU_2:
+                    detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
+                                                  detail::SignFunctor() );
+                    break;
+                case MILU_VARIANT::MILU_3:
+                    detail::milu0_decomposition ( *ILU, detail::AbsFunctor(),
+                                                  detail::SignFunctor() );
+                    break;
+                case MILU_VARIANT::MILU_4:
+                    detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
+                                                  detail::IsPositiveFunctor() );
+                    break;
+                default:
+                    bilu0_decomposition( *ILU );
+                    break;
+                }
             }
             else {
                 // create ILU-n decomposition
                 ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
-                bilu_decomposition( A, iluIteration, *ILU );
+                std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
+                if ( ordering_.empty() )
+                {
+                    reorderer.reset(new detail::NoReorderer());
+                    inverseReorderer.reset(new detail::NoReorderer());
+                }
+                else
+                {
+                    reorderer.reset(new detail::RealReorderer(ordering_));
+                    inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
+                }
+
+                milun_decomposition( A, iluIteration, milu, *ILU, *reorderer, *inverseReorderer );
             }
         }
         catch ( Dune::MatrixBlockError error )
@@ -475,11 +942,70 @@ protected:
         detail::convertToCRS( *ILU, lower_, upper_, inv_ );
     }
 
+    /// \brief Reorder D if needed and return a reference to it.
+    Range& reorderD(const Range& d)
+    {
+        if ( ordering_.empty())
+        {
+            // As d is non-const in the apply method of the
+            // solver casting away constness in this particular
+            // setting is not undefined. It is ugly though but due
+            // to the preconditioner interface of dune-istl.
+            return const_cast<Range&>(d);
+        }
+        else
+        {
+            reorderedD_.resize(d.size());
+            std::size_t i = 0;
+            for(auto index: ordering_)
+            {
+                reorderedD_[index]=d[i++];
+            }
+            return reorderedD_;
+        }
+    }
+
+    /// \brief Reorder V if needed and return a reference to it.
+    Domain& reorderV(Domain& v)
+    {
+        if ( ordering_.empty())
+        {
+            return v;
+        }
+        else
+        {
+            reorderedV_.resize(v.size());
+            std::size_t i = 0;
+            for(auto index: ordering_)
+            {
+                reorderedV_[index]=v[i++];
+            }
+            return reorderedV_;
+        }
+    }
+
+    void reorderBack(const Range& reorderedV, Range& v)
+    {
+        if ( !ordering_.empty() )
+        {
+            std::size_t i = 0;
+            for(auto index: ordering_)
+            {
+                v[i++] = reorderedV[index];
+            }
+        }
+    }
 protected:
     //! \brief The ILU0 decomposition of the matrix.
     CRS lower_;
     CRS upper_;
     std::vector< block_type > inv_;
+    //! \brief the reordering of the unknowns
+    std::vector< std::size_t > ordering_;
+    //! \brief The reordered right hand side
+    Range reorderedD_;
+    //! \brief The reordered left hand side.
+    Domain reorderedV_;
 
     const ParallelInfo* comm_;
     //! \brief The relaxation factor to use.
diff --git a/tests/test_graphcoloring.cpp b/tests/test_graphcoloring.cpp
new file mode 100644
index 000000000..22b2d34dd
--- /dev/null
+++ b/tests/test_graphcoloring.cpp
@@ -0,0 +1,101 @@
+#include <config.h>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/paamg/graph.hh>
+
+#include <opm/autodiff/GraphColoring.hpp>
+
+#define BOOST_TEST_MODULE GraphColoringTest
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+///! \brief check that all indices are represented in the new ordering.
+void checkAllIndices(const std::vector<std::size_t>& ordering)
+{
+    std::vector<int> counters(ordering.size(), 0);
+    for(auto index: ordering)
+    {
+        ++counters[index];
+    }
+
+    for(auto count: counters)
+    {
+        BOOST_CHECK(count==1);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(TestWelschPowell)
+{
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+    using Graph = Dune::Amg::MatrixGraph<Matrix>;
+    int N = 10;
+    Matrix matrix(N*N, N*N, 5, 0.4, Matrix::implicit);
+    for( int j = 0; j < N; j++)
+    {
+        for(int i = 0; i < N; i++)
+        {
+            auto index = j*10+i;
+            matrix.entry(index,index) = 1;
+
+            if ( i > 0 )
+            {
+                matrix.entry(index,index-1) = 1;
+            }
+            if ( i  < N - 1)
+            {
+                matrix.entry(index,index+1) = 1;
+            }
+
+            if ( j > 0 )
+            {
+                matrix.entry(index,index-N) = 1;
+            }
+            if ( j  < N - 1)
+            {
+                matrix.entry(index,index+N) = 1;
+            }
+
+        }
+    }
+    matrix.compress();
+
+    Graph graph(matrix);
+    auto colorsTuple = Opm::colorVerticesWelshPowell(graph);
+    const auto& colors = std::get<0>(colorsTuple);
+    const auto& verticesPerColor = std::get<2>(colorsTuple);
+    auto noColors = std::get<1>(colorsTuple);
+    auto firstCornerColor = colors[0];
+    BOOST_CHECK(noColors == 2);
+
+    // Check for checkerboard coloring
+
+    for( int j = 0, index = 0; j < N; j++)
+    {
+        auto expectedColor = firstCornerColor;
+
+        for(int i = 0; i < N; i++)
+        {
+            BOOST_CHECK(colors[index]==expectedColor);
+            index++;
+            expectedColor = (expectedColor + 1) % 2;
+        }
+        firstCornerColor=(firstCornerColor + 1) % 2;
+    }
+    auto newOrder = Opm::reorderVerticesPreserving(colors, noColors, verticesPerColor,
+                                                   graph);
+    std::vector<std::size_t> colorIndex(noColors, 0);
+    std::partial_sum(verticesPerColor.begin(),
+                    verticesPerColor.begin()+verticesPerColor.size()-1,
+                    colorIndex.begin()+1);
+
+    for (auto vertex : graph)
+    {
+        BOOST_CHECK(colorIndex[colors[vertex]]++ == newOrder[vertex]);
+    }
+    checkAllIndices(newOrder);
+    newOrder = Opm::reorderVerticesSpheres(colors, noColors, verticesPerColor,
+                                           graph, 0);
+    checkAllIndices(newOrder);
+}
diff --git a/tests/test_milu.cpp b/tests/test_milu.cpp
new file mode 100644
index 000000000..4c3ebda73
--- /dev/null
+++ b/tests/test_milu.cpp
@@ -0,0 +1,193 @@
+#include<config.h>
+
+#define BOOST_TEST_MODULE MILU0Test
+
+#include<vector>
+#include<memory>
+
+#include<dune/istl/bcrsmatrix.hh>
+#include<dune/istl/bvector.hh>
+#include<dune/common/fmatrix.hh>
+#include<dune/common/fvector.hh>
+#include<opm/autodiff/ParallelOverlappingILU0.hpp>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+template<class M>
+void test_milu0(M& A)
+{
+    auto ILU = A;
+    
+    std::shared_ptr<std::vector<typename M::block_type>> diagonal = nullptr;
+    diagonal.reset(new std::vector<typename M::block_type>());
+
+    Opm::detail::milu0_decomposition(ILU, diagonal.get());
+#ifdef DEBUG
+    if ( A.N() < 11)
+    {
+        Dune::printmatrix(std::cout, ILU, "ILU", "row");
+        std::cout << "Diagonal: ";
+
+        for (const auto& d : *diagonal)
+            std::cout << d << " ";
+        std::cout<<std::endl;
+    }
+#endif
+    Dune::BlockVector<Dune::FieldVector<typename M::field_type, M::block_type::rows>> e(A.N()), x1(A.N()), x2(A.N()), t(A.N());
+    e = 1;
+    A.mv(e, x1);
+    // Compute L*U*e;
+    // t=Ue
+    
+    for ( auto irow = ILU.begin(), iend = ILU.end(); irow != iend; ++irow)
+    {
+        auto i = irow.index();
+        (*diagonal)[i].mv(e[0], t[i]);
+        auto col = irow->find(i);
+        auto colend = irow->end();
+
+        if ( col == colend )
+        {
+            OPM_THROW(std::logic_error, "Missing diagonal entry for row " << irow.index());
+        }
+        for (++col ;col != colend; ++col)
+        {
+            col->umv(e[0], t[i]);
+        }
+    }
+            
+    for ( auto irow = ILU.begin(), iend = ILU.end(); irow != iend; ++irow)
+    {
+        auto i = irow.index();
+        x2[i] = 0;
+        for (auto col = irow->begin(); col.index() < irow.index(); ++col)
+        {
+            col->umv(t[col.index()], x2[i]);
+        }
+        x2[i] += t[i];
+    }
+    auto diff = x2;
+    diff-=x1;
+    for ( std::size_t i = 0, end = A.N(); i < end; ++i)
+    {
+        auto point_difference = diff[i].two_norm();
+#ifdef DEBUG
+        std::cout<<"index "<<i<<" size "<<diff.size()<<" difference"<<point_difference<<std::endl;
+#endif
+        BOOST_CHECK(point_difference < 1e-12);
+    }
+
+    // Test that (LU)^-1Ae=e
+    A.mv(e, x1);
+    bilu_backsolve(ILU, x2, x1);
+    diff = x2;
+    diff -= e;
+
+    for ( std::size_t i = 0, end = A.N(); i < end; ++i)
+    {
+#ifdef DEBUG
+        auto point_difference = diff[i].two_norm();
+        std::cout<<"index "<<i<<" size "<<diff.size()<<" point_difference "<<point_difference<<std::endl;
+#endif
+        BOOST_CHECK_CLOSE(x2[i].two_norm(), e[i].two_norm(), 1e-12);
+    }
+}
+
+template<class B, class Alloc>
+void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
+{
+  typedef typename Dune::BCRSMatrix<B,Alloc> Matrix;
+  A.setSize(N*N, N*N, N*N*5);
+  A.setBuildMode(Matrix::row_wise);
+
+  for (typename Dune::BCRSMatrix<B,Alloc>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
+    int x = i.index()%N; // x coordinate in the 2d field
+    int y = i.index()/N;  // y coordinate in the 2d field
+
+    if(y>0)
+      // insert lower neighbour
+      i.insert(i.index()-N);
+    if(x>0)
+      // insert left neighbour
+      i.insert(i.index()-1);
+
+    // insert diagonal value
+    i.insert(i.index());
+
+    if(x<N-1)
+      //insert right neighbour
+      i.insert(i.index()+1);
+    if(y<N-1)
+      // insert upper neighbour
+      i.insert(i.index()+N);
+  }
+}
+
+
+template<class B, class Alloc>
+void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
+{
+  typedef typename Dune::BCRSMatrix<B,Alloc>::field_type FieldType;
+
+  setupSparsityPattern(A,N);
+
+  B diagonal(static_cast<FieldType>(0)), bone(static_cast<FieldType>(0)),
+  beps(static_cast<FieldType>(0));
+  for(typename B::RowIterator b = diagonal.begin(); b !=  diagonal.end(); ++b)
+    b->operator[](b.index())=4;
+
+
+  for(typename B::RowIterator b=bone.begin(); b !=  bone.end(); ++b)
+    b->operator[](b.index())=-1.0;
+
+
+  for (typename Dune::BCRSMatrix<B,Alloc>::RowIterator i = A.begin(); i != A.end(); ++i) {
+    int x = i.index()%N; // x coordinate in the 2d field
+    int y = i.index()/N;  // y coordinate in the 2d field
+    
+    i->operator[](i.index())=diagonal;
+    
+    if(y>0)
+        i->operator[](i.index()-N)=bone;
+
+    if(y<N-1)
+        i->operator[](i.index()+N)=bone;
+
+    if(x>0)
+        i->operator[](i.index()-1)=bone;
+
+    if(x < N-1)
+        i->operator[](i.index()+1)=bone;
+  }
+}
+
+template<int bsize>
+void test()
+{
+    std::size_t N = 32;
+    Dune::BCRSMatrix<Dune::FieldMatrix<double, bsize, bsize> > A;
+    setupLaplacian(A, N);
+    test_milu0(A);
+#ifdef DEBUG
+    std::cout<< "Tested block size "<< bsize<<std::endl;
+#endif
+}
+
+BOOST_AUTO_TEST_CASE(MILULaplace1)
+{
+    test<1>();
+}
+
+BOOST_AUTO_TEST_CASE(MILULaplace2)
+{
+    test<2>();
+}
+BOOST_AUTO_TEST_CASE(MILULaplace3)
+{
+    test<3>();
+}
+BOOST_AUTO_TEST_CASE(MILULaplace4)
+{
+    test<4>();
+}