diff --git a/CMakeLists.txt b/CMakeLists.txt
index 376337301bf..2e408bfdcdd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,6 +35,7 @@ option(USE_HYPRE "Use the Hypre library for linear solvers?" OFF)
 set(OPM_COMPILE_COMPONENTS "2;3;4;5;6;7;8;9;10" CACHE STRING "The components to compile support for")
 option(USE_OPENCL "Enable OpenCL support?" ON)
 option(USE_DEV_SIMULATOR_IN_TESTS "Use the development simulator binaries in tests?" OFF)
+option(USE_EXPERIMENTAL_ISTL_MATRIX "Use experimental IBCRSMatrix from dune-istl?" OFF)
 
 # Wrapper for opm_add_target_options that also adds
 # compile definitions and target links from the library
@@ -126,6 +127,9 @@ macro(opm-simulators_prereqs_hook)
     target_link_libraries(opmsimulators PUBLIC HDF5::HDF5)
     target_compile_definitions(opmsimulators PUBLIC HAVE_HDF5=1)
   endif()
+  if(USE_EXPERIMENTAL_ISTL_MATRIX)
+    target_compile_definitions(opmsimulators PUBLIC OPM_USE_EXPERIMENTAL_IBCRSMATRIX=1)
+  endif()
 
   if (OPM_ENABLE_PYTHON)
     # Be backwards compatible.
diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 40891042c65..2d0c31cb634 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -1095,6 +1095,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/amgcpr.hh
   opm/simulators/linalg/bicgstabsolver.hh
   opm/simulators/linalg/blacklist.hh
+  opm/simulators/linalg/BlockSparseMatrix.hpp
   opm/simulators/linalg/combinedcriterion.hh
   opm/simulators/linalg/convergencecriterion.hh
   opm/simulators/linalg/DILU.hpp
diff --git a/flowexperimental/comp/wells/CompWellEquations.hpp b/flowexperimental/comp/wells/CompWellEquations.hpp
index 2bb847e5717..dda50399c07 100644
--- a/flowexperimental/comp/wells/CompWellEquations.hpp
+++ b/flowexperimental/comp/wells/CompWellEquations.hpp
@@ -20,9 +20,10 @@
 #ifndef OPM_COMP_WELL_EQUATIONS_HPP
 #define OPM_COMP_WELL_EQUATIONS_HPP
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 namespace Opm {
@@ -48,11 +49,11 @@ class CompWellEquations
 
     // the matrix type for the diagonal matrix D
     using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq>;
-    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
+    using DiagMatWell = BlockSparseMatrix<DiagMatrixBlockWellType>;
 
     // the matrix type for the non-diagonal matrix B and C^T
     using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, numEq>;
-    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
+    using OffDiagMatWell = BlockSparseMatrix<OffDiagMatrixBlockWellType>;
 
     void init(const int num_conn, const std::vector<std::size_t>& cells);
 
diff --git a/opm/simulators/flow/GenericTemperatureModel.hpp b/opm/simulators/flow/GenericTemperatureModel.hpp
index 20719b0e9e6..4957bd73b50 100644
--- a/opm/simulators/flow/GenericTemperatureModel.hpp
+++ b/opm/simulators/flow/GenericTemperatureModel.hpp
@@ -28,7 +28,6 @@
 #ifndef OPM_GENERIC_TEMPERATURE_MODEL_HPP
 #define OPM_GENERIC_TEMPERATURE_MODEL_HPP
 
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/grid/common/CartesianIndexMapper.hpp>
 
@@ -37,6 +36,7 @@
 #include <opm/models/blackoil/blackoilmodel.hh>
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 
 #include <cstddef>
@@ -54,7 +54,7 @@ class GenericTemperatureModel
 public:
     // the jacobian matrix
     using MatrixBlockTemp = MatrixBlock<Scalar, 1, 1>;
-    using EnergyMatrix = Dune::BCRSMatrix<MatrixBlockTemp>;
+    using EnergyMatrix = BlockSparseMatrix<MatrixBlockTemp>;
     using EnergyVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static constexpr int dimWorld = Grid::dimensionworld;
diff --git a/opm/simulators/flow/GenericTracerModel.hpp b/opm/simulators/flow/GenericTracerModel.hpp
index a026856a0c7..14f844d5c52 100644
--- a/opm/simulators/flow/GenericTracerModel.hpp
+++ b/opm/simulators/flow/GenericTracerModel.hpp
@@ -28,14 +28,13 @@
 #ifndef OPM_GENERIC_TRACER_MODEL_HPP
 #define OPM_GENERIC_TRACER_MODEL_HPP
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <opm/grid/common/CartesianIndexMapper.hpp>
 
 #include <opm/input/eclipse/EclipseState/Phase.hpp>
 
 #include <opm/models/blackoil/blackoilmodel.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/wells/WellTracerRate.hpp>
 
@@ -56,7 +55,7 @@ template<class Grid, class GridView, class DofMapper, class Stencil, class Fluid
 class GenericTracerModel {
 public:
     using TracerVectorSingle = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
-    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
+    using TracerMatrix = BlockSparseMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
     using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar, 2>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static constexpr int dimWorld = Grid::dimensionworld;
diff --git a/opm/simulators/flow/NonlinearSolver.hpp b/opm/simulators/flow/NonlinearSolver.hpp
index 5d51e634866..c7439171383 100644
--- a/opm/simulators/flow/NonlinearSolver.hpp
+++ b/opm/simulators/flow/NonlinearSolver.hpp
@@ -22,7 +22,6 @@
 #define OPM_NON_LINEAR_SOLVER_HPP
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>
diff --git a/opm/simulators/linalg/AmgxPreconditioner.hpp b/opm/simulators/linalg/AmgxPreconditioner.hpp
index 43e39ff2341..8828d146c3c 100644
--- a/opm/simulators/linalg/AmgxPreconditioner.hpp
+++ b/opm/simulators/linalg/AmgxPreconditioner.hpp
@@ -28,7 +28,6 @@
 #include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <amgx_c.h>
 
diff --git a/opm/simulators/linalg/BlockSparseMatrix.hpp b/opm/simulators/linalg/BlockSparseMatrix.hpp
new file mode 100644
index 00000000000..9afc922cb1f
--- /dev/null
+++ b/opm/simulators/linalg/BlockSparseMatrix.hpp
@@ -0,0 +1,50 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ */
+#ifndef OPM_BLOCKSPARSEMATRIX_HEADER_INCLUDED
+#define OPM_BLOCKSPARSEMATRIX_HEADER_INCLUDED
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+#include <dune/istl/ibcrsmatrix.hh>
+#else
+#include <dune/istl/bcrsmatrix.hh>
+#endif
+
+namespace Opm {
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+// Use the exprimental IBCRSMatrix, which has a better performance under distributed workloads.
+template<class B, class A = std::allocator<B>>
+using BlockSparseMatrix = Dune::IBCRSMatrix<B, uint_least32_t, A>;
+#else
+// Use the classic BCRSMatrix.
+template<class B, class A = std::allocator<B>>
+using BlockSparseMatrix = Dune::BCRSMatrix<B, A>;
+#endif
+
+} // namespace Opm
+
+#endif
diff --git a/opm/simulators/linalg/DILU.hpp b/opm/simulators/linalg/DILU.hpp
index 7e4ce5f7a6d..43130d918f3 100644
--- a/opm/simulators/linalg/DILU.hpp
+++ b/opm/simulators/linalg/DILU.hpp
@@ -24,7 +24,6 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/unused.hh>
 #include <dune/common/version.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/simulators/linalg/GraphColoring.hpp>
 
diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index 9910256c63c..5b4ab0b9874 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -27,6 +27,7 @@
 #include <opm/simulators/linalg/ilufirstelement.hh>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/Preconditioner2InverseOperator.hpp>
 #include <opm/simulators/linalg/WellOperators.hpp>
@@ -355,7 +356,7 @@ namespace Dune
 template<class Scalar, int N>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar, N>>;
 template<class Scalar, int N>
-using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, N, N>>;
+using OBM = Opm::BlockSparseMatrix<Opm::MatrixBlock<Scalar, N, N>>;
 
 // Sequential operators.
 template<class Scalar, int N>
diff --git a/opm/simulators/linalg/HyprePreconditioner.hpp b/opm/simulators/linalg/HyprePreconditioner.hpp
index c5c1421abac..4fb405313b7 100644
--- a/opm/simulators/linalg/HyprePreconditioner.hpp
+++ b/opm/simulators/linalg/HyprePreconditioner.hpp
@@ -34,7 +34,6 @@
 #endif
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <HYPRE.h>
 #include <HYPRE_krylov.h>
diff --git a/opm/simulators/linalg/ISTLSolver.cpp b/opm/simulators/linalg/ISTLSolver.cpp
index c7d7031a25b..1924df3fd06 100644
--- a/opm/simulators/linalg/ISTLSolver.cpp
+++ b/opm/simulators/linalg/ISTLSolver.cpp
@@ -29,6 +29,7 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <fmt/format.h>
@@ -164,7 +165,7 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
 }
 
 template<class Scalar, int Dim>
-using BM = Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>;
+using BM = BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>;
 template<class Scalar, int Dim>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
 
diff --git a/opm/simulators/linalg/ISTLSolverGpuBridge.cpp b/opm/simulators/linalg/ISTLSolverGpuBridge.cpp
index c00b5ce1110..ddaa38f7740 100644
--- a/opm/simulators/linalg/ISTLSolverGpuBridge.cpp
+++ b/opm/simulators/linalg/ISTLSolverGpuBridge.cpp
@@ -29,6 +29,7 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <fmt/format.h>
@@ -237,7 +238,7 @@ copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
 }
 
 template<class Scalar, int Dim>
-using BM = Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>;
+using BM = BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>;
 template<class Scalar, int Dim>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
 
diff --git a/opm/simulators/linalg/MILU.cpp b/opm/simulators/linalg/MILU.cpp
index c1a32556f6c..85309823200 100644
--- a/opm/simulators/linalg/MILU.cpp
+++ b/opm/simulators/linalg/MILU.cpp
@@ -23,13 +23,13 @@
 
 #include <dune/common/version.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/ilu.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <array>
 #include <algorithm>
@@ -267,8 +267,8 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
     INSTANTIATE_ILUN(__VA_ARGS__)
 
 #define INSTANTIATE_DIM(T,Dim)                                         \
-    INSTANTIATE_FULL(T,Dune::BCRSMatrix<MatrixBlock<T,Dim,Dim>>)       \
-    INSTANTIATE_FULL(T,Dune::BCRSMatrix<Dune::FieldMatrix<T,Dim,Dim>>)
+    INSTANTIATE_FULL(T,BlockSparseMatrix<MatrixBlock<T,Dim,Dim>>)       \
+    INSTANTIATE_FULL(T,BlockSparseMatrix<Dune::FieldMatrix<T,Dim,Dim>>)
 
 #define INSTANTIATE_TYPE(T)                 \
     template T identityFunctor(const T&);   \
diff --git a/opm/simulators/linalg/MatrixMarketSpecializations.hpp b/opm/simulators/linalg/MatrixMarketSpecializations.hpp
index 2743b4410fd..12e69fb7499 100644
--- a/opm/simulators/linalg/MatrixMarketSpecializations.hpp
+++ b/opm/simulators/linalg/MatrixMarketSpecializations.hpp
@@ -44,6 +44,18 @@ namespace MatrixMarketImpl
         }
     };
 
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+    template <typename T, typename I, int i, int j, typename A>
+    struct mm_header_printer<IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>>
+    {
+        static void print(std::ostream& os)
+        {
+            os << "%%MatrixMarket matrix coordinate ";
+            os << mm_numeric_type<T>::str() << " general" << std::endl;
+        }
+    };
+#endif
+
     template <typename T, int i, int j, typename A>
     struct mm_block_structure_header<BCRSMatrix<Opm::MatrixBlock<T,i,j>, A>>
     {
@@ -54,6 +66,19 @@ namespace MatrixMarketImpl
             os << i << " " << j << std::endl;
         }
     };
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+    template <typename T, typename I, int i, int j, typename A>
+    struct mm_block_structure_header<IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>>
+    {
+        using M = IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>;
+        static void print(std::ostream& os, const M&)
+        {
+            os << "% ISTL_STRUCT blocked ";
+            os << i << " " << j << std::endl;
+        }
+    };
+#endif
 } // namespace MatrixMarketImpl
 
 namespace MatrixMarketImpl
diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
index 63145c73e8b..9c630d183e6 100644
--- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
+++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
@@ -27,7 +27,6 @@
 #include <opm/common/ErrorMacros.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/paamg/amg.hh>
 
 #include <fstream>
diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.cpp b/opm/simulators/linalg/ParallelOverlappingILU0.cpp
index 47679303bf1..d77128c394c 100644
--- a/opm/simulators/linalg/ParallelOverlappingILU0.cpp
+++ b/opm/simulators/linalg/ParallelOverlappingILU0.cpp
@@ -22,17 +22,19 @@
 
 #include <opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <dune/istl/owneroverlapcopy.hh>
 
 namespace Opm
 {
 
 #define INSTANTIATE_PAR(T, Dim, ...)                                                     \
-  template class ParallelOverlappingILU0<Dune::BCRSMatrix<MatrixBlock<T,Dim,Dim>>,       \
+  template class ParallelOverlappingILU0<BlockSparseMatrix<MatrixBlock<T,Dim,Dim>>, \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          __VA_ARGS__>;                                   \
-  template class ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::FieldMatrix<T,Dim,Dim>>, \
+  template class ParallelOverlappingILU0<BlockSparseMatrix<Dune::FieldMatrix<T,Dim,Dim>>, \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          __VA_ARGS__>;
diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp
index deb84fbc409..c455e890772 100644
--- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp
+++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp
@@ -25,6 +25,7 @@
 
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/DILU.hpp>
 #include <opm/simulators/linalg/ExtraSmoothers.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
@@ -198,21 +199,21 @@ PreconditionerFactory<Operator, Comm>::addCreator(const std::string& type, ParCr
 using CommSeq = Dune::Amg::SequentialInformation;
 
 template<class Scalar, int Dim>
-using OpFSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
+using OpFSeq = Dune::MatrixAdapter<BlockSparseMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 template<class Scalar, int Dim>
-using OpBSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpBSeq = Dune::MatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 
 template<class Scalar, int Dim, bool overlap>
-using OpW = WellModelMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpW = WellModelMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 
 template<class Scalar, int Dim, bool overlap>
-using OpWG = WellModelGhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpWG = WellModelGhostLastMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                              overlap>;
@@ -221,24 +222,24 @@ using OpWG = WellModelGhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar
 using CommPar = Dune::OwnerOverlapCopyCommunication<int, int>;
 
 template<class Scalar, int Dim>
-using OpFPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
+using OpFPar = Dune::OverlappingSchwarzOperator<BlockSparseMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 CommPar>;
 
 template<class Scalar, int Dim>
-using OpBPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpBPar = Dune::OverlappingSchwarzOperator<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 CommPar>;
 template<class Scalar, int Dim>
-using OpGLFPar = Opm::GhostLastMatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,Dim,Dim>>,
+using OpGLFPar = Opm::GhostLastMatrixAdapter<BlockSparseMatrix<Dune::FieldMatrix<Scalar,Dim,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              CommPar>;
 
 template<class Scalar, int Dim>
-using OpGLBPar = Opm::GhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>,
+using OpGLBPar = Opm::GhostLastMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              CommPar>;
diff --git a/opm/simulators/linalg/PressureBhpTransferPolicy.hpp b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp
index 3d6e78bae92..04b238dabbb 100644
--- a/opm/simulators/linalg/PressureBhpTransferPolicy.hpp
+++ b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp
@@ -23,6 +23,7 @@
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/twolevelmethodcpr.hh>
 
@@ -77,7 +78,7 @@ namespace Opm
 
     namespace Details
     {
-        template<class Scalar> using PressureMatrixType = Dune::BCRSMatrix<MatrixBlock<Scalar, 1, 1>>;
+        template<class Scalar> using PressureMatrixType = BlockSparseMatrix<MatrixBlock<Scalar, 1, 1>>;
         template<class Scalar> using PressureVectorType = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
         template<class Scalar> using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType<Scalar>,
                                                                                  PressureVectorType<Scalar>,
diff --git a/opm/simulators/linalg/PressureTransferPolicy.hpp b/opm/simulators/linalg/PressureTransferPolicy.hpp
index 60d19904a42..dfff9416a83 100644
--- a/opm/simulators/linalg/PressureTransferPolicy.hpp
+++ b/opm/simulators/linalg/PressureTransferPolicy.hpp
@@ -24,13 +24,14 @@
 
 #include <opm/simulators/linalg/twolevelmethodcpr.hh>
 #include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/linalg/WellOperators.hpp>
 
 #include <cstddef>
 
 namespace Opm { namespace Details {
-    template<class Scalar> using PressureMatrixType = Dune::BCRSMatrix<MatrixBlock<Scalar, 1, 1>>;
+    template<class Scalar> using PressureMatrixType = BlockSparseMatrix<MatrixBlock<Scalar, 1, 1>>;
     template<class Scalar> using PressureVectorType = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
     template<class Scalar> using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType<Scalar>,
                                                                              PressureVectorType<Scalar>,
diff --git a/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp b/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp
index aeb116e96c3..889ba6cdb81 100644
--- a/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_gpu_mpi.hpp
@@ -17,7 +17,6 @@
 #ifndef OPM_STANDARDPRECONDITIONERS_GPU_MPI_HEADER
 #define OPM_STANDARDPRECONDITIONERS_GPU_MPI_HEADER
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <string>
 
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
diff --git a/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp
index 49c0f1289e3..99dcf4ed4cb 100644
--- a/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp
@@ -23,8 +23,6 @@
 #include <opm/simulators/linalg/gpuistl/detail/gpu_preconditioner_utils.hpp>
 #endif
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <type_traits>
 
 
diff --git a/opm/simulators/linalg/StandardPreconditioners_serial.hpp b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
index e18c007ef56..2a598b468c5 100644
--- a/opm/simulators/linalg/StandardPreconditioners_serial.hpp
+++ b/opm/simulators/linalg/StandardPreconditioners_serial.hpp
@@ -29,6 +29,8 @@
 #endif
 #endif
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <functional>
 #include <memory>
 #include <type_traits>
@@ -273,8 +275,7 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typen
             const double w = prm.get<double>("relaxation", 1.0);
             using block_type = typename V::block_type;
             using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
-            using matrix_type_to =
-                typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
+            using matrix_type_to = BlockSparseMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
             using GpuILU0 = typename gpuistl::
                 GpuSeqILU0<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
             using Adapter = typename gpuistl::PreconditionerAdapter<VTo, VTo, GpuILU0>;
@@ -337,7 +338,7 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation, typen
 
             using block_type = typename V::block_type;
             using VTo = Dune::BlockVector<Dune::FieldVector<float, block_type::dimension>>;
-            using matrix_type_to = typename Dune::BCRSMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
+            using matrix_type_to = BlockSparseMatrix<Dune::FieldMatrix<float, block_type::dimension, block_type::dimension>>;
             using GpuDILU = typename gpuistl::GpuDILU<matrix_type_to, gpuistl::GpuVector<float>, gpuistl::GpuVector<float>>;
             using MatrixOwner = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<gpuistl::GpuVector<float>,
                 gpuistl::GpuVector<float>, GpuDILU, matrix_type_to>;
diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp
index 44d4a15cc8c..6db661e9368 100644
--- a/opm/simulators/linalg/WellOperators.hpp
+++ b/opm/simulators/linalg/WellOperators.hpp
@@ -24,11 +24,11 @@
 
 #include <dune/common/parallel/communication.hh>
 #include <dune/istl/operators.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <dune/common/shared_ptr.hh>
 #include <dune/istl/paamg/smoother.hh>
 
@@ -57,7 +57,7 @@ class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
 {
 public:
     using field_type = typename X::field_type;
-    using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
+    using PressureMatrix = BlockSparseMatrix<MatrixBlock<field_type, 1, 1>>;
     virtual void addWellPressureEquations(PressureMatrix& jacobian,
                                           const X& weights,
                                           const bool use_well_weights) const = 0;
@@ -228,7 +228,7 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
     using domain_type = X;
     using range_type = Y;
     using field_type = typename X::field_type;
-    using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
+    using PressureMatrix = BlockSparseMatrix<MatrixBlock<field_type, 1, 1>>;
 
     Dune::SolverCategory::Category category() const override
     {
@@ -302,7 +302,7 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X
     using domain_type = X;
     using range_type = Y;
     using field_type = typename X::field_type;
-    using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
+    using PressureMatrix = BlockSparseMatrix<MatrixBlock<field_type, 1, 1>>;
 #if HAVE_MPI
     using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
 #else
diff --git a/opm/simulators/linalg/foreignoverlapfrombcrsmatrix.hh b/opm/simulators/linalg/foreignoverlapfrombcrsmatrix.hh
index 14e0c608c7c..878bea399a7 100644
--- a/opm/simulators/linalg/foreignoverlapfrombcrsmatrix.hh
+++ b/opm/simulators/linalg/foreignoverlapfrombcrsmatrix.hh
@@ -33,7 +33,6 @@
 #include <opm/models/parallel/mpibuffer.hh>
 
 #include <dune/grid/common/datahandleif.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scalarproducts.hh>
 #include <dune/istl/operators.hh>
 
diff --git a/opm/simulators/linalg/globalindices.hh b/opm/simulators/linalg/globalindices.hh
index 61fe0a89aed..73f9be9c8de 100644
--- a/opm/simulators/linalg/globalindices.hh
+++ b/opm/simulators/linalg/globalindices.hh
@@ -28,7 +28,6 @@
 #define EWOMS_GLOBAL_INDICES_HH
 
 #include <dune/grid/common/datahandleif.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scalarproducts.hh>
 #include <dune/istl/operators.hh>
 
diff --git a/opm/simulators/linalg/gpubridge/CprCreation.hpp b/opm/simulators/linalg/gpubridge/CprCreation.hpp
index c3a3bf53cde..7480255777c 100644
--- a/opm/simulators/linalg/gpubridge/CprCreation.hpp
+++ b/opm/simulators/linalg/gpubridge/CprCreation.hpp
@@ -24,6 +24,7 @@
 #include <dune/istl/paamg/matrixhierarchy.hh>
 #include <dune/istl/umfpack.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpubridge/Matrix.hpp>
 #include <opm/simulators/linalg/gpubridge/Preconditioner.hpp>
 
@@ -55,7 +56,7 @@ class CprCreation
 
     BlockedMatrix<Scalar> *mat = nullptr;    // input matrix, blocked
 
-    using DuneMat = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1> >;
+    using DuneMat = BlockSparseMatrix<Dune::FieldMatrix<Scalar, 1, 1> >;
     using DuneVec = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
     using MatrixOperator = Dune::MatrixAdapter<DuneMat, DuneVec, DuneVec>;
     using DuneAmg = Dune::Amg::MatrixHierarchy<MatrixOperator, Dune::Amg::SequentialInformation>;
diff --git a/opm/simulators/linalg/gpubridge/GpuBridge.cpp b/opm/simulators/linalg/gpubridge/GpuBridge.cpp
index 8d49c22213f..2ceaf6ef923 100644
--- a/opm/simulators/linalg/gpubridge/GpuBridge.cpp
+++ b/opm/simulators/linalg/gpubridge/GpuBridge.cpp
@@ -19,8 +19,8 @@
 
 #include <config.h>
 #include <opm/common/TimingMacros.hpp>
-#include "dune/istl/bcrsmatrix.hh"
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/ErrorMacros.hpp>
@@ -359,9 +359,9 @@ initWellContributions([[maybe_unused]] WellContributions<Scalar>& wellContribs,
 
 // the tests use Dune::FieldMatrix, Flow uses Opm::MatrixBlock
 #define INSTANTIATE_GPU_FUNCTIONS(T,n)                                     \
-    template class GpuBridge<Dune::BCRSMatrix<MatrixBlock<T,n,n>>,         \
+    template class GpuBridge<Opm::BlockSparseMatrix<MatrixBlock<T,n,n>>,   \
                    Dune::BlockVector<Dune::FieldVector<T,n>>,n>;           \
-    template class GpuBridge<Dune::BCRSMatrix<Dune::FieldMatrix<T,n,n>>,   \
+    template class GpuBridge<Opm::BlockSparseMatrix<Dune::FieldMatrix<T,n,n>>, \
                              Dune::BlockVector<Dune::FieldVector<T,n>>,n>;
 
 #define INSTANTIATE_TYPE(T)        \
diff --git a/opm/simulators/linalg/gpubridge/GpuBridge.hpp b/opm/simulators/linalg/gpubridge/GpuBridge.hpp
index 8a67534ed55..b46458cc578 100644
--- a/opm/simulators/linalg/gpubridge/GpuBridge.hpp
+++ b/opm/simulators/linalg/gpubridge/GpuBridge.hpp
@@ -70,8 +70,8 @@ class GpuBridge
 
     /// Solve linear system, A*x = b
     /// \warning Values of A might get overwritten!
-    /// \param[in] bridgeMat       matrix A, should be of type Dune::BCRSMatrix
-    /// \param[in] jacMat          matrix A, but modified for the preconditioner, should be of type Dune::BCRSMatrix
+    /// \param[in] bridgeMat       matrix A, should be a matrix like Dune::BCRSMatrix
+    /// \param[in] jacMat          matrix A, but modified for the preconditioner, should a matrix like Dune::BCRSMatrix
     /// \param[in] numJacobiBlocks number of jacobi blocks in jacMat
     /// \param[in] b               vector b, should be of type Dune::BlockVector
     /// \param[in] wellContribs    contains all WellContributions, to apply them separately, instead of adding them to matrix A
diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.cpp b/opm/simulators/linalg/gpuistl/GpuDILU.cpp
index ca19d9607ae..56d544ed0e1 100644
--- a/opm/simulators/linalg/gpuistl/GpuDILU.cpp
+++ b/opm/simulators/linalg/gpuistl/GpuDILU.cpp
@@ -18,11 +18,11 @@
 */
 #include <config.h>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <fmt/core.h>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <opm/simulators/linalg/GraphColoring.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuDILU.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
@@ -535,10 +535,10 @@ GpuDILU<M, X, Y, l>::tuneThreadBlockSizes()
 
 } // namespace Opm::gpuistl
 #define INSTANTIATE_CUDILU_DUNE(realtype, blockdim)                                                                    \
-    template class ::Opm::gpuistl::GpuDILU<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,          \
+    template class ::Opm::gpuistl::GpuDILU<Opm::BlockSparseMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,          \
                                            ::Opm::gpuistl::GpuVector<realtype>,                                        \
                                            ::Opm::gpuistl::GpuVector<realtype>>;                                       \
-    template class ::Opm::gpuistl::GpuDILU<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,           \
+    template class ::Opm::gpuistl::GpuDILU<Opm::BlockSparseMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,           \
                                            ::Opm::gpuistl::GpuVector<realtype>,                                        \
                                            ::Opm::gpuistl::GpuVector<realtype>>
 
diff --git a/opm/simulators/linalg/gpuistl/GpuJac.cpp b/opm/simulators/linalg/gpuistl/GpuJac.cpp
index 741dc162d4e..c3d2cdb7a7a 100644
--- a/opm/simulators/linalg/gpuistl/GpuJac.cpp
+++ b/opm/simulators/linalg/gpuistl/GpuJac.cpp
@@ -21,7 +21,6 @@
 #include <opm/simulators/linalg/gpuistl/GpuJac.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
diff --git a/opm/simulators/linalg/gpuistl/GpuSeqILU0.cpp b/opm/simulators/linalg/gpuistl/GpuSeqILU0.cpp
index 2b63c361f77..be5d153739f 100644
--- a/opm/simulators/linalg/gpuistl/GpuSeqILU0.cpp
+++ b/opm/simulators/linalg/gpuistl/GpuSeqILU0.cpp
@@ -22,10 +22,10 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/is_gpu_operator.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/gpu_constants.hpp>
@@ -324,10 +324,10 @@ GpuSeqILU0<M, X, Y, l>::updateILUConfiguration()
 }
 } // namespace Opm::gpuistl
 #define INSTANTIATE_GPUSEQILU0_DUNE(realtype, blockdim)                                                                  \
-    template class ::Opm::gpuistl::GpuSeqILU0<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,         \
+    template class ::Opm::gpuistl::GpuSeqILU0<Opm::BlockSparseMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,         \
                                             ::Opm::gpuistl::GpuVector<realtype>,                                         \
                                             ::Opm::gpuistl::GpuVector<realtype>>;                                        \
-    template class ::Opm::gpuistl::GpuSeqILU0<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,          \
+    template class ::Opm::gpuistl::GpuSeqILU0<Opm::BlockSparseMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,          \
                                             ::Opm::gpuistl::GpuVector<realtype>,                                         \
                                             ::Opm::gpuistl::GpuVector<realtype>>
 
diff --git a/opm/simulators/linalg/gpuistl/GpuSparseMatrix.cpp b/opm/simulators/linalg/gpuistl/GpuSparseMatrix.cpp
index 6104408b3e2..13d061c2ec7 100644
--- a/opm/simulators/linalg/gpuistl/GpuSparseMatrix.cpp
+++ b/opm/simulators/linalg/gpuistl/GpuSparseMatrix.cpp
@@ -29,7 +29,6 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 #include <cuda.h>
diff --git a/opm/simulators/linalg/gpuistl/GpuSparseMatrixGeneric.cpp b/opm/simulators/linalg/gpuistl/GpuSparseMatrixGeneric.cpp
index ffac016badb..85add2ac45e 100644
--- a/opm/simulators/linalg/gpuistl/GpuSparseMatrixGeneric.cpp
+++ b/opm/simulators/linalg/gpuistl/GpuSparseMatrixGeneric.cpp
@@ -26,7 +26,6 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 #include <cuda.h>
diff --git a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp
index d56b65f3cc7..5c3e645682b 100644
--- a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp
+++ b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp
@@ -18,11 +18,11 @@
 */
 #include <config.h>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <fmt/core.h>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <opm/simulators/linalg/GraphColoring.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
 #include <opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp>
@@ -432,10 +432,10 @@ OpmGpuILU0<M, X, Y, l>::tuneThreadBlockSizes()
 
 } // namespace Opm::gpuistl
 #define INSTANTIATE_GPUILU_DUNE(realtype, blockdim)                                                                    \
-    template class ::Opm::gpuistl::OpmGpuILU0<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,         \
+    template class ::Opm::gpuistl::OpmGpuILU0<Opm::BlockSparseMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>,         \
                                             ::Opm::gpuistl::GpuVector<realtype>,                                         \
                                             ::Opm::gpuistl::GpuVector<realtype>>;                                        \
-    template class ::Opm::gpuistl::OpmGpuILU0<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,          \
+    template class ::Opm::gpuistl::OpmGpuILU0<Opm::BlockSparseMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>,          \
                                             ::Opm::gpuistl::GpuVector<realtype>,                                         \
                                             ::Opm::gpuistl::GpuVector<realtype>>
 
diff --git a/opm/simulators/linalg/gpuistl/PreconditionerConvertFieldTypeAdapter.hpp b/opm/simulators/linalg/gpuistl/PreconditionerConvertFieldTypeAdapter.hpp
index 615b816283d..915fec74ac3 100644
--- a/opm/simulators/linalg/gpuistl/PreconditionerConvertFieldTypeAdapter.hpp
+++ b/opm/simulators/linalg/gpuistl/PreconditionerConvertFieldTypeAdapter.hpp
@@ -19,9 +19,9 @@
 #ifndef OPM_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
 #define OPM_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
 #include <cusparse.h>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/preconditioner.hh>
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/CuMatrixDescription.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/CuSparseHandle.hpp>
@@ -55,10 +55,10 @@ namespace Opm::gpuistl
 //!
 //! using XDouble = Dune::BlockVector<Dune::FieldVector<double, 2>>;
 //! using MDouble = Dune::FieldMatrix<double, 2, 2>;
-//! using SpMatrixDouble = Dune::BCRSMatrix<MDouble>;
+//! using SpMatrixDouble = BlockSparseMatrix<MDouble>;
 //! using XFloat = Dune::BlockVector<Dune::FieldVector<float, 2>>;
 //! using MFloat = Dune::FieldMatrix<float, 2, 2>;
-//! using SpMatrixFloat = Dune::BCRSMatrix<MFloat>;
+//! using SpMatrixFloat = BlockSparseMatrix<MFloat>;
 //!
 //! template<class ParallelInfo>
 //! void applyILU0AsFloat(const MDouble& matrix, const XDouble& x, XDouble& y) {
@@ -107,8 +107,7 @@ class PreconditionerConvertFieldTypeAdapter : public Dune::PreconditionerWithUpd
 
     using XTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
     using YTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
-    using matrix_type_to =
-        typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type_to, block_type::dimension, block_type::dimension>>;
+    using matrix_type_to = BlockSparseMatrix<Dune::FieldMatrix<field_type_to, block_type::dimension, block_type::dimension>>;
 
     //! \brief Constructor.
     //!
diff --git a/opm/simulators/linalg/gpuistl/SolverAdapter.hpp b/opm/simulators/linalg/gpuistl/SolverAdapter.hpp
index 5f33d24e13e..96310a8825d 100644
--- a/opm/simulators/linalg/gpuistl/SolverAdapter.hpp
+++ b/opm/simulators/linalg/gpuistl/SolverAdapter.hpp
@@ -19,7 +19,6 @@
 #ifndef OPM_SOLVERADAPTER_HPP
 #define OPM_SOLVERADAPTER_HPP
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/operators.hh>
 #include <dune/istl/paamg/pinfo.hh>
 #include <dune/istl/schwarz.hh>
diff --git a/opm/simulators/linalg/gpuistl/detail/gpu_preconditioner_utils.hpp b/opm/simulators/linalg/gpuistl/detail/gpu_preconditioner_utils.hpp
index e0749841f0b..7d603ef490b 100644
--- a/opm/simulators/linalg/gpuistl/detail/gpu_preconditioner_utils.hpp
+++ b/opm/simulators/linalg/gpuistl/detail/gpu_preconditioner_utils.hpp
@@ -21,8 +21,7 @@
 #define OPM_GPUISTL_DETAIL_GPU_PRECONDITIONER_UTILS_HEADER
 
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
-
-#include <dune/istl/bcrsmatrix.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 
 namespace Opm::gpuistl::detail {
@@ -44,7 +43,7 @@ namespace Opm::gpuistl::detail {
  * used in the constructor, **not** in the update function or the apply function.
  */
 template<class Operator, class BlockType>
-Dune::BCRSMatrix<BlockType> makeCPUMatrix(const Operator& op) {
+BlockSparseMatrix<BlockType> makeCPUMatrix(const Operator& op) {
     // TODO: Make this more efficient. Maybe we can simply copy the memory areas directly?
     //       Do note that this function is anyway going away when we have a GPU
     //       constructor for the preconditioners, so it is not a priority.
@@ -57,7 +56,7 @@ Dune::BCRSMatrix<BlockType> makeCPUMatrix(const Operator& op) {
     const auto numberOfNonZeroes = gpuMatrix.nonzeroes();
     const auto N = gpuMatrix.N();
 
-    Dune::BCRSMatrix<BlockType> matrix(N, N, numberOfNonZeroes, Dune::BCRSMatrix<BlockType>::row_wise);
+    BlockSparseMatrix<BlockType> matrix(N, N, numberOfNonZeroes, BlockSparseMatrix<BlockType>::row_wise);
     for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) {
         for (auto j = rowIndices[row.index()]; j != rowIndices[row.index() + 1]; ++j) {
             const auto columnIndex = columnIndices[j];
diff --git a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_utilities.hpp b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_utilities.hpp
index 5bb5e6b2474..e12fee8161e 100644
--- a/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_utilities.hpp
+++ b/opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_utilities.hpp
@@ -20,6 +20,7 @@
 #define OPM_GPUISTL_GPUSPARSE_MATRIX_UTILITIES_HPP
 
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/cusparse_safe_call.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/safe_conversion.hpp>
 
@@ -231,13 +232,13 @@ validateVectorMatrixSizes(size_t vectorSize, size_t matrixBlockSize, size_t matr
  */
 #define INSTANTIATE_SPARSE_MATRIX_DUNE_OPERATIONS(CLASS_NAME, realtype, blockdim)                                      \
     template CLASS_NAME<realtype> CLASS_NAME<realtype>::fromMatrix(                                                    \
-        const Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>&, bool);                               \
+        const BlockSparseMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>&, bool);                               \
     template CLASS_NAME<realtype> CLASS_NAME<realtype>::fromMatrix(                                                    \
-        const Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>&, bool);                                \
+        const BlockSparseMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>&, bool);                                \
     template void CLASS_NAME<realtype>::updateNonzeroValues(                                                           \
-        const Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>&, bool);                               \
+        const BlockSparseMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>&, bool);                               \
     template void CLASS_NAME<realtype>::updateNonzeroValues(                                                           \
-        const Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>&, bool)
+        const BlockSparseMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>&, bool)
 
 /**
  * @brief Macro for generating template instantiations for a range of block sizes
diff --git a/opm/simulators/linalg/istlsparsematrixadapter.hh b/opm/simulators/linalg/istlsparsematrixadapter.hh
index cd6fd977d99..244f3fd4183 100644
--- a/opm/simulators/linalg/istlsparsematrixadapter.hh
+++ b/opm/simulators/linalg/istlsparsematrixadapter.hh
@@ -27,10 +27,11 @@
 #ifndef EWOMS_ISTL_SPARSE_MATRIX_ADAPTER_HH
 #define EWOMS_ISTL_SPARSE_MATRIX_ADAPTER_HH
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/version.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 namespace Opm {
 namespace Linear {
 
@@ -43,7 +44,7 @@ class IstlSparseMatrixAdapter
 {
 public:
     //! \brief Implementation of matrix
-    using IstlMatrix = Dune::BCRSMatrix<MatrixBlockType, AllocatorType>;
+    using IstlMatrix = BlockSparseMatrix<MatrixBlockType, AllocatorType>;
 
     //! \brief block type forming the matrix entries
     using MatrixBlock = typename IstlMatrix::block_type;
diff --git a/opm/simulators/linalg/matrixblock.hh b/opm/simulators/linalg/matrixblock.hh
index ee65037bfe9..4c3c522962a 100644
--- a/opm/simulators/linalg/matrixblock.hh
+++ b/opm/simulators/linalg/matrixblock.hh
@@ -289,6 +289,23 @@ public:
         : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
     {}
 };
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+template <typename T, typename I, typename A, int n, int m>
+class UMFPack<IBCRSMatrix<Opm::MatrixBlock<T, n, m>, I, A> >
+    : public UMFPack<IBCRSMatrix<FieldMatrix<T, n, m>, I, A> >
+{
+    using Base = UMFPack<IBCRSMatrix<FieldMatrix<T, n, m>, I, A> >;
+    using Matrix = IBCRSMatrix<FieldMatrix<T, n, m>, I, A>;
+
+public:
+    using RealMatrix = IBCRSMatrix<Opm::MatrixBlock<T, n, m>, I, A>;
+
+    UMFPack(const RealMatrix& matrix, int verbose, bool)
+        : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
+    {}
+};
+#endif
 #endif
 
 #if HAVE_SUPERLU
@@ -309,6 +326,23 @@ public:
         : Base(reinterpret_cast<const Matrix&>(matrix), verb, reuse)
     {}
 };
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+template <typename T, typename I, typename A, int n, int m>
+class SuperLU<IBCRSMatrix<Opm::MatrixBlock<T, n, m>, I, A> >
+    : public SuperLU<IBCRSMatrix<FieldMatrix<T, n, m>, I, A> >
+{
+    using Base = SuperLU<IBCRSMatrix<FieldMatrix<T, n, m>, I, A> >;
+    using Matrix = IBCRSMatrix<FieldMatrix<T, n, m>, I, A>;
+
+public:
+    using RealMatrix = IBCRSMatrix<Opm::MatrixBlock<T, n, m>, I, A>;
+
+    SuperLU(const RealMatrix& matrix, int verb, bool reuse=true)
+        : Base(reinterpret_cast<const Matrix&>(matrix), verb, reuse)
+    {}
+};
+#endif
 #endif
 
 template<typename T, int n, int m>
diff --git a/opm/simulators/linalg/parallelamgbackend.hh b/opm/simulators/linalg/parallelamgbackend.hh
index 457e9731e2c..bff7671cb38 100644
--- a/opm/simulators/linalg/parallelamgbackend.hh
+++ b/opm/simulators/linalg/parallelamgbackend.hh
@@ -33,7 +33,6 @@
 #include "combinedcriterion.hh"
 #include "istlsparsematrixadapter.hh"
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/paamg/amg.hh>
 #include <dune/istl/paamg/pinfo.hh>
 #include <dune/istl/owneroverlapcopy.hh>
diff --git a/opm/simulators/linalg/parallelbasebackend.hh b/opm/simulators/linalg/parallelbasebackend.hh
index 52ede634dcf..e1e6bd0397f 100644
--- a/opm/simulators/linalg/parallelbasebackend.hh
+++ b/opm/simulators/linalg/parallelbasebackend.hh
@@ -43,6 +43,7 @@
 #include <opm/simulators/linalg/linalgparameters.hh>
 #include <opm/simulators/linalg/linalgproperties.hh>
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/overlappingbcrsmatrix.hh>
 #include <opm/simulators/linalg/overlappingblockvector.hh>
 #include <opm/simulators/linalg/overlappingoperator.hh>
@@ -401,7 +402,7 @@ private:
     static constexpr int numEq = getPropValue<TypeTag, Properties::NumEq>();
     using LinearSolverScalar = GetPropType<TypeTag, Properties::LinearSolverScalar>;
     using MatrixBlock = Opm::MatrixBlock<LinearSolverScalar, numEq, numEq>;
-    using NonOverlappingMatrix = Dune::BCRSMatrix<MatrixBlock>;
+    using NonOverlappingMatrix = Opm::BlockSparseMatrix<MatrixBlock>;
 
 public:
     using type = Opm::Linear::OverlappingBCRSMatrix<NonOverlappingMatrix>;
diff --git a/opm/simulators/linalg/superlubackend.hh b/opm/simulators/linalg/superlubackend.hh
index f33b36e3617..14eb91c830b 100644
--- a/opm/simulators/linalg/superlubackend.hh
+++ b/opm/simulators/linalg/superlubackend.hh
@@ -40,6 +40,7 @@
 #include <opm/models/utils/parametersystem.hh>
 
 #include <opm/simulators/linalg/linalgproperties.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 namespace Opm::Properties::TTag {
 
@@ -151,7 +152,7 @@ public:
         using DoubleEqVector = Dune::FieldVector<double, numEq>;
         using DoubleEqMatrix = Dune::FieldMatrix<double, numEq, numEq>;
         using DoubleVector = Dune::BlockVector<DoubleEqVector>;
-        using DoubleMatrix = Dune::BCRSMatrix<DoubleEqMatrix>;
+        using DoubleMatrix = BlockSparseMatrix<DoubleEqMatrix>;
 
         // copy the inputs into the double precision data structures
         DoubleVector bDouble(b);
diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 14f7db21c8e..7e9fb4aedb8 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -25,7 +25,6 @@
 
 #include <dune/common/fmatrix.hh>
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmatrix.hh>
 
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -300,7 +299,7 @@ template<class Scalar> class WellContributions;
 
             const WellInterface<TypeTag>& getWell(const std::string& well_name) const;
 
-            using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
+            using PressureMatrix = BlockSparseMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
 
             void addWellPressureEquations(PressureMatrix& jacobian,
                                           const BVector& weights,
diff --git a/opm/simulators/wells/MSWellHelpers.cpp b/opm/simulators/wells/MSWellHelpers.cpp
index 29c6acb1520..e17e3ebf005 100644
--- a/opm/simulators/wells/MSWellHelpers.cpp
+++ b/opm/simulators/wells/MSWellHelpers.cpp
@@ -33,6 +33,8 @@
 #include <opm/material/densead/Evaluation.hpp>
 #include <opm/material/densead/Math.hpp>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
@@ -360,7 +362,7 @@ ValueType emulsionViscosity(const ValueType& water_fraction,
 template<class Scalar, int Dim>
 using Vec = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
 template<class Scalar, int M, int N = M>
-using Mat = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,M,N>>;
+using Mat = BlockSparseMatrix<Dune::FieldMatrix<Scalar,M,N>>;
 
 #define INSTANTIATE_PARALLELLMSWELLB(T, M, N)                                                                                                                                                                                                                                      \
     template class ParallellMSWellB<Mat<T,M,N>>;                                                                                                                                                                                                 \
diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp
index 729f948e800..e710e1abc95 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.cpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.cpp
@@ -467,7 +467,7 @@ sumDistributed(Parallel::Communication comm)
     template void MultisegmentWellEquations<T,BlackOilDefaultFluidSystemIndices,numWellEq,numEq>::                               \
         extract(Linear::IstlSparseMatrixAdapter<MatrixBlock<T,numEq,numEq>>&) const;           \
     template void MultisegmentWellEquations<T, BlackOilDefaultFluidSystemIndices, numWellEq,numEq>::                               \
-        extractCPRPressureMatrix(Dune::BCRSMatrix<MatrixBlock<T,1,1>>&,                        \
+        extractCPRPressureMatrix(BlockSparseMatrix<MatrixBlock<T,1,1>>&,                        \
                                  const MultisegmentWellEquations<T,BlackOilDefaultFluidSystemIndices,numWellEq,numEq>::BVector&, \
                                  const int,                                                    \
                                  const bool,                                                   \
diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp
index 5e44001de50..2d187314ddb 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.hpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.hpp
@@ -22,12 +22,12 @@
 #ifndef OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
 #define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
 #include <opm/simulators/wells/MSWellHelpers.hpp>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 #include <memory>
@@ -65,11 +65,11 @@ class MultisegmentWellEquations
 
     // the matrix type for the diagonal matrix D
     using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar,numWellEq,numWellEq>;
-    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
+    using DiagMatWell = BlockSparseMatrix<DiagMatrixBlockWellType>;
 
     // the matrix type for the non-diagonal matrix B and C^T
     using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar,numWellEq,numEq>;
-    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
+    using OffDiagMatWell = BlockSparseMatrix<OffDiagMatrixBlockWellType>;
 
     MultisegmentWellEquations(const MultisegmentWellGeneric<Scalar, IndexTraits>& well, const ParallelWellInfo<Scalar>& parallel_well_info);
 
diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp
index 9efa27445f3..0b7e1b9c510 100644
--- a/opm/simulators/wells/StandardWellEquations.cpp
+++ b/opm/simulators/wells/StandardWellEquations.cpp
@@ -422,7 +422,7 @@ sumDistributed(Parallel::Communication comm)
     template void StandardWellEquations<T,BlackOilDefaultFluidSystemIndices,N>::                                        \
         extract(Linear::IstlSparseMatrixAdapter<MatrixBlock<T,N,N>>&) const;          \
     template void StandardWellEquations<T,BlackOilDefaultFluidSystemIndices,N>::                                        \
-        extractCPRPressureMatrix(Dune::BCRSMatrix<MatrixBlock<T,1,1>>&,               \
+        extractCPRPressureMatrix(BlockSparseMatrix<MatrixBlock<T,1,1>>&,               \
                                  const typename StandardWellEquations<T,BlackOilDefaultFluidSystemIndices,N>::BVector&, \
                                  const int,                                           \
                                  const bool,                                          \
diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp
index 2295e07c407..1e814fe1229 100644
--- a/opm/simulators/wells/StandardWellEquations.hpp
+++ b/opm/simulators/wells/StandardWellEquations.hpp
@@ -23,12 +23,12 @@
 #ifndef OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
 #define OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 #include <opm/simulators/wells/WellHelpers.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <dune/common/dynmatrix.hh>
 #include <dune/common/dynvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 namespace Opm
@@ -56,11 +56,11 @@ class StandardWellEquations
 
     // the matrix type for the diagonal matrix D
     using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
-    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
+    using DiagMatWell = BlockSparseMatrix<DiagMatrixBlockWellType>;
 
     // the matrix type for the non-diagonal matrix B and C^T
     using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
-    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
+    using OffDiagMatWell = BlockSparseMatrix<OffDiagMatrixBlockWellType>;
 
     // block vector type
     using BVector = Dune::BlockVector<Dune::FieldVector<Scalar,numEq>>;
diff --git a/opm/simulators/wells/WellHelpers.hpp b/opm/simulators/wells/WellHelpers.hpp
index 4df4e1a4126..77fd8889b3f 100644
--- a/opm/simulators/wells/WellHelpers.hpp
+++ b/opm/simulators/wells/WellHelpers.hpp
@@ -23,9 +23,10 @@
 #ifndef OPM_WELLHELPERS_HEADER_INCLUDED
 #define OPM_WELLHELPERS_HEADER_INCLUDED
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/common/dynmatrix.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <cstdint>
 
 namespace Opm {
@@ -53,7 +54,7 @@ class ParallelStandardWellB
 {
 public:
     using Block = Dune::DynamicMatrix<Scalar>;
-    using Matrix = Dune::BCRSMatrix<Block>;
+    using Matrix = BlockSparseMatrix<Block>;
 
     ParallelStandardWellB(const Matrix& B,
                           const ParallelWellInfo<Scalar>& parallel_well_info);
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index f591ed6b17e..29c21c5f73c 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -41,6 +41,7 @@ namespace Opm {
 
 #include <opm/models/blackoil/blackoilproperties.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/linalgproperties.hh>
 
 #include <opm/simulators/wells/BlackoilWellModel.hpp>
@@ -57,7 +58,6 @@ namespace Opm {
 #include <opm/simulators/utils/DeferredLogger.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmatrix.hh>
 
 #include <opm/material/densead/Evaluation.hpp>
@@ -96,7 +96,7 @@ class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Propertie
     using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
     using Eval = typename Base::Eval;
     using BVector = Dune::BlockVector<VectorBlockType>;
-    using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
+    using PressureMatrix = BlockSparseMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
 
     using WellStateType = WellState<Scalar, IndexTraits>;
     using SingleWellStateType = SingleWellState<Scalar, IndexTraits>;
diff --git a/tests/AmgxPreconditionerTestHelper.hpp b/tests/AmgxPreconditionerTestHelper.hpp
index a549fe8ccb9..6ac4cf72a87 100644
--- a/tests/AmgxPreconditionerTestHelper.hpp
+++ b/tests/AmgxPreconditionerTestHelper.hpp
@@ -24,12 +24,12 @@
 #include <boost/test/unit_test.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/operators.hh>
 #include <dune/istl/solvers.hh>
 
 #include <opm/simulators/linalg/AmgxPreconditioner.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
@@ -93,7 +93,7 @@ createMatrix(int N)
     } else {
         // For GPU matrix, create CPU matrix first then convert
         using ScalarType = typename MatrixType::field_type;
-        using CpuMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<ScalarType, 1, 1>>;
+        using CpuMatrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<ScalarType, 1, 1>>;
         CpuMatrix cpu_matrix;
         setupLaplace2d(N, cpu_matrix);
         return MatrixType::fromMatrix(cpu_matrix);
@@ -214,7 +214,7 @@ testAmgxPreconditioner()
         using Vector = Opm::gpuistl::GpuVector<Scalar>;
         testAmgxPreconditionerGeneric<Matrix, Vector>();
     } else {
-        using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
+        using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
         using Vector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
         testAmgxPreconditionerGeneric<Matrix, Vector>();
     }
@@ -229,7 +229,7 @@ testAmgxPreconditionerUpdate()
         using Vector = Opm::gpuistl::GpuVector<Scalar>;
         testAmgxPreconditionerUpdateGeneric<Matrix, Vector>();
     } else {
-        using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
+        using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
         using Vector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
         testAmgxPreconditionerUpdateGeneric<Matrix, Vector>();
     }
diff --git a/tests/HypreTestHelper.hpp b/tests/HypreTestHelper.hpp
index 52bafa3b0f3..f9d87649f11 100644
--- a/tests/HypreTestHelper.hpp
+++ b/tests/HypreTestHelper.hpp
@@ -24,12 +24,12 @@
 #include <boost/test/unit_test.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/operators.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/HyprePreconditioner.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/hypreinterface/HypreDataStructures.hpp>
@@ -412,7 +412,7 @@ void testMatrixTransfer(bool use_gpu_backend)
     HypreInterface::initialize(use_gpu_backend);
 
     // Create test matrix - start with CPU matrix
-    using CpuMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using CpuMatrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     CpuMatrix cpu_matrix;
     cpu_matrix.setBuildMode(CpuMatrix::row_wise);
     cpu_matrix.setSize(3, 3, 3);
@@ -429,7 +429,7 @@ void testMatrixTransfer(bool use_gpu_backend)
 
     // Convert to target matrix type
     MatrixType matrix = [&]() {
-        if constexpr (std::is_same_v<MatrixType, Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>>) {
+        if constexpr (std::is_same_v<MatrixType, Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>>) {
             return std::move(cpu_matrix);
         }
 #if HAVE_CUDA || HAVE_HIP
@@ -444,7 +444,7 @@ void testMatrixTransfer(bool use_gpu_backend)
     std::vector<HYPRE_BigInt> rows;
     std::vector<HYPRE_BigInt> cols;
 
-    if constexpr (std::is_same_v<MatrixType, Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>>) {
+    if constexpr (std::is_same_v<MatrixType, Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>>) {
         setupMatrixHelperArrays(matrix, ncols, rows, cols);
     }
 #if HAVE_CUDA || HAVE_HIP
diff --git a/tests/gpuistl/test_AmgxInterface.cpp b/tests/gpuistl/test_AmgxInterface.cpp
index 11df7b3d70b..c609de6a51f 100644
--- a/tests/gpuistl/test_AmgxInterface.cpp
+++ b/tests/gpuistl/test_AmgxInterface.cpp
@@ -26,9 +26,9 @@
 #include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <cstddef>
 #include <random>
@@ -127,7 +127,7 @@ BOOST_AUTO_TEST_CASE(TestMatrixTransfer)
     const int blockSize = 1;
 
     using FieldMatrix = Dune::FieldMatrix<double, blockSize, blockSize>;
-    using BCRSMatrix = Dune::BCRSMatrix<FieldMatrix>;
+    using BCRSMatrix = Opm::BlockSparseMatrix<FieldMatrix>;
 
     BCRSMatrix dune_matrix(N, N, BCRSMatrix::row_wise);
 
diff --git a/tests/gpuistl/test_GpuDILU.cpp b/tests/gpuistl/test_GpuDILU.cpp
index 86e25ed6a53..89b344e8e51 100644
--- a/tests/gpuistl/test_GpuDILU.cpp
+++ b/tests/gpuistl/test_GpuDILU.cpp
@@ -22,8 +22,8 @@
 
 #include <boost/test/unit_test.hpp>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <memory>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/DILU.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuDILU.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
@@ -40,8 +40,8 @@ using FM1x1 = Dune::FieldMatrix<T, 1, 1>;
 using FM2x2 = Dune::FieldMatrix<T, 2, 2>;
 using B1x1Vec = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 using B2x2Vec = Dune::BlockVector<Dune::FieldVector<double, 2>>;
-using Sp1x1BlockMatrix = Dune::BCRSMatrix<FM1x1>;
-using Sp2x2BlockMatrix = Dune::BCRSMatrix<FM2x2>;
+using Sp1x1BlockMatrix = Opm::BlockSparseMatrix<FM1x1>;
+using Sp2x2BlockMatrix = Opm::BlockSparseMatrix<FM2x2>;
 using CuMatrix = Opm::gpuistl::GpuSparseMatrixWrapper<T>;
 using CuIntVec = Opm::gpuistl::GpuVector<int>;
 using CuFloatingPointVec = Opm::gpuistl::GpuVector<T>;
diff --git a/tests/gpuistl/test_GpuJac.cpp b/tests/gpuistl/test_GpuJac.cpp
index 2b5a8ad1d1f..8f83ead07f3 100644
--- a/tests/gpuistl/test_GpuJac.cpp
+++ b/tests/gpuistl/test_GpuJac.cpp
@@ -22,7 +22,7 @@
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <cuda_runtime.h>
-#include <dune/istl/bcrsmatrix.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuJac.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
@@ -48,7 +48,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize2, T, NumericTypes)
     constexpr int blocksize = 2;
     const int nonZeroes = 3;
     using M = Dune::FieldMatrix<T, blocksize, blocksize>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     using Vector = Dune::BlockVector<Dune::FieldVector<T, blocksize>>;
     using GpuMatrix = Opm::gpuistl::GpuSparseMatrixWrapper<T>;
     using GpuJac = Opm::gpuistl::GpuJac<GpuMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
@@ -105,7 +105,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize1, T, NumericTypes)
     constexpr int blocksize = 1;
     const int nonZeroes = 8;
     using M = Dune::FieldMatrix<T, blocksize, blocksize>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     using Vector = Dune::BlockVector<Dune::FieldVector<T, blocksize>>;
     using GpuMatrix = Opm::gpuistl::GpuSparseMatrixWrapper<T>;
     using GpuJac = Opm::gpuistl::GpuJac<GpuMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
diff --git a/tests/gpuistl/test_GpuOwnerOverlapCopy.cpp b/tests/gpuistl/test_GpuOwnerOverlapCopy.cpp
index 4e30406d463..b24c8a8ecea 100644
--- a/tests/gpuistl/test_GpuOwnerOverlapCopy.cpp
+++ b/tests/gpuistl/test_GpuOwnerOverlapCopy.cpp
@@ -23,7 +23,6 @@
 
 #include <boost/test/unit_test.hpp>
 #include <dune/common/parallel/mpihelper.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <memory>
 #include <opm/simulators/linalg/gpuistl/GpuOwnerOverlapCopy.hpp>
diff --git a/tests/gpuistl/test_GpuPressureTransferPolicy.cpp b/tests/gpuistl/test_GpuPressureTransferPolicy.cpp
index 471b6778937..dfb6ae070f3 100644
--- a/tests/gpuistl/test_GpuPressureTransferPolicy.cpp
+++ b/tests/gpuistl/test_GpuPressureTransferPolicy.cpp
@@ -25,10 +25,10 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/operators.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PressureTransferPolicy.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
@@ -51,7 +51,7 @@ struct TestFixture {
     static constexpr int N = 10;
 
     using BlockType = Dune::FieldMatrix<Scalar, blockSize, blockSize>;
-    using MatrixType = Dune::BCRSMatrix<BlockType>;
+    using MatrixType = Opm::BlockSparseMatrix<BlockType>;
     using VectorType = Dune::BlockVector<Dune::FieldVector<Scalar, blockSize>>;
     using CommInfo = Dune::Amg::SequentialInformation;
 
diff --git a/tests/gpuistl/test_GpuSeqILU0.cpp b/tests/gpuistl/test_GpuSeqILU0.cpp
index d396157503d..7898728c734 100644
--- a/tests/gpuistl/test_GpuSeqILU0.cpp
+++ b/tests/gpuistl/test_GpuSeqILU0.cpp
@@ -25,8 +25,8 @@
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <dune/common/parallel/mpihelper.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/preconditioners.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSeqILU0.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
 #include <opm/simulators/linalg/gpuistl/PreconditionerAdapter.hpp>
@@ -61,7 +61,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(TestFiniteDifference1D, T, NumericTypes)
     const int N = 5;
     const int nonZeroes = N * 3 - 2;
     using M = Dune::FieldMatrix<T, 1, 1>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     using Vector = Dune::BlockVector<Dune::FieldVector<T, 1>>;
     using GpuILU0 = Opm::gpuistl::GpuSeqILU0<SpMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
 
@@ -156,7 +156,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(TestFiniteDifferenceBlock2, T, NumericTypes)
     const int N = 5;
     const int nonZeroes = N * 3 - 2;
     using M = Dune::FieldMatrix<T, 2, 2>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     using Vector = Dune::BlockVector<Dune::FieldVector<T, 2>>;
     using GpuILU0 = Opm::gpuistl::GpuSeqILU0<SpMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
 
diff --git a/tests/gpuistl/test_GpuSparseMatrix.cu b/tests/gpuistl/test_GpuSparseMatrix.cu
index 936b7556b57..a8fc8d9645a 100644
--- a/tests/gpuistl/test_GpuSparseMatrix.cu
+++ b/tests/gpuistl/test_GpuSparseMatrix.cu
@@ -26,11 +26,10 @@
 #include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
 #include <opm/simulators/linalg/gpuistl/gpu_smart_pointer.hpp>
 #include <opm/simulators/linalg/istlsparsematrixadapter.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <opm/common/utility/pointerArithmetic.hpp>
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <boost/mpl/range_c.hpp>
 #include <boost/test/data/monomorphic.hpp>
 #include <boost/test/data/test_case.hpp>
@@ -63,7 +62,7 @@ BOOST_AUTO_TEST_CASE(TestConstruction1D)
     const int N = 5;
     const int nonZeroes = N * 3 - 2;
     using M = Dune::FieldMatrix<double, 1, 1>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
 
     SpMatrix B(N, N, nonZeroes, SpMatrix::row_wise);
     for (auto row = B.createbegin(); row != B.createend(); ++row) {
@@ -131,7 +130,7 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrs)
     const int N = 3;
     const int nonZeroes = 6;
     using M = Dune::FieldMatrix<double, 1, 1>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
 
     /*
     Fill in a matrix with this sparsity pattern:
@@ -218,7 +217,7 @@ BOOST_AUTO_TEST_CASE(TestGetDiagPtrsBlockedNonScalar)
     const int N = 3;
     const int nonZeroes = 6;
     using M = Dune::FieldMatrix<double, 2, 2>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
 
     /*
     Fill in a blockmatrix with this sparsity pattern:
@@ -319,7 +318,7 @@ void runRandomSparsityMatrixTest()
     std::uniform_real_distribution<double> distribution(0.0, 1.0);
     const std::size_t N = 300;
     using M = Dune::FieldMatrix<double, dim, dim>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
 
     std::vector<std::vector<std::size_t>> nonzerocols(N);
diff --git a/tests/gpuistl/test_GpuSparseTable.cu b/tests/gpuistl/test_GpuSparseTable.cu
index b7b8c2ff8f8..aead89df18f 100644
--- a/tests/gpuistl/test_GpuSparseTable.cu
+++ b/tests/gpuistl/test_GpuSparseTable.cu
@@ -27,8 +27,6 @@
 #include <opm/simulators/linalg/gpuistl/gpu_smart_pointer.hpp>
 #include <opm/simulators/linalg/istlsparsematrixadapter.hh>
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <boost/mpl/range_c.hpp>
 #include <boost/test/data/monomorphic.hpp>
 #include <boost/test/data/test_case.hpp>
diff --git a/tests/gpuistl/test_HypreInterfaceCPU.cpp b/tests/gpuistl/test_HypreInterfaceCPU.cpp
index 9c6279d0d4e..a8444ffd5d6 100644
--- a/tests/gpuistl/test_HypreInterfaceCPU.cpp
+++ b/tests/gpuistl/test_HypreInterfaceCPU.cpp
@@ -30,9 +30,9 @@
 
 #include <opm/simulators/linalg/hypreinterface/HypreInterface.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/gpu_type_detection.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 using namespace Opm::gpuistl;
@@ -45,7 +45,7 @@ BOOST_FIXTURE_TEST_CASE(TestTypeTraits, HypreTestFixture)
 {
     BOOST_CHECK(!is_gpu_type<std::vector<double>>::value);
 
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
     BOOST_CHECK(!is_gpu_type<Matrix>::value);
     BOOST_CHECK(!is_gpu_type<Vector>::value);
@@ -65,7 +65,7 @@ BOOST_FIXTURE_TEST_CASE(TestVectorTransfer_CpuInputCpuBackend, HypreTestFixture)
 
 BOOST_FIXTURE_TEST_CASE(TestMatrixTransfer_CpuInputCpuBackend, HypreTestFixture)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     testMatrixTransfer<Matrix>(false);
 }
 
diff --git a/tests/gpuistl/test_HypreInterfaceGPU.cpp b/tests/gpuistl/test_HypreInterfaceGPU.cpp
index 732f132bdf4..a7e0a4c9f2b 100644
--- a/tests/gpuistl/test_HypreInterfaceGPU.cpp
+++ b/tests/gpuistl/test_HypreInterfaceGPU.cpp
@@ -30,10 +30,10 @@
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
 #include <opm/simulators/linalg/hypreinterface/HypreInterface.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/parallel/mpihelper.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 
@@ -57,7 +57,7 @@ BOOST_FIXTURE_TEST_CASE(TestVectorTransfer_CpuInputGpuBackend, HypreTestFixture)
 
 BOOST_FIXTURE_TEST_CASE(TestMatrixTransfer_CpuInputGpuBackend, HypreTestFixture)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     testMatrixTransfer<Matrix>(true);
 }
 
diff --git a/tests/gpuistl/test_converttofloatadapter.cpp b/tests/gpuistl/test_converttofloatadapter.cpp
index bc0c49e7a43..e138afdd1b7 100644
--- a/tests/gpuistl/test_converttofloatadapter.cpp
+++ b/tests/gpuistl/test_converttofloatadapter.cpp
@@ -25,19 +25,19 @@
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <dune/common/parallel/mpihelper.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/preconditioners.hh>
 #include <limits>
 #include <memory>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/PreconditionerConvertFieldTypeAdapter.hpp>
 
 
 using XDouble = Dune::BlockVector<Dune::FieldVector<double, 2>>;
 using MDouble = Dune::FieldMatrix<double, 2, 2>;
-using SpMatrixDouble = Dune::BCRSMatrix<MDouble>;
+using SpMatrixDouble = Opm::BlockSparseMatrix<MDouble>;
 using XFloat = Dune::BlockVector<Dune::FieldVector<float, 2>>;
 using MFloat = Dune::FieldMatrix<float, 2, 2>;
-using SpMatrixFloat = Dune::BCRSMatrix<MFloat>;
+using SpMatrixFloat = Opm::BlockSparseMatrix<MFloat>;
 namespace
 {
 class TestPreconditioner : Dune::PreconditionerWithUpdate<XFloat, XFloat>
diff --git a/tests/gpuistl/test_cuSparse_matrix_operations.cpp b/tests/gpuistl/test_cuSparse_matrix_operations.cpp
index ea65f328d4d..f5364e6a41d 100644
--- a/tests/gpuistl/test_cuSparse_matrix_operations.cpp
+++ b/tests/gpuistl/test_cuSparse_matrix_operations.cpp
@@ -22,7 +22,7 @@
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <cuda_runtime.h>
-#include <dune/istl/bcrsmatrix.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
 #include <opm/simulators/linalg/gpuistl/PreconditionerAdapter.hpp>
@@ -38,7 +38,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(FlattenAndInvertDiagonalWith3By3Blocks, T, Numeric
     const size_t N = 2;
     const int nonZeroes = 3;
     using M = Dune::FieldMatrix<T, blocksize, blocksize>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     /*
         create this sparse matrix
         | |1 2 3| | 1  0  0| |
@@ -123,7 +123,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(FlattenAndInvertDiagonalWith2By2Blocks, T, Numeric
     const size_t N = 2;
     const int nonZeroes = 3;
     using M = Dune::FieldMatrix<T, blocksize, blocksize>;
-    using SpMatrix = Dune::BCRSMatrix<M>;
+    using SpMatrix = Opm::BlockSparseMatrix<M>;
     /*
         create this sparse matrix
         | |  1 2| | 1  0| |
diff --git a/tests/gpuistl/test_cuVector_operations.cpp b/tests/gpuistl/test_cuVector_operations.cpp
index 17bb33c6f06..7e1f4bcc074 100644
--- a/tests/gpuistl/test_cuVector_operations.cpp
+++ b/tests/gpuistl/test_cuVector_operations.cpp
@@ -22,7 +22,6 @@
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <cuda_runtime.h>
-#include <dune/istl/bcrsmatrix.hh>
 #include <opm/simulators/linalg/gpuistl/GpuJac.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
 #include <opm/simulators/linalg/gpuistl/PreconditionerAdapter.hpp>
diff --git a/tests/gpuistl/test_preconditioner_factory_gpu.cpp b/tests/gpuistl/test_preconditioner_factory_gpu.cpp
index 7e9ff74ab0c..e9572224f4b 100644
--- a/tests/gpuistl/test_preconditioner_factory_gpu.cpp
+++ b/tests/gpuistl/test_preconditioner_factory_gpu.cpp
@@ -27,13 +27,14 @@
 
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory_impl.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 using ScalarT = double;
 constexpr static int Dim = 1;
 
 using CommSeq = Dune::Amg::SequentialInformation;
 
-using MatrixTypeCPU = Dune::BCRSMatrix<Dune::FieldMatrix<ScalarT, Dim, Dim>>;
+using MatrixTypeCPU = Opm::BlockSparseMatrix<Dune::FieldMatrix<ScalarT, Dim, Dim>>;
 using VectorCPU = Dune::BlockVector<Dune::FieldVector<ScalarT, 1>>;
 using CpuOperatorType = Dune::MatrixAdapter<MatrixTypeCPU, VectorCPU, VectorCPU>;
 using FactoryTypeCpu = Opm::PreconditionerFactory<CpuOperatorType, CommSeq>;
diff --git a/tests/gpuistl/test_solver_adapter.cpp b/tests/gpuistl/test_solver_adapter.cpp
index c04d9de743c..0ffb89109a3 100644
--- a/tests/gpuistl/test_solver_adapter.cpp
+++ b/tests/gpuistl/test_solver_adapter.cpp
@@ -22,12 +22,13 @@
 
 #include <boost/test/unit_test.hpp>
 #include <dune/istl/solvers.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/gpuistl/SolverAdapter.hpp>
 
 static const constexpr int dim = 3;
-using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>>;
+using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, dim, dim>>;
 using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
 using Moperator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
 using PrecondFactory = Opm::PreconditionerFactory<Moperator, Dune::Amg::SequentialInformation>;
diff --git a/tests/test_HyprePreconditionerCPU.cpp b/tests/test_HyprePreconditionerCPU.cpp
index 720378ecfe3..e49245bb5e5 100644
--- a/tests/test_HyprePreconditionerCPU.cpp
+++ b/tests/test_HyprePreconditionerCPU.cpp
@@ -31,9 +31,10 @@
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 
 using namespace HypreTestHelpers;
 
@@ -42,7 +43,7 @@ BOOST_GLOBAL_FIXTURE(MPIFixture);
 BOOST_FIXTURE_TEST_CASE(TestHyprePreconditioner_CpuInputCpuBackend, HypreTestFixture)
 {
     constexpr int N = 100; // 100x100 grid
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
     // Create matrix
diff --git a/tests/test_HyprePreconditionerGPU.cpp b/tests/test_HyprePreconditionerGPU.cpp
index cb176990daa..4c561787964 100644
--- a/tests/test_HyprePreconditionerGPU.cpp
+++ b/tests/test_HyprePreconditionerGPU.cpp
@@ -47,7 +47,7 @@ BOOST_GLOBAL_FIXTURE(MPIFixture);
 BOOST_FIXTURE_TEST_CASE(TestHyprePreconditioner_CpuInputGpuBackend, HypreTestFixture)
 {
     constexpr int N = 100; // 100x100 grid
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
     // Create matrix
@@ -62,7 +62,7 @@ BOOST_FIXTURE_TEST_CASE(TestHyprePreconditioner_GpuInputGpuBackend, HypreTestFix
     using namespace Opm::gpuistl;
 
     constexpr int N = 100; // 100x100 grid
-    using CpuMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using CpuMatrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using GpuMatrixType = GpuSparseMatrixWrapper<double>;
     using GpuVectorType = GpuVector<double>;
 
@@ -81,7 +81,7 @@ BOOST_FIXTURE_TEST_CASE(TestHyprePreconditioner_GpuInputCpuBackend, HypreTestFix
     using namespace Opm::gpuistl;
 
     constexpr int N = 100; // 100x100 grid
-    using CpuMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using CpuMatrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using GpuMatrixType = GpuSparseMatrixWrapper<double>;
     using GpuVectorType = GpuVector<double>;
 
diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp
index 39f906c4c71..4d365204de3 100644
--- a/tests/test_blackoil_amg.cpp
+++ b/tests/test_blackoil_amg.cpp
@@ -24,6 +24,7 @@
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 
 #include <dune/common/parallel/mpihelper.hh>
@@ -32,7 +33,6 @@
 #include <dune/common/parallel/indexset.hh>
 #include <dune/common/parallel/plocalindex.hh>
 #include <dune/common/parallel/communication.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/schwarz.hh>
 
@@ -270,7 +270,7 @@ void runBlackoilAmgLaplace()
 {
     constexpr int BS=2, N=100;
     typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
-    typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+    typedef Opm::BlockSparseMatrix<MatrixBlock> BCRSMat;
     typedef Dune::FieldVector<double,BS> VectorBlock;
     typedef Dune::BlockVector<VectorBlock> Vector;
     typedef int GlobalId;
diff --git a/tests/test_csrToCscOffsetMap.cpp b/tests/test_csrToCscOffsetMap.cpp
index 34287b91e63..65b74da22ed 100644
--- a/tests/test_csrToCscOffsetMap.cpp
+++ b/tests/test_csrToCscOffsetMap.cpp
@@ -6,14 +6,14 @@
 #include <boost/version.hpp>
 
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpubridge/opencl/openclBISAI.hpp>
 
 BOOST_AUTO_TEST_CASE(testcsrtocscoffsetmap){
 
-    using Matrix = Dune::BCRSMatrix<double>;
+    using Matrix = Opm::BlockSparseMatrix<double>;
 
     Matrix matrix;
     {
diff --git a/tests/test_cusparseSolver.cpp b/tests/test_cusparseSolver.cpp
index 79917f32be3..de15b54f85b 100644
--- a/tests/test_cusparseSolver.cpp
+++ b/tests/test_cusparseSolver.cpp
@@ -25,10 +25,10 @@
 
 #include <opm/simulators/linalg/gpubridge/GpuBridge.hpp>
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
@@ -45,7 +45,7 @@ class DeviceInitException : public std::logic_error
 };
 
 template <int bz>
-using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 template <int bz>
 using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
diff --git a/tests/test_dilu.cpp b/tests/test_dilu.cpp
index 22f572adc6c..095047a87b9 100644
--- a/tests/test_dilu.cpp
+++ b/tests/test_dilu.cpp
@@ -16,11 +16,11 @@
 
 #include <config.h>
 #include <opm/simulators/linalg/DILU.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <boost/mpl/list.hpp>
 #include <boost/test/unit_test.hpp>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/preconditioners.hh>
 
 using NumericTypes = boost::mpl::list<double, float>;
@@ -43,7 +43,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes)
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 4;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -111,7 +111,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2, T, NumericTypes)
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 3;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -179,7 +179,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes)
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 4;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -288,7 +288,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect1, T, NumericTypes)
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 3;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -391,7 +391,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect2, T, NumericTypes)
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 3;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -497,7 +497,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect3x3, T, NumericTypes)
     const int N = 3;
     constexpr int bz = 3;
     const int nonZeroes = 5;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -620,7 +620,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes)
     const int N = 3;
     constexpr int bz = 3;
     const int nonZeroes = 5;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
@@ -795,7 +795,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsEqualToDuneSeqILUApply, T, NumericTy
     const int N = 2;
     constexpr int bz = 2;
     const int nonZeroes = 2;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
     Matrix A(N, N, nonZeroes, Matrix::row_wise);
diff --git a/tests/test_extractMatrix.cpp b/tests/test_extractMatrix.cpp
index 29ee8fe408d..b1d905266bc 100644
--- a/tests/test_extractMatrix.cpp
+++ b/tests/test_extractMatrix.cpp
@@ -23,12 +23,12 @@
 #include <boost/test/unit_test.hpp>
 
 #include <opm/simulators/linalg/extractMatrix.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 using B = Dune::FieldMatrix<double, 2, 2>;
-using M = Dune::BCRSMatrix<B>;
+using M = Opm::BlockSparseMatrix<B>;
 using V = Dune::BlockVector<Dune::FieldVector<double, 2>>;
 
 M build3x3BlockMatrix()
diff --git a/tests/test_flexiblesolver.cpp b/tests/test_flexiblesolver.cpp
index ad6ef0cf05f..8101b10ca42 100644
--- a/tests/test_flexiblesolver.cpp
+++ b/tests/test_flexiblesolver.cpp
@@ -24,11 +24,11 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 
 #include <fstream>
@@ -39,7 +39,7 @@ template <int bz>
 Dune::BlockVector<Dune::FieldVector<double, bz>>
 testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
 {
-    using Matrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Opm::MatrixBlock<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
     Matrix matrix;
     {
@@ -47,7 +47,7 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        using M = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
         readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
     }
     Vector rhs;
diff --git a/tests/test_ghostlastmatrixadapter.cpp b/tests/test_ghostlastmatrixadapter.cpp
index 63ac6897205..529afdc13b5 100644
--- a/tests/test_ghostlastmatrixadapter.cpp
+++ b/tests/test_ghostlastmatrixadapter.cpp
@@ -23,6 +23,7 @@
 #include <boost/test/unit_test.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 #include <opm/simulators/linalg/WellOperators.hpp>
 
@@ -32,7 +33,6 @@
 #include <dune/common/parallel/indexset.hh>
 #include <dune/common/parallel/plocalindex.hh>
 #include <dune/common/parallel/communication.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/schwarz.hh>
 
@@ -265,7 +265,7 @@ int main(int argc, char** argv)
 
     constexpr int BS=2, N=100;
     using MatrixBlock = Dune::FieldMatrix<double,BS,BS>;
-    using BCRSMat = Dune::BCRSMatrix<MatrixBlock>;
+    using BCRSMat = Opm::BlockSparseMatrix<MatrixBlock>;
     using VectorBlock = Dune::FieldVector<double,BS>;
     using Vector = Dune::BlockVector<VectorBlock>;
     using Communication = Dune::OwnerOverlapCopyCommunication<int>;
diff --git a/tests/test_graphcoloring.cpp b/tests/test_graphcoloring.cpp
index 335b49d0d22..15e1915c078 100644
--- a/tests/test_graphcoloring.cpp
+++ b/tests/test_graphcoloring.cpp
@@ -1,9 +1,9 @@
 #include <config.h>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/paamg/graph.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/GraphColoring.hpp>
 
 #define BOOST_TEST_MODULE GraphColoringTest
@@ -26,7 +26,7 @@ void checkAllIndices(const std::vector<std::size_t>& ordering)
 
 BOOST_AUTO_TEST_CASE(TestWelschPowell)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+    using Matrix = Opm::BlockSparseMatrix<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);
@@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms3x3Matrix)
     const int nonZeroes = 5;
 
     // creating some shorthand typenames
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 
     Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise);
     for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) {
@@ -170,7 +170,7 @@ BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Simple)
     const int nonZeroes = 11;
 
     // creating some shorthand typenames
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 
     Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise);
     for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) {
@@ -250,7 +250,7 @@ BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Tridiagonal)
     const int nonZeroes = 13;
 
     // creating some shorthand typenames
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 
     Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise);
     for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) {
@@ -325,7 +325,7 @@ BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Complex)
     const int nonZeroes = 15;
 
     // creating some shorthand typenames
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 
     Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise);
     for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) {
diff --git a/tests/test_milu.cpp b/tests/test_milu.cpp
index 605076d5858..2db12912bbe 100644
--- a/tests/test_milu.cpp
+++ b/tests/test_milu.cpp
@@ -5,11 +5,11 @@
 #include<vector>
 #include<memory>
 
-#include<dune/istl/bcrsmatrix.hh>
 #include<dune/istl/bvector.hh>
 #include<dune/common/version.hh>
 #include<dune/common/fmatrix.hh>
 #include<dune/common/fvector.hh>
+#include<opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include<opm/simulators/linalg/ParallelOverlappingILU0.hpp>
 
 #include <opm/common/ErrorMacros.hpp>
@@ -104,13 +104,13 @@ void test_milu0(M& A)
 }
 
 template<class B, class Alloc>
-void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
+void setupSparsityPattern(Opm::BlockSparseMatrix<B,Alloc>& A, int N)
 {
-  typedef typename Dune::BCRSMatrix<B,Alloc> Matrix;
-  A.setSize(N*N, N*N, N*N*5);
+    typedef typename Opm::BlockSparseMatrix<B,Alloc> Matrix;
   A.setBuildMode(Matrix::row_wise);
+  A.setSize(N*N, N*N, N*N*5);
 
-  for (typename Dune::BCRSMatrix<B,Alloc>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
+    for (typename Opm::BlockSparseMatrix<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
 
@@ -135,9 +135,9 @@ void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
 
 
 template<class B, class Alloc>
-void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
+void setupLaplacian(Opm::BlockSparseMatrix<B,Alloc>& A, int N)
 {
-  typedef typename Dune::BCRSMatrix<B,Alloc>::field_type FieldType;
+    typedef typename Opm::BlockSparseMatrix<B,Alloc>::field_type FieldType;
 
   setupSparsityPattern(A,N);
 
@@ -150,7 +150,7 @@ void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
     b->operator[](b.index())=-1.0;
 
 
-  for (typename Dune::BCRSMatrix<B,Alloc>::RowIterator i = A.begin(); i != A.end(); ++i) {
+    for (typename Opm::BlockSparseMatrix<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
 
@@ -174,7 +174,7 @@ template<int bsize>
 void test()
 {
     std::size_t N = 32;
-    Dune::BCRSMatrix<Dune::FieldMatrix<double, bsize, bsize> > A;
+    Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bsize, bsize> > A;
     setupLaplacian(A, N);
     test_milu0(A);
 #ifdef DEBUG
diff --git a/tests/test_openclSolver.cpp b/tests/test_openclSolver.cpp
index 8dada17f118..478a79de3f6 100644
--- a/tests/test_openclSolver.cpp
+++ b/tests/test_openclSolver.cpp
@@ -25,10 +25,10 @@
 
 #include <opm/simulators/linalg/gpubridge/GpuBridge.hpp>
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
@@ -43,7 +43,7 @@ class PlatformInitException : public std::logic_error
 };
 
 template <int bz>
-using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 template <int bz>
 using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
diff --git a/tests/test_preconditionerfactory.cpp b/tests/test_preconditionerfactory.cpp
index 66b38008788..21a627b664a 100644
--- a/tests/test_preconditionerfactory.cpp
+++ b/tests/test_preconditionerfactory.cpp
@@ -28,11 +28,11 @@
 #include <opm/simulators/linalg/PreconditionerFactory_impl.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/solvers.hh>
 
@@ -68,7 +68,7 @@ template <int bz>
 Dune::BlockVector<Dune::FieldVector<double, bz>>
 testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
 {
-    using Matrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
+    using Matrix = Opm::BlockSparseMatrix<Opm::MatrixBlock<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
     Matrix matrix;
     {
@@ -76,7 +76,7 @@ testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        using M = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
         readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
     }
     Vector rhs;
@@ -171,7 +171,7 @@ BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
 
 
 template <int bz>
-using M = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
+using M = Opm::BlockSparseMatrix<Opm::MatrixBlock<double, bz, bz>>;
 template <int bz>
 using V = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 template <int bz>
@@ -292,7 +292,7 @@ testPrecRepeating(const Opm::PropertyTree& prm, const std::string& matrix_filena
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        using M = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wstrict-aliasing"
         readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
diff --git a/tests/test_rocalutionSolver.cpp b/tests/test_rocalutionSolver.cpp
index 3da8146f48b..665eb29398f 100644
--- a/tests/test_rocalutionSolver.cpp
+++ b/tests/test_rocalutionSolver.cpp
@@ -25,11 +25,11 @@
 
 #include <opm/simulators/linalg/gpubridge/GpuBridge.hpp>
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <rocalution/rocalution.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
@@ -38,7 +38,7 @@
 #include <boost/property_tree/ptree.hpp>
 
 template <int bz>
-using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 template <int bz>
 using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
diff --git a/tests/test_rocsparseSolver.cpp b/tests/test_rocsparseSolver.cpp
index 56d78659d9f..db83db4438d 100644
--- a/tests/test_rocsparseSolver.cpp
+++ b/tests/test_rocsparseSolver.cpp
@@ -29,10 +29,10 @@
 
 #include <opm/simulators/linalg/gpubridge/GpuBridge.hpp>
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
@@ -47,7 +47,7 @@ class HIPInitException : public std::logic_error
 };
 
 template <int bz>
-using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+using Matrix = Opm::BlockSparseMatrix<Dune::FieldMatrix<double, bz, bz>>;
 template <int bz>
 using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 
diff --git a/tests/test_solvetransposed3x3.cpp b/tests/test_solvetransposed3x3.cpp
index 3be615e8681..b8c4b8f85ef 100644
--- a/tests/test_solvetransposed3x3.cpp
+++ b/tests/test_solvetransposed3x3.cpp
@@ -29,8 +29,6 @@
 #include <boost/test/tools/floating_point_comparison.hpp>
 #endif
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <opm/simulators/linalg/gpubridge/opencl/openclCPR.hpp>
 #include <opm/simulators/linalg/gpubridge/Misc.hpp>