diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index ceea210fa..a399903d8 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -150,6 +150,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
diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp
index e458fb0dd..43339d9b5 100644
--- a/opm/autodiff/ParallelOverlappingILU0.hpp
+++ b/opm/autodiff/ParallelOverlappingILU0.hpp
@@ -81,6 +81,79 @@ namespace Opm
 {
     namespace detail
     {
+
+    template<class M>
+    void milu0_decomposition(M& A,
+                             std::vector<typename M::block_type>* diagonal = nullptr)
+    {
+        using block = typename M::block_type;
+
+        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();
+
+            block sum_dropped;
+            sum_dropped = 0.0;
+
+            // 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
+                    {
+                        sum_dropped += modifier;
+                        ++a_kj;
+                    }
+                }
+            }
+
+            if ( a_ik.index() != irow.index() )
+                OPM_THROW(std::logic_error, "Matrix is missing diagonal for row " << irow.index());
+
+            *a_ik -= sum_dropped;
+
+            if ( diagonal )
+            {
+                diagonal->push_back(*a_ik);
+            }
+            a_ik->invert();   // compute inverse of diagonal block
+        }
+    }
+
+
       //! 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 )
diff --git a/tests/test_milu.cpp b/tests/test_milu.cpp
new file mode 100644
index 000000000..1b192f69c
--- /dev/null
+++ b/tests/test_milu.cpp
@@ -0,0 +1,165 @@
+#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());
+    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_ASSERT(point_difference < 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);
+    std::cout<< "Tested block size "<< bsize<<std::endl;
+}
+
+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>();
+}