From dcb3a20f7ab626682c3d982df087eda1a7cd1491 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Fri, 2 Dec 2016 15:34:31 +0100
Subject: [PATCH] Specialize direct solver for Dune::MatrixBlock.

In upstream they expect to be used with FielMatrix as
the block type only. Thus it it is impossible to use AMG
with a direct coarse solver within OPM. With this patch
it would be possible if we did not use field type float.
---
 opm/autodiff/ISTLSolver.hpp | 44 +++++++++++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp
index 591ecd231..5810dc51b 100644
--- a/opm/autodiff/ISTLSolver.hpp
+++ b/opm/autodiff/ISTLSolver.hpp
@@ -122,6 +122,50 @@ struct MatrixDimension< MatrixBlock< Scalar, n, m > >
 {
 };
 
+
+#if HAVE_UMFPACK
+
+/// \brief UMFPack specialization for MatrixBlock to make AMG happy
+///
+/// Without this the empty default implementation would be used.
+template<typename T, typename A, int n, int m>
+class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
+    : public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
+{
+    typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
+    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
+
+public:
+    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
+
+    UMFPack(const RealMatrix& matrix, int verbose, bool)
+        : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
+    {}
+};
+#endif
+
+#if HAVE_SUPERLU
+
+/// \brief SuperLU specialization for MatrixBlock to make AMG happy
+///
+/// Without this the empty default implementation would be used.
+template<typename T, typename A, int n, int m>
+class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
+    : public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
+{
+    typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
+    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
+
+public:
+    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
+
+    SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
+        : Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
+    {}
+};
+#endif
+
+
 } // end namespace Dune
 
 namespace Opm