diff --git a/CMakeLists.txt b/CMakeLists.txt
index 52867ea6..e8a3c37c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -110,6 +110,12 @@ if(NOT PETSC_FOUND)
                    ${PROJECT_SOURCE_DIR}/src/LinAlg/PETScSolParams.C)
 endif()
 
+if(NOT ISTL_FOUND)
+  list(REMOVE_ITEM IFEM_SRCS
+                   ${PROJECT_SOURCE_DIR}/src/LinAlg/ISTLMatrix.C
+                   ${PROJECT_SOURCE_DIR}/src/LinAlg/ISTLSolParams.C)
+endif()
+
 if(LRSPLINE_FOUND OR LRSpline_FOUND)
   file(GLOB LR_SRCS ${PROJECT_SOURCE_DIR}/src/ASM/LR/*.C)
   list(APPEND IFEM_SRCS ${LR_SRCS})
diff --git a/cmake/Modules/FindIFEMDeps.cmake b/cmake/Modules/FindIFEMDeps.cmake
index 2cc2b1c9..e167d691 100644
--- a/cmake/Modules/FindIFEMDeps.cmake
+++ b/cmake/Modules/FindIFEMDeps.cmake
@@ -211,6 +211,16 @@ IF(IFEM_USE_OPENMP)
   ENDIF(OPENMP_FOUND)
 ENDIF(IFEM_USE_OPENMP)
 
+# ISTL
+if(IFEM_USE_ISTL)
+  find_package(ISTL)
+  if(ISTL_FOUND)
+    list(APPEND IFEM_DEPLIBS ${ISTL_LIBRARIES})
+    list(APPEND IFEM_DEPINCLUDES ${ISTL_INCLUDE_DIRS})
+    list(APPEND IFEM_DEFINITIONS -DHAS_ISTL=1 ${ISTL_DEFINITIONS})
+  endif()
+endif()
+
 # Portability issues
 include(CheckFunctionExists)
 set(CMAKE_REQUIRED_DEFINITIONS)
diff --git a/cmake/Modules/FindISTL.cmake b/cmake/Modules/FindISTL.cmake
new file mode 100644
index 00000000..2aee35aa
--- /dev/null
+++ b/cmake/Modules/FindISTL.cmake
@@ -0,0 +1,38 @@
+find_package(PkgConfig)
+
+set(OLD_PKG $ENV{PKG_CONFIG_PATH})
+set(ENV{PKG_CONFIG_PATH} $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig)
+pkg_check_modules(ISTL dune-istl)
+set(ENV{PKG_CONFIG_PATH} ${OLD_PKG})
+
+list(APPEND ISTL_DEFINITIONS -DHAVE_NULLPTR=1 -DISTL_VERSION="${ISTL_VERSION}")
+string(REPLACE "." ";" VERSION_LIST "${ISTL_VERSION}")
+if(VERSION_LIST)
+  list(GET VERSION_LIST 0 DUNE_ISTL_VERSION_MAJOR)
+  list(GET VERSION_LIST 1 DUNE_ISTL_VERSION_MINOR)
+  list(GET VERSION_LIST 2 DUNE_ISTL_VERSION_PATCH)
+  list(APPEND ISTL_DEFINITIONS -DDUNE_ISTL_VERSION_MAJOR=${DUNE_ISTL_VERSION_MAJOR}
+                               -DDUNE_ISTL_VERSION_MINOR=${DUNE_ISTL_VERSION_MINOR}
+                               -DDUNE_ISTL_VERSION_REVISION=${DUNE_ISTL_VERSION_PATCH})
+endif()
+
+if(SUPERLU_FOUND)
+  list(APPEND ISTL_DEFINITIONS -DHAVE_SUPERLU=1 -DSUPERLU_POST_2005_VERSION=1 -DSUPERLU_NTYPE=1 -DSUPERLU_MIN_VERSION_4_3=1)
+endif()
+
+if(IFEM_USE_UMFPACK)
+  find_package(SuiteSparse COMPONENTS umfpack)
+endif()
+
+if(SuiteSparse_UMFPACK_FOUND)
+  list(APPEND ISTL_DEFINITIONS -DHAVE_UMFPACK=1)
+  list(APPEND ISTL_INCLUDE_DIRS ${UMFPACK_INCLUDE_DIR})
+  list(APPEND ISTL_LIBRARIES ${UMFPACK_LIBRARY})
+endif()
+
+
+list(APPEND ISTL_INCLUDE_DIRS ${ISTL_INCLUDEDIR})
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(ISTL DEFAULT_MSG
+                                  ISTL_INCLUDE_DIRS ISTL_LIBRARIES)
diff --git a/cmake/Modules/FindSuiteSparse.cmake b/cmake/Modules/FindSuiteSparse.cmake
new file mode 100644
index 00000000..0b77b836
--- /dev/null
+++ b/cmake/Modules/FindSuiteSparse.cmake
@@ -0,0 +1,294 @@
+# - Find Tim Davis' SuiteSparse collection of sparse matrix libraries
+#
+# Synopsis:
+#   find_package (SuiteSparse COMPONENTS <list-of-components>)
+#
+# Components are:
+#   amd              Approximate Minimum Degree ordering
+#   camd             Constrained Approximate Minimum Degree ordering
+#   colamd           COLumn Approximate Minimum Degree ordering
+#   ccolamd          Constrained COLumn Approximate Minimum Degree ordering
+#   cholmod          Supernodal sparse Cholesky factorization and update
+#   umfpack          Unsymmetric MultiFrontal sparse LU factorization
+#
+# The following variables will be set:
+#
+#   SuiteSparse_FOUND            True if all dependencies are satisfied
+#   SuiteSparse_Xxx_FOUND        True if module Xxx is found
+#   HAVE_SUITESPARSE_Xxx_H       Binary value indicating presence of header
+#   SuiteSparse_INCLUDE_DIRS     Paths containing the SuiteSparse header files
+#   SuiteSparse_LIBRARIES        Name of the libraries which must be linked
+#   SuiteSparse_DEFINITIONS      Defines that must be passed to the compiler
+#   SuiteSparse_LINKER_FLAGS     Options that must be passed when linking
+#
+# The following options can be set to configure the module:
+#
+#   SUITESPARSE_USE_STATIC       Link with a static library, even if a
+#                                dynamic library is also present. Note that
+#                                setting this to OFF does not ensure that a
+#                                shared library will be used.
+#
+# See <http://www.cise.ufl.edu/research/sparse/SuiteSparse>.
+
+# Copyright (C) 2012 Uni Research AS
+# This file is licensed under the GNU General Public License v3.0
+
+function (try_compile_umfpack varname)
+  include (CMakePushCheckState)
+  include (CheckCSourceCompiles)
+  cmake_push_check_state ()
+  set (CMAKE_REQUIRED_INCLUDES ${UMFPACK_INCLUDE_DIRS})
+  set (CMAKE_REQUIRED_LIBRARIES ${UMFPACK_LIBRARY} ${ARGN} ${SuiteSparse_EXTRA_LIBS})
+  check_c_source_compiles (
+        "#include <umfpack.h>
+int main (void) {
+  void *Symbolic, *Numeric;
+  double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL];
+  umfpack_dl_defaults(Control);
+  umfpack_dl_symbolic(0, 0, 0, 0, 0,
+                      &Symbolic, Control, Info);
+  umfpack_dl_numeric (0, 0, 0,
+                      Symbolic, &Numeric, Control, Info);
+  umfpack_dl_free_symbolic(&Symbolic);
+  umfpack_dl_solve(UMFPACK_A, 0, 0, 0, 0, 0,
+                   Numeric, Control, Info);
+  umfpack_dl_free_numeric(&Numeric);
+  umfpack_timer ();
+  return 0;
+}" ${varname})
+  cmake_pop_check_state ()
+  set (${varname} "${${varname}}" PARENT_SCOPE)
+endfunction (try_compile_umfpack varname)
+
+# variables to pass on to other packages
+if (FIND_QUIETLY)
+  set (SuiteSparse_QUIET "QUIET")
+else (FIND_QUIETLY)
+  set (SuiteSparse_QUIET "")
+endif (FIND_QUIETLY)
+
+# we need to link to BLAS and LAPACK
+if (NOT BLAS_FOUND)
+  find_package (BLAS ${SuiteSparse_QUIET} REQUIRED)
+endif (NOT BLAS_FOUND)
+if (NOT LAPACK_FOUND)
+  find_package (LAPACK ${SuiteSparse_QUIET} REQUIRED)
+endif (NOT LAPACK_FOUND)
+
+# we also need the math part of the runtime library
+find_library (MATH_LIBRARY NAMES "m")
+set (SuiteSparse_EXTRA_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MATH_LIBRARY})
+
+# if we don't get any further clues about where to look, then start
+# roaming around the system
+set (_no_default_path "")
+
+# search system directories by default
+set (SuiteSparse_SEARCH_PATH)
+
+# pick up paths from the environment if specified there; these replace the
+# pre-defined paths so that we don't accidentially pick up old stuff
+if (NOT $ENV{SuiteSparse_DIR} STREQUAL "")
+  set (SuiteSparse_SEARCH_PATH "$ENV{SuiteSparse_DIR}")
+endif (NOT $ENV{SuiteSparse_DIR} STREQUAL "")
+if (SuiteSparse_DIR)
+  set (SuiteSparse_SEARCH_PATH "${SuiteSparse_DIR}")
+endif (SuiteSparse_DIR)
+# CMake uses _DIR suffix as default for config-mode files; it is unlikely
+# that we are building SuiteSparse ourselves; use _ROOT suffix to specify
+# location to pre-canned binaries
+if (NOT $ENV{SuiteSparse_ROOT} STREQUAL "")
+  set (SuiteSparse_SEARCH_PATH "$ENV{SuiteSparse_ROOT}")
+endif (NOT $ENV{SuiteSparse_ROOT} STREQUAL "")
+if (SuiteSparse_ROOT)
+  set (SuiteSparse_SEARCH_PATH "${SuiteSparse_ROOT}")
+endif (SuiteSparse_ROOT)
+# most commonly, we use the uppercase version of this variable
+if (SUITESPARSE_ROOT)
+  set (SuiteSparse_SEARCH_PATH "${SUITESPARSE_ROOT}")
+endif (SUITESPARSE_ROOT)
+
+# if we have specified a search path, then confine ourselves to that
+if (SuiteSparse_SEARCH_PATH)
+  set (_no_default_path "NO_DEFAULT_PATH")
+endif (SuiteSparse_SEARCH_PATH)
+
+# transitive closure of dependencies; after this SuiteSparse_MODULES is the
+# full list of modules that must be found to satisfy the user's link demands
+set (SuiteSparse_MODULES ${SuiteSparse_FIND_COMPONENTS})
+list (FIND SuiteSparse_MODULES "umfpack" UMFPACK_DESIRED)
+if (NOT UMFPACK_DESIRED EQUAL -1)
+  list (APPEND SuiteSparse_MODULES amd cholmod)
+endif (NOT UMFPACK_DESIRED EQUAL -1)
+list (FIND SuiteSparse_MODULES "cholmod" CHOLMOD_DESIRED)
+if (NOT CHOLMOD_DESIRED EQUAL -1)
+  list (APPEND SuiteSparse_MODULES amd camd colamd)
+endif (NOT CHOLMOD_DESIRED EQUAL -1)
+if (SuiteSparse_MODULES)
+  list (REMOVE_DUPLICATES SuiteSparse_MODULES)
+endif (SuiteSparse_MODULES)
+
+# if someone else already have found all the packages for us, then don't do anything
+set (SuiteSparse_EVERYTHING_FOUND TRUE)
+foreach (module IN LISTS SuiteSparse_MODULES)
+  string (TOUPPER ${module} MODULE)
+  if (NOT SuiteSparse_${MODULE}_FOUND)
+	set (SuiteSparse_EVERYTHING_FOUND FALSE)
+	break ()
+  endif (NOT SuiteSparse_${MODULE}_FOUND)
+endforeach (module)
+if (SuiteSparse_EVERYTHING_FOUND)
+  return ()
+endif (SuiteSparse_EVERYTHING_FOUND)
+
+# only search in architecture-relevant directory
+if (CMAKE_SIZEOF_VOID_P)
+  math (EXPR _BITS "8 * ${CMAKE_SIZEOF_VOID_P}")
+endif (CMAKE_SIZEOF_VOID_P)
+
+# if we are told to link SuiteSparse statically, add these parts
+# to the name so we always match only that particular type of lib
+option (SUITESPARSE_USE_STATIC "Link SuiteSparse statically" OFF)
+mark_as_advanced (SUITESPARSE_USE_STATIC)
+if (SUITESPARSE_USE_STATIC)
+  set (_pref_ "${CMAKE_STATIC_LIBRARY_PREFIX}")
+  set (_suff_ "${CMAKE_STATIC_LIBRARY_SUFFIX}")
+else (SUITESPARSE_USE_STATIC)
+  set (_pref_ "")
+  set (_suff_ "")
+endif (SUITESPARSE_USE_STATIC)
+
+# if SuiteSparse >= 4.0 we must also link with libsuitesparseconfig
+# assume that this is the case if we find the library; otherwise just
+# ignore it (older versions don't have a file named like this)
+find_library (config_LIBRARY
+  NAMES "${_pref_}suitesparseconfig${_suff_}"
+  PATHS ${SuiteSparse_SEARCH_PATH}
+  PATH_SUFFIXES ".libs" "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}" "lib/ufsparse"
+  ${_no_default_path}
+  )
+if (config_LIBRARY)
+  list (APPEND SuiteSparse_EXTRA_LIBS ${config_LIBRARY})
+  # POSIX.1-2001 REALTIME portion require us to link this library too for
+  # clock_gettime() which is used by suitesparseconfig
+  if ("${CMAKE_SYSTEM_NAME}" MATCHES "Linux")
+	list (APPEND SuiteSparse_EXTRA_LIBS "-lrt")
+  endif ("${CMAKE_SYSTEM_NAME}" MATCHES "Linux")
+endif (config_LIBRARY)
+
+# search filesystem for each of the module individually
+foreach (module IN LISTS SuiteSparse_MODULES)
+  string (TOUPPER ${module} MODULE)
+  # search for files which implements this module
+  find_path (${MODULE}_INCLUDE_DIR
+	NAMES ${module}.h
+	PATHS ${SuiteSparse_SEARCH_PATH}
+	PATH_SUFFIXES "include" "include/suitesparse" "include/ufsparse" "${MODULE}/Include"
+	${_no_default_path}
+	)
+
+  find_library (${MODULE}_LIBRARY
+	NAMES "${_pref_}${module}${_suff_}"
+	PATHS ${SuiteSparse_SEARCH_PATH}
+	PATH_SUFFIXES "lib/.libs" "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}" "lib/ufsparse" "${MODULE}/Lib"
+	${_no_default_path}
+	)
+  # start out by including the module itself; other dependencies will be added later
+  set (${MODULE}_INCLUDE_DIRS ${${MODULE}_INCLUDE_DIR})
+  set (${MODULE}_LIBRARIES ${${MODULE}_LIBRARY})
+endforeach (module)
+
+# insert any inter-modular dependencies here
+if (CHOLMOD_LIBRARY)
+  list (APPEND CHOLMOD_LIBRARIES ${AMD_LIBRARIES} ${COLAMD_LIBRARIES})
+  # optional libraries; don't insert any -NOT_FOUND paths
+  if (CAMD_LIBRARY)
+	list (APPEND CHOLMOD_LIBRARIES ${CAMD_LIBRARIES})
+  endif (CAMD_LIBRARY)
+  if (CCOLAMD_LIBRARY)
+	list (APPEND CHOLMOD_LIBRARIES ${CCOLAMD_LIBRARIES})
+  endif (CCOLAMD_LIBRARY)
+  list (REVERSE CHOLMOD_LIBRARIES)
+  # always remove the *first* library from the list
+  list (REMOVE_DUPLICATES CHOLMOD_LIBRARIES)
+  list (REVERSE CHOLMOD_LIBRARIES)
+endif (CHOLMOD_LIBRARY)
+if (UMFPACK_LIBRARY)
+  set (UMFPACK_EXTRA_LIBS)
+  # test if umfpack is usable with only amd and not cholmod
+  try_compile_umfpack (HAVE_UMFPACK_WITHOUT_CHOLMOD ${AMD_LIBRARIES})
+  if (HAVE_UMFPACK_WITHOUT_CHOLMOD)
+	list (APPEND UMFPACK_EXTRA_LIBS ${AMD_LIBRARIES})
+  else (HAVE_UMFPACK_WITHOUT_CHOLMOD)
+	if (CHOLMOD_LIBRARIES)
+	  try_compile_umfpack (HAVE_UMFPACK_WITH_CHOLMOD ${CHOLMOD_LIBRARIES})
+	  if (HAVE_UMFPACK_WITH_CHOLMOD)
+		list (APPEND UMFPACK_EXTRA_LIBS ${CHOLMOD_LIBRARIES})
+	  else (HAVE_UMFPACK_WITH_CHOLMOD)
+		set (UMFPACK_EXTRA_LIBS "-NOTFOUND")
+	  endif (HAVE_UMFPACK_WITH_CHOLMOD)
+	else (CHOLMOD_LIBRARIES)
+	  # if we don't have cholmod, then we certainly cannot have umfpack with cholmod
+	  set (UMFPACK_EXTRA_LIBS "-NOTFOUND")
+	endif (CHOLMOD_LIBRARIES)
+  endif (HAVE_UMFPACK_WITHOUT_CHOLMOD)
+  list (APPEND UMFPACK_LIBRARIES ${UMFPACK_EXTRA_LIBS})
+  list (REVERSE UMFPACK_LIBRARIES)
+  list (REMOVE_DUPLICATES UMFPACK_LIBRARIES)
+  list (REVERSE UMFPACK_LIBRARIES)
+endif (UMFPACK_LIBRARY)
+
+# don't reset these sets; if two packages request SuiteSparse with
+# different modules, we want the sets to be merged
+#set (SuiteSparse_LIBRARIES "")
+#set (SuiteSparse_INCLUDE_DIRS "")
+
+# determine which modules were found based on whether all dependencies
+# were satisfied; create a list of ALL modules (specified) that was found
+# (to be included in one swoop in CMakeLists.txt)
+set (SuiteSparse_FOUND TRUE)
+foreach (module IN LISTS SuiteSparse_FIND_COMPONENTS)
+  string (TOUPPER ${module} MODULE)
+  set (SuiteSparse_${MODULE}_FOUND TRUE)
+  foreach (file IN LISTS ${MODULE}_INCLUDE_DIRS ${MODULE}_LIBRARIES)
+	if (NOT EXISTS ${file})
+	  set (SuiteSparse_${MODULE}_FOUND FALSE)
+	endif (NOT EXISTS ${file})
+  endforeach (file)
+  if (NOT SuiteSparse_${MODULE}_FOUND)
+	set (SuiteSparse_FOUND FALSE)
+	# use empty string instead of zero, so it can be tested with #ifdef
+	# as well as #if in the source code
+	set (HAVE_SUITESPARSE_${MODULE}_H "" CACHE INT "Is ${module} header present?")
+  else (NOT SuiteSparse_${MODULE}_FOUND)
+	set (HAVE_SUITESPARSE_${MODULE}_H 1 CACHE INT "Is ${module} header present?")
+	list (APPEND SuiteSparse_LIBRARIES "${${MODULE}_LIBRARIES}")
+	list (APPEND SuiteSparse_LINKER_FLAGS "${${MODULE}_LINKER_FLAGS}")
+	list (APPEND SuiteSparse_INCLUDE_DIRS "${${MODULE}_INCLUDE_DIRS}")
+  endif (NOT SuiteSparse_${MODULE}_FOUND)
+  mark_as_advanced (HAVE_SUITESPARSE_${MODULE}_H)
+  mark_as_advanced (${MODULE}_INCLUDE_DIR)
+  mark_as_advanced (${MODULE}_LIBRARY)
+endforeach (module)
+
+if (SuiteSparse_INCLUDE_DIRS)
+  list (REMOVE_DUPLICATES SuiteSparse_INCLUDE_DIRS)
+endif (SuiteSparse_INCLUDE_DIRS)
+if (SuiteSparse_LIBRARIES)
+  list (REVERSE SuiteSparse_LIBRARIES)
+  list (REMOVE_DUPLICATES SuiteSparse_LIBRARIES)
+  list (REVERSE SuiteSparse_LIBRARIES)
+endif (SuiteSparse_LIBRARIES)
+
+# print a message to indicate status of this package
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (SuiteSparse
+  DEFAULT_MSG
+  SuiteSparse_LIBRARIES
+  SuiteSparse_INCLUDE_DIRS
+  )
+
+# add these after checking to not pollute the message output (checking for
+# BLAS and LAPACK is REQUIRED so if they are not found, we'll have failed
+# already; suitesparseconfig is "optional" anyway)
+list (APPEND SuiteSparse_LIBRARIES ${SuiteSparse_EXTRA_LIBS})
diff --git a/cmake/Modules/IFEMOptions.cmake b/cmake/Modules/IFEMOptions.cmake
index bc14ce91..6b9cdb1c 100644
--- a/cmake/Modules/IFEMOptions.cmake
+++ b/cmake/Modules/IFEMOptions.cmake
@@ -4,11 +4,13 @@ OPTION(IFEM_USE_SUPERLU        "Compile with SuperLU support?"        ON)
 OPTION(IFEM_USE_SUPERLU_MT     "Compile with SuperLU_MT support?"     OFF)
 OPTION(IFEM_USE_PETSC          "Compile with PETSc support?"          ON)
 OPTION(IFEM_USE_MPI            "Compile with MPI support?"            OFF)
+OPTION(IFEM_USE_ISTL           "Compile with dune-istl support?"      ON)
 OPTION(IFEM_USE_SLEPC          "Compile with SLEPc support?"          OFF)
 OPTION(IFEM_USE_SPR            "Compile with SPR support?"            OFF)
 OPTION(IFEM_USE_SAMG           "Compile with SAMG support?"           OFF)
 OPTION(IFEM_USE_HDF5           "Compile with HDF5 support?"           ON)
 OPTION(IFEM_USE_VTFWRITER      "Compile with VTFWriter support?"      ON)
+OPTION(IFEM_USE_UMFPACK        "Compile with UMFPACK support?"        ON)
 OPTION(IFEM_AS_SUBMODULE       "Compile IFEM as a submodule of apps?" OFF)
 OPTION(IFEM_WHOLE_PROG_OPTIM   "Compile IFEM with link-time optimizations?" OFF)
 OPTION(IFEM_TEST_MEMCHECK      "Run tests through valgrind?"          OFF)
diff --git a/cmake/Scripts/IFEMTesting.cmake b/cmake/Scripts/IFEMTesting.cmake
index 1a574806..1e11e4f4 100644
--- a/cmake/Scripts/IFEMTesting.cmake
+++ b/cmake/Scripts/IFEMTesting.cmake
@@ -53,6 +53,12 @@ macro(IFEM_add_unittests IFEM_PATH)
   file(GLOB TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C;${IFEM_PATH}/src/SIM/Test/*.C)
   if(NOT PETSC_FOUND)
     list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestPETScMatrix.C)
+    list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLPETScMatrix.C)
+  endif()
+
+  if(NOT ISTL_FOUND)
+    list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLMatrix.C)
+    list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLPETScMatrix.C)
   endif()
 
   IFEM_add_test_app("${TEST_SOURCES}"
@@ -66,6 +72,9 @@ macro(IFEM_add_unittests IFEM_PATH)
     if(PETSC_FOUND)
       list(APPEND TEST_SRCS_MPI ${IFEM_PATH}/src/LinAlg/Test/MPI/TestPETScMatrix.C)
     endif()
+    if(ISTL_FOUND)
+      list(APPEND TEST_SRCS_MPI ${IFEM_PATH}/src/LinAlg/Test/MPI/TestISTLMatrix.C)
+    endif()
     add_executable(IFEM-MPI-test EXCLUDE_FROM_ALL
                    ${IFEM_PATH}/src/IFEM-test.C ${TEST_SRCS_MPI})
     target_link_libraries(IFEM-MPI-test ${IFEM_LIBRARIES} ${IFEM_DEPLIBS} gtest)
diff --git a/src/ASM/AlgEqSystem.C b/src/ASM/AlgEqSystem.C
index e24667aa..4ffa8e78 100644
--- a/src/ASM/AlgEqSystem.C
+++ b/src/ASM/AlgEqSystem.C
@@ -52,6 +52,8 @@ bool AlgEqSystem::init (SystemMatrix::Type mtype, const LinSolParams* spar,
   SystemVector::Type vtype = SystemVector::STD;
   if (mtype == SystemMatrix::PETSC)
     vtype = SystemVector::PETSC;
+  if (mtype == SystemMatrix::ISTL)
+    vtype = SystemVector::ISTL;
 
   for (i = 0; i < b.size(); i++)
     if (!b[i])
diff --git a/src/IFEM.C b/src/IFEM.C
index 3a094d16..c16aad58 100644
--- a/src/IFEM.C
+++ b/src/IFEM.C
@@ -114,6 +114,13 @@ int IFEM::Init (int arg_c, char** arg_v, const char* title)
   std::cout <<"disabled";
 #endif
 
+  std::cout <<"\n       ISTL support: ";
+#if HAS_ISTL
+  std::cout <<"enabled (v"<< ISTL_VERSION <<")";
+#else
+  std::cout <<"disabled";
+#endif
+
   std::cout <<"\n        VTF support: ";
 #if HAS_VTFAPI == 2
   std::cout <<"enabled (v2)";
diff --git a/src/LinAlg/ISTLMatrix.C b/src/LinAlg/ISTLMatrix.C
new file mode 100644
index 00000000..f5442d76
--- /dev/null
+++ b/src/LinAlg/ISTLMatrix.C
@@ -0,0 +1,253 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ISTLMatrix.C
+//!
+//! \date Mar 2 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Representation of the system matrix in ISTL format.
+//!
+//==============================================================================
+
+#include "ISTLMatrix.h"
+#include "LinSolParams.h"
+#include "LinAlgInit.h"
+#include "ProcessAdm.h"
+#include "SIMenums.h"
+#include "ASMstruct.h"
+#include "Profiler.h"
+#include "SAMpatch.h"
+#include "DomainDecomposition.h"
+
+
+ISTLVector::ISTLVector(const ProcessAdm& padm) : adm(padm)
+{
+  LinAlgInit::increfs();
+}
+
+
+ISTLVector::ISTLVector(const ProcessAdm& padm, size_t n) : adm(padm)
+{
+  x.resize(n);
+  LinAlgInit::increfs();
+}
+
+
+ISTLVector::ISTLVector(const ProcessAdm& padm, const Real* values, size_t n) : adm(padm)
+{
+  x.resize(n);
+  this->restore(values);
+  LinAlgInit::increfs();
+}
+
+
+ISTLVector::ISTLVector(const ISTLVector& vec) : adm(vec.adm)
+{
+  x = vec.x;
+  LinAlgInit::increfs();
+}
+
+
+ISTLVector::~ISTLVector()
+{
+  LinAlgInit::decrefs();
+}
+
+
+void ISTLVector::init(Real value)
+{
+  StdVector::init(value);
+  x = value;
+}
+
+size_t ISTLVector::dim() const
+{
+  return x.size();
+}
+
+
+void ISTLVector::redim(size_t n)
+{
+  x.resize(n);
+  StdVector::redim(n);
+}
+
+
+bool ISTLVector::beginAssembly()
+{
+  for (size_t i = 0; i < size(); ++i)
+    x[i] = (*this)[i];
+
+  return true;
+}
+
+
+bool ISTLVector::endAssembly()
+{
+  return true;
+}
+
+
+Real ISTLVector::L1norm() const
+{
+  return x.one_norm();
+}
+
+
+Real ISTLVector::L2norm() const
+{
+  return x.two_norm();
+}
+
+
+Real ISTLVector::Linfnorm() const
+{
+  return x.infinity_norm();
+}
+
+
+ISTLMatrix::ISTLMatrix (const ProcessAdm& padm, const LinSolParams& spar,
+                        LinAlg::LinearSystemType ltype) :
+ SparseMatrix(SUPERLU, 1), adm(padm),
+ solParams(spar, adm), linsysType(ltype)
+{
+  LinAlgInit::increfs();
+
+  setParams = true;
+  nLinSolves = 0;
+}
+
+
+ISTLMatrix::ISTLMatrix (const ISTLMatrix& B) :
+  adm(B.adm), solParams(B.solParams.get(), B.adm)
+{
+  A = B.A;
+
+  LinAlgInit::increfs();
+
+  setParams = true;
+}
+
+
+ISTLMatrix::~ISTLMatrix ()
+{
+  LinAlgInit::decrefs();
+}
+
+
+void ISTLMatrix::initAssembly (const SAM& sam, bool b)
+{
+  SparseMatrix::initAssembly(sam, b);
+  SparseMatrix::preAssemble(sam, b);
+
+  std::vector<std::set<int>> dofc;
+  sam.getDofCouplings(dofc);
+
+  // Set correct number of rows and columns for matrix.
+  size_t sum = 0;
+  for (const auto& it : dofc)
+    sum += it.size();
+
+  A.setSize(rows(), cols(), sum);
+  A.setBuildMode(ISTL::Mat::random);
+
+  for (size_t i = 0; i < dofc.size(); ++i)
+    A.setrowsize(i,dofc[i].size());
+  A.endrowsizes();
+
+  for (size_t i = 0; i < dofc.size(); ++i)
+    for (const auto& it : dofc[i])
+      A.addindex(i, it-1);
+
+  A.endindices();
+
+  A = 0;
+}
+
+bool ISTLMatrix::beginAssembly()
+{
+  for (size_t j = 0; j < cols(); ++j)
+    for (int i = IA[j]; i < IA[j+1]; ++i)
+      A[JA[i]][j] = SparseMatrix::A[i];
+
+  return true;
+}
+
+
+bool ISTLMatrix::endAssembly()
+{
+  return true;
+}
+
+
+void ISTLMatrix::init ()
+{
+  SparseMatrix::init();
+
+  // Set all matrix elements to zero
+  A = 0;
+}
+
+
+
+bool ISTLMatrix::solve (SystemVector& B, bool newLHS, Real*)
+{
+  if (!pre)
+    std::tie(solver, pre, op) = solParams.setupPC(A);
+
+  ISTLVector* Bptr = dynamic_cast<ISTLVector*>(&B);
+  if (!Bptr || !solver || !pre)
+    return false;
+
+  try {
+    Dune::InverseOperatorResult r;
+    ISTL::Vec b(Bptr->getVector());
+    Bptr->getVector() = 0;
+    solver->apply(Bptr->getVector(), b, r);
+  } catch (Dune::ISTLError& e) {
+    std::cerr << "ISTL exception " << e << std::endl;
+    return false;
+  }
+
+  for (size_t i = 0; i < rows(); ++i)
+    (*Bptr)(i+1) = Bptr->getVector()[i];
+
+  return true;
+}
+
+
+bool ISTLMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS)
+{
+  if (!pre)
+    std::tie(solver, pre, op) = solParams.setupPC(A);
+
+  const ISTLVector* Bptr = dynamic_cast<const ISTLVector*>(&b);
+  if (!Bptr || ! solver || !pre)
+    return false;
+
+  ISTLVector* Xptr = dynamic_cast<ISTLVector*>(&x);
+  if (!Xptr)
+    return false;
+
+  try {
+    Dune::InverseOperatorResult r;
+    solver->apply(Xptr->getVector(),
+                  const_cast<ISTL::Vec&>(Bptr->getVector()), r);
+  } catch (Dune::ISTLError& e) {
+    std::cerr << "ISTL exception " << e << std::endl;
+    return false;
+  }
+
+  for (size_t i = 0; i < rows(); ++i)
+    (*Xptr)(i+1) = Xptr->getVector()[i];
+
+  return true;
+}
+
+
+Real ISTLMatrix::Linfnorm () const
+{
+  return A.infinity_norm();
+}
diff --git a/src/LinAlg/ISTLMatrix.h b/src/LinAlg/ISTLMatrix.h
new file mode 100644
index 00000000..1e68d6a7
--- /dev/null
+++ b/src/LinAlg/ISTLMatrix.h
@@ -0,0 +1,163 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ISTLMatrix.h
+//!
+//! \date Mar 2 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Representation of the system matrix in ISTL format with interface
+//! to ISTL routines for solving linear equation systems.
+//!
+//==============================================================================
+
+#ifndef _ISTL_MATRIX_H
+#define _ISTL_MATRIX_H
+
+#include "SystemMatrix.h"
+#include "SparseMatrix.h"
+#include "ISTLSolParams.h"
+#include "LinAlgenums.h"
+#include <memory>
+
+class LinSolParams;
+class SAMpatch;
+
+
+/*!
+  \brief Class for representing the system vector in ISTL format.
+*/
+
+class ISTLVector : public StdVector
+{
+public:
+  //! \brief Constructor creating an empty vector.
+  ISTLVector(const ProcessAdm& padm);
+  //! \brief Constructor creating a vector of length \a n.
+  ISTLVector(const ProcessAdm& padm, size_t n);
+  //! \brief Constructor creating a vector from an array.
+  ISTLVector(const ProcessAdm& padm, const Real* values, size_t n);
+  //! \brief Copy constructor.
+  ISTLVector(const ISTLVector& vec);
+  //! \brief Destructor.
+  virtual ~ISTLVector();
+#endif
+
+  //! \brief Returns the vector type.
+  virtual Type getType() const { return ISTL; }
+
+  //! \brief Initializes the vector to a given scalar value.
+  virtual void init(Real value = Real(0));
+
+  //! \brief Returns the dimension of the system vector.
+  virtual size_t dim() const;
+
+  //! \brief Sets the dimension of the system vector.
+  virtual void redim(size_t n);
+
+  //! \brief Begins communication step needed in parallel vector assembly.
+  //! \details Must be called together with endAssembly after vector assembly
+  //! is completed on each processor and before the linear system is solved.
+  virtual bool beginAssembly();
+
+  //! \brief Ends communication step needed in parallel vector assembly.
+  //! \details Must be called together with beginAssembly after vector assembly
+  //! is completed on each processor and before the linear system is solved.
+  virtual bool endAssembly();
+
+  //! \brief L1-norm of vector.
+  virtual Real L1norm() const;
+
+  //! \brief L2-norm of vector.
+  virtual Real L2norm() const;
+
+  //! \brief Linfinity-norm of vector.
+  virtual Real Linfnorm() const;
+
+  //! \brief Return associated process administrator
+  const ProcessAdm& getAdm() const { return adm; }
+  typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vec;
+
+  Vec& getVector() { return x; }
+  const Vec& getVector() const { return x; }
+
+protected:
+  ISTL::Vec x; //!< ISTL vector
+  const ProcessAdm& adm;  //!< Process administrator
+};
+
+
+/*!
+  \brief Class for representing the system matrix in ISTL format.
+*/
+
+class ISTLMatrix : public SparseMatrix
+{
+public:
+  //! \brief Constructor.
+  ISTLMatrix(const ProcessAdm& padm, const LinSolParams& spar,
+             LinAlg::LinearSystemType ltype);
+  //! \brief Copy constructor.
+  ISTLMatrix(const ISTLMatrix& A);
+  //! \brief The destructor frees the dynamically allocated arrays.
+  virtual ~ISTLMatrix();
+
+  //! \brief Returns the matrix type.
+  virtual Type getType() const { return ISTL; }
+
+  //! \brief Returns the dimension of the system matrix.
+  virtual size_t dim(int = 1) const { return 0; }
+
+  //! \brief Creates a copy of the system matrix and returns a pointer to it.
+  virtual SystemMatrix* copy() const { return new ISTLMatrix(*this); }
+
+  //! \brief Initializes the element assembly process.
+  //! \details Must be called once before the element assembly loop.
+  //! The PETSc data structures are initialized and the all symbolic operations
+  //! that are needed before the actual assembly can start are performed.
+  //! \param[in] sam Auxiliary data describing the FE model topology, etc.
+  virtual void initAssembly(const SAM& sam, bool);
+
+  //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
+  virtual void init();
+
+  //! \brief Begins communication step needed in parallel matrix assembly.
+  //! \details Must be called together with endAssembly after matrix assembly
+  //! is completed on each processor and before the linear system is solved.
+  virtual bool beginAssembly();
+  //! \brief Ends communication step needed in parallel matrix assembly.
+  //! \details Must be called together with beginAssembly after matrix assembly
+  //! is completed on each processor and before the linear system is solved.
+  virtual bool endAssembly();
+
+  //! \brief Solves the linear system of equations for a given right-hand-side.
+  //! \param B Right-hand-side vector on input, solution vector on output
+  //! \param[in] newLHS \e true if the left-hand-side matrix has been updated
+  virtual bool solve(SystemVector& B, bool newLHS, Real*);
+
+  //! \brief Solves the linear system of equations for a given right-hand-side.
+  //! \param[in] B Right-hand-side vector
+  //! \param[out] x Solution vector
+  //! \param[in] newLHS \e true if the left-hand-side matrix has been updated
+  virtual bool solve(const SystemVector& B, SystemVector& x, bool newLHS);
+
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual Real Linfnorm() const;
+
+  //! \brief Returns the ISTL matrix (for assignment).
+  virtual ISTL::Mat& getMatrix() { return A; }
+  //! \brief Returns the ISTL matrix (for read access).
+  virtual const ISTL::Mat& getMatrix() const { return A; }
+
+protected:
+  ISTL::Mat A; //!< The actual ISTL matrix
+  std::unique_ptr<ISTL::Operator> op; //!< The matrix adapter
+  std::unique_ptr<ISTL::InverseOperator> solver; //!< Solver to use
+  std::unique_ptr<ISTL::Preconditioner> pre; //!< Preconditioner to use
+  const ProcessAdm&   adm;             //!< Process administrator
+  ISTLSolParams       solParams;       //!< Linear solver parameters
+  bool                setParams;       //!< If linear solver parameters are set
+  int                 nLinSolves;      //!< Number of linear solves
+  LinAlg::LinearSystemType linsysType; //!< Linear system type
+};
diff --git a/src/LinAlg/ISTLSolParams.C b/src/LinAlg/ISTLSolParams.C
new file mode 100644
index 00000000..d096bf6b
--- /dev/null
+++ b/src/LinAlg/ISTLSolParams.C
@@ -0,0 +1,510 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ISTLSolParams.C
+//!
+//! \date Mar 10 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Linear solver parameters for ISTL.
+//!
+//==============================================================================
+
+#include "ISTLSolParams.h"
+#include "ASMstruct.h"
+#include "LinSolParams.h"
+#include "ProcessAdm.h"
+#include "SAMpatch.h"
+#include <dune/istl/overlappingschwarz.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/matrixmatrix.hh>
+#include <dune/common/version.hh>
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+#include <dune/istl/paamg/fastamg.hh>
+#include <dune/istl/paamg/twolevelmethod.hh>
+#endif
+
+
+namespace ISTL {
+
+BlockPreconditioner::BlockPreconditioner(const ISTL::Mat& A,
+                                         const DomainDecomposition& dd_,
+                                         const std::string& schurType) :
+  dd(dd_)
+{
+  blocks.resize(dd.getNoBlocks()*dd.getNoBlocks()+1);
+  size_t k = 0;
+  for (size_t i = 0; i < dd.getNoBlocks(); ++i)
+    for (size_t j = 0; j < dd.getNoBlocks(); ++j)
+      extractBlock(blocks[k++], A, dd.getBlockEqs(i), dd.getBlockEqs(j), dd.getG2LEQ(j));
+
+  // Scale blocks[1] with inverse diagonal of blocks[0]
+  blocks.back() = blocks[1];
+  for (auto row = blocks.back().begin(); row != blocks.back().end(); ++row)
+    for (auto it = row->begin(); it != row->end(); ++it)
+      *it /= blocks[0][row.index()][row.index()];
+
+  // blocks[1] = D - C*diag(A)^-1*B
+  ISTL::Mat SP;
+  Dune::matMultMat(SP, blocks[2], blocks.back());
+  blocks[2] = blocks[1];
+  subtractMatrices(blocks[1], blocks[3], SP);
+
+  // clear memory for unused blocks
+  blocks.resize(schurType == "diag" ? 2 : 3);
+
+  blockPre.resize(dd.getNoBlocks());
+  blockOp.resize(dd.getNoBlocks());
+}
+
+
+void BlockPreconditioner::pre(ISTL::Vec& x, ISTL::Vec& b)
+{
+  ISTL::Vec tempx, tempb;
+  for (size_t block = 0; block < dd.getNoBlocks(); ++block) {
+    ISTL::Vec tempx, tempb;
+    tempx.resize(dd.getBlockEqs(block).size());
+    tempb.resize(dd.getBlockEqs(block).size());
+    size_t i = 0;
+    for (auto& it : dd.getBlockEqs(block)) {
+      tempx[i] = x[it-1];
+      tempb[i++] = b[it-1];
+    }
+
+    blockPre[block]->pre(tempx, tempb);
+
+    i = 0;
+    for (auto& it : dd.getBlockEqs(block)) {
+      x[it-1] = tempx[i];
+      b[it-1] = tempb[i++];
+    }
+  }
+}
+
+
+void BlockPreconditioner::apply(ISTL::Vec& v, const ISTL::Vec& d)
+{
+  // backwards due to upper schur complement
+  ISTL::Vec tempx, tempb;
+  for (int block = dd.getNoBlocks()-1; block >= 0; --block) {
+    tempb.resize(dd.getBlockEqs(block).size());
+    size_t i = 0;
+    for (auto& it : dd.getBlockEqs(block))
+      tempb[i++] = d[it-1];
+
+    if (block == 0 && blocks.size() > 2)
+      blocks.back().mmv(tempx, tempb);
+
+    tempx.resize(dd.getBlockEqs(block).size());
+    tempx = 0;
+    blockPre[block]->apply(tempx, tempb);
+
+    i = 0;
+    for (auto& it : dd.getBlockEqs(block))
+      v[it-1] = tempx[i++];
+  }
+}
+
+
+void BlockPreconditioner::post(ISTL::Vec& x)
+{
+  ISTL::Vec tempx, tempb;
+  for (size_t block = 0; block < dd.getNoBlocks(); ++block) {
+    ISTL::Vec tempx;
+    tempx.resize(dd.getBlockEqs(block).size());
+    size_t i = 0;
+    for (auto& it : dd.getBlockEqs(block))
+      tempx[i++] = x[it-1];
+
+    blockPre[block]->post(tempx);
+
+    i = 0;
+    for (auto& it : dd.getBlockEqs(block))
+      x[it-1] = tempx[i++];
+  }
+}
+
+
+void BlockPreconditioner::extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
+                                       const std::set<int>& eqs_row1,
+                                       const std::set<int>& eqs_col1,
+                                       const std::map<int,int>& eqs_col_g2l)
+{
+  std::vector<int> eqs_row(eqs_row1.begin(), eqs_row1.end());
+  std::vector<int> eqs_col(eqs_col1.begin(), eqs_col1.end());
+  size_t sum=0;
+  std::vector<std::set<int>> adj;
+  adj.resize(eqs_row.size());
+  size_t i = 0;
+
+  for (auto& it3 : eqs_row) {
+    auto it = A.begin()+it3-1;
+    for (auto it2 = it->begin(); it2 != it->end(); ++it2) {
+      auto itc = eqs_col_g2l.find(it2.index()+1);
+      if (itc != eqs_col_g2l.end()) {
+        adj[i].insert(itc->second-1);
+        ++sum;
+      }
+    }
+    ++i;
+  }
+
+  B.setSize(eqs_row.size(), eqs_col.size(), sum);
+  B.setBuildMode(ISTL::Mat::random);
+
+  for (size_t i = 0; i < adj.size(); i++)
+    B.setrowsize(i,adj[i].size());
+  B.endrowsizes();
+
+  for (size_t i = 0; i < adj.size(); i++)
+    for (auto& it : adj[i])
+      B.addindex(i, it);
+
+  B.endindices();
+  B = 0;
+  for (auto it = B.begin(); it != B.end(); ++it)
+    for (auto it2 = it->begin(); it2 != it->end(); ++it2)
+      *it2 = A[eqs_row[it.index()]-1][eqs_col[it2.index()]-1];
+}
+
+
+void BlockPreconditioner::subtractMatrices(ISTL::Mat& A, const ISTL::Mat& B,
+                                           const ISTL::Mat& C)
+{
+  size_t sum=0;
+  std::vector<std::set<int>> adj;
+  adj.resize(B.N());
+  for (auto row = B.begin(); row != B.end(); ++row) {
+    for (auto it = row->begin(); it != row->end(); ++it)
+      adj[row.index()].insert(it.index());
+  }
+  for (auto row = C.begin(); row != C.end(); ++row) {
+    for (auto it = row->begin(); it != row->end(); ++it)
+      adj[row.index()].insert(it.index());
+    sum += adj[row.index()].size();
+  }
+
+  A.setSize(B.N(), B.M(), sum);
+  A.setBuildMode(ISTL::Mat::random);
+
+  for (size_t i = 0; i < adj.size(); i++)
+    A.setrowsize(i,adj[i].size());
+  A.endrowsizes();
+
+  for (size_t i = 0; i < adj.size(); i++)
+    for (auto& it : adj[i])
+      A.addindex(i, it);
+
+  A.endindices();
+  A = 0;
+  for (auto row = B.begin(); row != B.end(); ++row)
+    for (auto it = row->begin(); it != row->end(); ++it)
+      A[row.index()][it.index()] = *it;
+
+  for (auto row = C.begin(); row != C.end(); ++row)
+    for (auto it = row->begin(); it != row->end(); ++it)
+      A[row.index()][it.index()] -= *it;
+}
+
+
+} // namespace ISTL
+
+
+/*! \brief Helper template for setting up a solver with the appropriate
+           preconditioner type.
+    \details We cannot instance using dynamic polymorphism, the solver need
+             access to the real type for the preconditioner. We can however
+             call the solver in the interface class scope afterwards.
+ */
+template<class Prec>
+static ISTL::InverseOperator* setupWithPreType(const LinSolParams& solParams,
+                                               ISTL::Operator& op,
+                                               Prec& pre)
+{
+  std::string type = solParams.getStringValue("type");
+  double rtol = solParams.getDoubleValue("rtol");
+  int maxits = solParams.getIntValue("maxits");
+  int verbosity = solParams.getIntValue("verbosity");
+  int restart = solParams.getIntValue("gmres_restart_iterations");
+  if (type == "bcgs")
+    return new Dune::BiCGSTABSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
+  else if (type == "cg")
+    return new Dune::CGSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
+  else if (type == "minres")
+    return new Dune::MINRESSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
+  else if (type == "gmres")
+    return new Dune::RestartedGMResSolver<ISTL::Vec>(op, pre, rtol, restart,
+                                                     maxits, verbosity);
+
+  return nullptr;
+}
+
+
+/*! \brief Helper template for setting up an AMG preconditioner with a given smoother */
+
+template<class Smoother>
+static ISTL::Preconditioner* setupAMG(const LinSolParams& params,
+                                      size_t block, ISTL::Operator& op,
+                                      std::unique_ptr<ISTL::InverseOperator>* solver)
+{
+  // The coupling metric used in the AMG
+  typedef Dune::Amg::FirstDiagonal CouplingMetric;
+  // The coupling criterion used in the AMG
+  typedef Dune::Amg::SymmetricCriterion<ISTL::Mat, CouplingMetric> CritBase;
+  // The coarsening criterion used in the AMG
+  typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
+  Criterion crit;
+  typedef typename Dune::Amg::AMG<ISTL::Operator, ISTL::Vec, Smoother> AMG;
+  typename AMG::SmootherArgs args;
+  args.relaxationFactor = 1.0;
+  args.iterations = std::max(1, params.getBlock(block).getIntValue("multigrid_no_smooth"));
+
+  if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
+    crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
+
+  if (params.getBlock(block).hasValue("multigrid_no_smooth")) {
+    int val = params.getBlock(block).getIntValue("multigrid_no_smooth");
+    crit.setNoPreSmoothSteps(val);
+    crit.setNoPostSmoothSteps(val);
+  }
+
+  if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
+    crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
+
+  auto result = new AMG(op, crit, args);
+  if (solver)
+    solver->reset(setupWithPreType(params, op, *result));
+
+  return result;
+}
+
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+/*! \brief Helper template for setting up a AMG preconditioner
+ *!        with fine smoother differing from the other levels */
+
+template<class Smoother, class FineSmoother>
+static ISTL::Preconditioner* setupAMG2_full(const LinSolParams& params, size_t block,
+                                            ISTL::Operator& op,
+                                            std::unique_ptr<ISTL::InverseOperator>* solver,
+                                            FineSmoother* fsmooth)
+{
+  typedef Dune::Amg::FirstDiagonal CouplingMetric;
+  typedef Dune::Amg::SymmetricCriterion<ISTL::Mat,CouplingMetric> CritBase;
+  typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
+  typedef Dune::Amg::AggregationLevelTransferPolicy<ISTL::Operator,Criterion> TransferPolicy;
+  typedef Dune::Amg::OneStepAMGCoarseSolverPolicy<ISTL::Operator,Smoother,Criterion> CoarsePolicy;
+  typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+  typedef typename Dune::Amg::TwoLevelMethod<ISTL::Operator, CoarsePolicy, FineSmoother> AMG2;
+
+  SmootherArgs args;
+  args.relaxationFactor = 1.0;
+
+  Criterion crit;
+  if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
+    crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
+
+  if (params.getBlock(block).hasValue("multigrid_no_smooth")) {
+    int val = params.getBlock(block).getIntValue("multigrid_no_smooth");
+    crit.setNoPreSmoothSteps(val);
+    crit.setNoPostSmoothSteps(val);
+  }
+
+  if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
+    crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
+
+  CoarsePolicy coarsePolicy(args, crit);
+  TransferPolicy policy(crit);
+  Dune::shared_ptr<FineSmoother> fsp(fsmooth);
+  auto result = new AMG2(op, fsp, policy, coarsePolicy);
+  if (solver)
+    solver->reset(setupWithPreType(params, op, *result));
+
+  return result;
+}
+
+
+template<class FineSmoother>
+static ISTL::Preconditioner* setupAMG2_smoother(const LinSolParams& params, size_t block,
+                                                ISTL::Operator& op,
+                                                std::unique_ptr<ISTL::InverseOperator>* solver,
+                                                FineSmoother* fsmooth)
+{
+  std::string smoother = params.getBlock(block).getStringValue("multigrid_smoother");
+  if (smoother == "ilu")
+    return setupAMG2_full<ISTL::ILU0>(params, block, op, solver, fsmooth);
+  else if (smoother == "sor")
+    return setupAMG2_full<ISTL::SOR>(params, block, op, solver, fsmooth);
+  else if (smoother == "ssor")
+    return setupAMG2_full<ISTL::SSOR>(params, block, op, solver, fsmooth);
+  else if (smoother == "jacobi")
+    return setupAMG2_full<ISTL::GJ>(params, block, op, solver, fsmooth);
+  else {
+    std::cerr << "** ISTLSolParams ** Invalid smoother " << smoother << "." << std::endl;
+    return nullptr;
+  }
+}
+
+
+static ISTL::Preconditioner* setupAMG2(const LinSolParams& params, size_t block,
+                                       ISTL::Operator& op, std::unique_ptr<ISTL::InverseOperator>* solver)
+{
+  std::string smoother = params.getBlock(block).getStringValue("multigrid_finesmoother");
+  if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "ilu") {
+    auto fsmooth = new ISTL::ILU0(op.getmat(), 1.0);
+    return setupAMG2_smoother(params, block, op, solver, fsmooth);
+  } else  if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "sor") {
+    auto fsmooth = new ISTL::SOR(op.getmat(), 1, 1.0);
+    return setupAMG2_smoother(params, block, op, solver, fsmooth);
+  } else  if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "ssor") {
+    auto fsmooth = new ISTL::SSOR(op.getmat(), 1, 1.0);
+    return setupAMG2_smoother(params, block, op, solver, fsmooth);
+  } else  if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "jacobi") {
+    auto fsmooth = new ISTL::GJ(op.getmat(), 1, 1.0);
+    return setupAMG2_smoother(params, block, op, solver, fsmooth);
+  } else {
+    std::cerr << "** ISTLSolParams ** Invalid fine smoother " << smoother << "." << std::endl;
+    return nullptr;
+  }
+}
+#endif
+
+
+/*! \brief Conditionally setup a KSP */
+template<class Type>
+static ISTL::Preconditioner*
+setupSolver(Type* pre,
+            ISTL::Operator& op,
+            const LinSolParams& solParams,
+            std::unique_ptr<ISTL::InverseOperator>* solver)
+{
+  if (solver)
+    solver->reset(setupWithPreType(solParams, op, *pre));
+
+  return pre;
+}
+
+ISTL::Preconditioner* ISTLSolParams::setupPCInternal(ISTL::Mat& A,
+                                                     ISTL::Operator& op,
+                                                     size_t block,
+                                                     std::unique_ptr<ISTL::InverseOperator>* solver)
+{
+  std::string prec = solParams.getBlock(block).getStringValue("pc");
+  if (prec == "ilu") {
+    int fill_level = solParams.getBlock(block).getIntValue("ilu_fill_level");
+    if (fill_level == 0)
+      return setupSolver(new ISTL::ILU0(A, 1.0), op, solParams, solver);
+    else
+      return setupSolver(new ISTL::ILUn(A, fill_level, 1.0),
+                         op, solParams, solver);
+  } else if (prec == "lu")
+    return setupSolver(new ISTL::LU(new ISTL::LUType(A)),
+                       op, solParams, solver);
+  else if (prec == "sor")
+    return setupSolver(new ISTL::SOR(A, 1, 1.0), op, solParams, solver);
+  else if (prec == "ssor")
+    return setupSolver(new ISTL::SSOR(A, 1, 1.0), op, solParams, solver);
+  else if (prec == "jacobi")
+    return setupSolver(new ISTL::GJ(A, 1, 1.0), op, solParams, solver);
+ else if (prec == "gs")
+    return setupSolver(new ISTL::GS(A, 1, 1.0), op, solParams, solver);
+  else if (prec == "asm" || prec == "asmlu") {
+    size_t nx = solParams.getBlock(block).getIntValue("asm_nx");
+    nx = std::max(1ul, nx);
+    size_t ny = solParams.getBlock(block).getIntValue("asm_ny");
+    ny = std::max(1ul, ny);
+    size_t nz = solParams.getBlock(block).getIntValue("asm_nz");
+    nz = std::max(1ul, nz);
+    int overlap = solParams.getBlock(block).getIntValue("asm_overlap");
+    auto locSubdDofs = adm.dd.getSubdomains(nx, ny, nz, overlap, block);
+
+    if (prec == "asmlu") {
+      ISTL::ASMLU::subdomain_vector ddofs(locSubdDofs.size());
+      for (size_t i = 0; i < locSubdDofs.size(); ++i)
+        for (const auto& it : locSubdDofs[i])
+          ddofs[i].insert(it-1);
+
+      return setupSolver(new ISTL::ASMLU(A, ddofs), op, solParams, solver);
+    } else {
+      ISTL::ASM::subdomain_vector ddofs(locSubdDofs.size());
+      for (size_t i = 0; i < locSubdDofs.size(); ++i) {
+        for (size_t i = 0; i < locSubdDofs.size(); ++i)
+          for (const auto& it : locSubdDofs[i])
+            ddofs[i].insert(it-1);
+      }
+
+      return setupSolver(new ISTL::ASM(A, ddofs), op, solParams, solver);
+    }
+  } else if (prec == "amg" || prec == "gamg") {
+    if (solParams.getBlock(block).getStringValue("multigrid_smoother") !=
+        solParams.getBlock(block).getStringValue("multigrid_finesmoother")) {
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+      return setupAMG2(solParams, block, op, solver);
+#else
+      std::cerr << "** ISTLSolParams ** Separate fine smoother not implemented." << std::endl;
+      return nullptr;
+#endif
+    } else {
+      std::string smoother = solParams.getBlock(block).getStringValue("multigrid_smoother");
+      if (smoother.empty()) {
+        std::cerr << "** ISTLSolParams ** No smoother defined for AMG, defaulting to ILU" << std::endl;
+        smoother = "ilu";
+      }
+      if (smoother == "ilu")
+        return setupAMG<ISTL::ILU0>(solParams, block, op, solver);
+      else if (smoother == "sor")
+        return setupAMG<ISTL::SOR>(solParams, block, op, solver);
+      else if (smoother == "ssor")
+        return setupAMG<ISTL::SSOR>(solParams, block, op, solver);
+      else if (smoother == "jacobi")
+        return setupAMG<ISTL::GJ>(solParams, block, op, solver);
+    }
+  }
+
+  return nullptr;
+}
+
+
+std::tuple<std::unique_ptr<ISTL::InverseOperator>,
+           std::unique_ptr<ISTL::Preconditioner>,
+           std::unique_ptr<ISTL::Operator>> ISTLSolParams::setupPC(ISTL::Mat& A)
+{
+  std::unique_ptr<ISTL::InverseOperator> solver;
+  std::unique_ptr<ISTL::Preconditioner> pre;
+  std::unique_ptr<ISTL::Operator> op;
+
+  if (solParams.getNoBlocks() > 2) {
+    std::cerr << "*** ISTLSolParams ** More than two blocks are not implemented." << std::endl;
+    return std::make_tuple(nullptr, nullptr, nullptr);
+  }
+
+  if (solParams.getNoBlocks() > 1) {
+    op.reset(new ISTL::Operator(A));
+    ISTL::BlockPreconditioner* bpre = new ISTL::BlockPreconditioner(A, adm.dd, solParams.getStringValue("schur"));
+    pre.reset(bpre);
+    for (size_t i = 0; i < adm.dd.getNoBlocks(); ++i)
+      bpre->getBlockPre(i).reset(setupPCInternal(bpre->getBlock(i),
+                                                 bpre->getBlockOp(i), i, nullptr));
+    solver.reset(setupWithPreType(solParams, *op, *bpre));
+  } else {
+#ifdef HAVE_MPI
+    if (adm.isParallel()) {
+        Dune::IndexInfoFromGrid<int, int> index;
+      for (size_t i = 0; i < adm.dd.getMLGEQ().size(); ++i)
+        index.addLocalIndex(std::make_tuple(adm.dd.getGlobalEq(i+1), i, 0));
+      Dune::OwnerOverlapCopyCommunication<int,int> comm(index, *adm.getCommunicator());
+      //ISTL::ParMatrixAdapter* nop = new ISTL::ParMatrixAdapter(A, comm);
+      //op.reset(nop);
+      //pre.reset(setupPCInternal(A, *nop, 0, &solver));
+    } else
+#endif
+    {
+      ISTL::Operator* nop = new ISTL::Operator(A);
+      op.reset(nop);
+      pre.reset(setupPCInternal(A, *nop, 0, &solver));
+    }
+  }
+
+  return std::make_tuple(std::move(solver), std::move(pre), std::move(op));
+}
diff --git a/src/LinAlg/ISTLSolParams.h b/src/LinAlg/ISTLSolParams.h
new file mode 100644
index 00000000..239defcc
--- /dev/null
+++ b/src/LinAlg/ISTLSolParams.h
@@ -0,0 +1,258 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ISTLSolParams.h
+//!
+//! \date Mar 10 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Linear solver parameters for ISTL matrices.
+//! \details Includes linear solver method, preconditioner
+//! and convergence criteria.
+//!
+//==============================================================================
+
+#ifndef _ISTL_SOLPARAMS_H
+#define _ISTL_SOLPARAMS_H
+
+#include "ISTLSupport.h"
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvercategory.hh>
+
+#include <iostream>
+#include <set>
+#include <string>
+#include <vector>
+
+
+class DomainDecomposition;
+class LinSolParams;
+class ProcessAdm;
+
+
+/*! This implements a Schur-decomposition based preconditioner for the
+ *  block system
+ *  [A   B]
+ *  [C   D]
+ *
+ *  The preconditioner is
+ *  [Apre B]
+ *  [     P]
+ *  Here Apre is some preconditioner for A and P some preconditioner for
+ *  S = D - C*diag(A)^-1*B. The B block may be dropped.
+!*/
+
+namespace ISTL {
+
+class BlockPreconditioner : public Preconditioner {
+public:
+  // define the category
+  enum {
+    //! \brief The category the preconditioner is part of.
+    category=Dune::SolverCategory::sequential
+  };
+
+  //! \brief Constructor
+  //! \param[in] A The system matrix
+  //! \param[in] dd Domain decomposition
+  BlockPreconditioner(const ISTL::Mat& A, const DomainDecomposition& dd_,
+                      const std::string& schurType);
+
+  //! \brief Destructor
+  virtual ~BlockPreconditioner()
+  {}
+
+  //! \brief Preprocess preconditioner
+  virtual void pre(ISTL::Vec& x, ISTL::Vec& b);
+
+  //! \brief Applies the preconditioner
+  //! \param[out] v The resulting vector
+  //! \param[in] d The vector to apply the preconditioner to
+  virtual void apply(ISTL::Vec& v, const ISTL::Vec& d);
+
+  //! \brief Post-process function
+  virtual void post(ISTL::Vec& x);
+
+  //! \brief Obtain reference to a block preconditioner
+  std::unique_ptr<ISTL::Preconditioner>& getBlockPre(size_t block)
+  {
+    return blockPre[block];
+  }
+
+  //! \brief Obtain reference to a block matrix
+  ISTL::Mat& getBlock(size_t block)
+  {
+    return blocks[block];
+  }
+
+  //! \brief Obtain matrix adaptor for a block matrix
+  ISTL::Operator& getBlockOp(size_t block)
+  {
+    if (!blockOp[block])
+      blockOp[block].reset(new ISTL::Operator(blocks[block]));
+
+    return *blockOp[block];
+  }
+
+protected:
+  //! \brief Build block from block equations.
+  static void extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
+                           const std::set<int>& eqs_row,
+                           const std::set<int>& eqs_col,
+                           const std::map<int,int>& eqs_col_g2l);
+
+  //! \brief Find A = B - C
+  static void subtractMatrices(ISTL::Mat& A,
+                               const ISTL::Mat& B,
+                               const ISTL::Mat& C);
+
+  std::vector<std::unique_ptr<ISTL::Preconditioner>> blockPre; //!< The preconditioners
+  std::vector<std::unique_ptr<ISTL::Operator>> blockOp; //!< The matrix adaptors for the blocks
+  std::vector<ISTL::Mat> blocks; //!< Matrix blocks
+  const DomainDecomposition& dd; //!< Domain decomposition
+};
+
+#ifdef HAVE_MPI
+/**
+ * \brief An overlapping schwarz operator.
+ *
+ * This operator represents a parallel matrix product using
+ * sequential data structures together with a parallel index set
+ * describing an overlapping domain decomposition and the communication.
+ * \tparam M The type of the sequential matrix to use,
+ * e.g. BCRSMatrix or another matrix type fulfilling the
+ * matrix interface of ISTL.
+ * \tparam X The type of the sequential vector to use for the left hand side,
+ * e.g. BlockVector or another type fulfilling the ISTL
+ * vector interface.
+ * \tparam Y The type of the sequential vector to use for the right hand side,
+ * e..g. BlockVector or another type fulfilling the ISTL
+ * vector interface.
+ * \tparam C The type of the communication object.
+ * This must either be OwnerOverlapCopyCommunication or a type
+ * implementing the same interface.
+ *
+ * This is a modified version which does summation of the vector after evaluation.
+ */
+template<class M, class X, class Y, class C>
+class OverlappingSchwarzOperator : public Dune::AssembledLinearOperator<M,X,Y>
+{
+public:
+  //! \brief The type of the matrix we operate on.
+  //!
+  //! E.g. BCRSMatrix or another matrix type fulfilling the
+  //! matrix interface of ISTL
+  typedef M matrix_type;
+  //! \brief The type of the domain.
+  //!
+  //! E.g. BlockVector or another type fulfilling the ISTL
+  //! vector interface.
+  typedef X domain_type;
+  //! \brief The type of the range.
+  //!
+  //! E.g. BlockVector or another type fulfilling the ISTL
+  //! vector interface.
+  typedef Y range_type;
+  //! \brief The field type of the range
+  typedef typename X::field_type field_type;
+  //! \brief The type of the communication object.
+  //!
+  //! This must either be OwnerOverlapCopyCommunication or a type
+  //! implementing the same interface.
+  typedef C communication_type;
+
+  enum {
+    //! \brief The solver category.
+    category=Dune::SolverCategory::overlapping
+  };
+
+  /**
+   * @brief constructor: just store a reference to a matrix.
+   *
+   * @param A The assembled matrix.
+   * @param com The communication object for syncing overlap and copy
+   * data points. (E.~g. OwnerOverlapCopyCommunication )
+   */
+  OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
+    : _A_(A), communication(com)
+  {}
+
+  //! apply operator to x:  \f$ y = A(x) \f$
+  virtual void apply (const X& x, Y& y) const
+  {
+//    Y y2(y.size());
+    _A_.umv(x,y);
+    communication.addOwnerOverlapToAll(y,y);
+  }
+
+  //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
+  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+  {
+    Y y2(y.size());
+    _A_.usmv(alpha, x, y2);
+    y = 0;
+    communication.addOwnerOverlapToAll(y2,y);
+  }
+
+  //! get the sequential assembled linear operator.
+  virtual const matrix_type& getmat () const
+  {
+    return _A_;
+  }
+
+private:
+  const matrix_type& _A_;
+  const communication_type& communication;
+};
+
+typedef OverlappingSchwarzOperator<ISTL::Mat,ISTL::Vec,ISTL::Vec,Dune::OwnerOverlapCopyCommunication<int,int>> ParMatrixAdapter;
+#endif
+
+}
+
+
+/*!
+  \brief Class for ISTL solver parameters.
+  \details It contains information about solver method, preconditioner
+  and convergence criteria.
+*/
+
+class ISTLSolParams
+{
+public:
+  //! \brief Constructor.
+  //! \param spar Linear solver parameters to use
+  //! \param padm Process administrator to use
+  ISTLSolParams(const LinSolParams& spar, const ProcessAdm& padm) :
+    solParams(spar), adm(padm)
+  {}
+
+
+  //! \brief Setup solver and preconditioner.
+  //! \param A Matrix to use
+  //! \return tuple with (solver, preconditioner, op)
+  std::tuple<std::unique_ptr<ISTL::InverseOperator>,
+             std::unique_ptr<ISTL::Preconditioner>,
+             std::unique_ptr<ISTL::Operator>>
+    setupPC(ISTL::Mat& A);
+
+  //! \brief Obtain linear solver parameters.
+  const LinSolParams& get() const { return solParams; }
+
+protected:
+  //! \brief Internal helper function for setting up a preconditioner.
+  //! \param A Matrix to construct preconditioner for
+  //! \param op Matrix adaptor to use (used with AMG)
+  //! \param block Block to read settings from
+  //! \param solver Solver object to instance. nullptr to do nothing.
+  ISTL::Preconditioner* setupPCInternal(ISTL::Mat& A,
+                                        ISTL::Operator& op,
+                                        size_t block,
+                                        std::unique_ptr<ISTL::InverseOperator>* solver);
+
+  const LinSolParams& solParams; //!< Reference to linear solver parameters.
+  const ProcessAdm& adm;      //!< Reference to process administrator.
+};
+
+#endif
diff --git a/src/LinAlg/ISTLSupport.h b/src/LinAlg/ISTLSupport.h
new file mode 100644
index 00000000..c9fafaad
--- /dev/null
+++ b/src/LinAlg/ISTLSupport.h
@@ -0,0 +1,76 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ISTLSupport.h
+//!
+//! \date May 12 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief IFEM ISTL support
+//!
+//==============================================================================
+
+#ifndef _ISTL_SUPPORT_H_
+#define _ISTL_SUPPORT_H_
+
+#include <vector>
+
+#ifdef HAS_ISTL
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/overlappingschwarz.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/schwarz.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#ifdef HAVE_SUPERLU
+#include <dune/istl/superlu.hh>
+#endif
+#ifdef HAVE_UMFPACK
+#include <dune/istl/umfpack.hh>
+#endif
+
+namespace ISTL
+{
+  typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>> Mat; //!< A sparse system matrix
+  typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vec;  //!< A vector
+  typedef Dune::MatrixAdapter<Mat,Vec,Vec> Operator;      //!< A serial matrix operator
+  typedef Dune::InverseOperator<Vec, Vec> InverseOperator;     //!< Linear system inversion abstraction
+  typedef Dune::Preconditioner<Vec,Vec> Preconditioner;        //!< Preconditioner abstraction
+
+  /*! \brief Wrapper template to avoid memory leaks */
+  template<template<class M> class Pre>
+  class IOp2Pre : public Dune::InverseOperator2Preconditioner<Pre<ISTL::Mat>, Dune::SolverCategory::sequential> {
+    typedef Dune::InverseOperator2Preconditioner<Pre<ISTL::Mat>, Dune::SolverCategory::sequential> SolverType;
+  public:
+    IOp2Pre(Pre<ISTL::Mat>* iop) : SolverType(*iop)
+    {
+      m_op.reset(iop);
+    }
+  protected:
+    std::unique_ptr<Pre<ISTL::Mat>> m_op;
+  };
+
+#if defined(HAVE_UMFPACK)
+  typedef Dune::UMFPack<ISTL::Mat> LUType;
+  typedef IOp2Pre<Dune::UMFPack> LU;
+#elif defined(HAVE_SUPERLU)
+  typedef Dune::SuperLU<ISTL::Mat> LUType;
+  typedef IOp2Pre<Dune::SuperLU> LU;
+#endif
+
+  // Preconditioner types
+  typedef Dune::SeqILU0<Mat, Vec, Vec> ILU0; //!< Sequential ILU(0)
+  typedef Dune::SeqILUn<Mat, Vec, Vec> ILUn; //!< Sequential ILU(n)
+  typedef  Dune::SeqSOR<Mat, Vec, Vec>  SOR; //!< Sequential SOR
+  typedef Dune::SeqSSOR<Mat, Vec,Vec> SSOR; //!< Sequential SSOR
+  typedef  Dune::SeqJac<Mat, Vec, Vec>    GJ; //!< Sequential Gauss-Jacobi
+  typedef   Dune::SeqGS<Mat, Vec, Vec>     GS; //!< Sequential Gauss-Seidel
+  typedef Dune::SeqOverlappingSchwarz<Mat, Vec,
+                                 Dune::AdditiveSchwarzMode,
+                                 LUType> ASMLU; //!< Sequential additive Schwarz w/ LU subdomain solvers
+  typedef Dune::SeqOverlappingSchwarz<Mat, Vec> ASM; //!< Sequential additive Schwarz w/ ILU subdomain solvers
+}
+#endif
+
+#endif
diff --git a/src/LinAlg/SystemMatrix.C b/src/LinAlg/SystemMatrix.C
index 68b14ec0..8e94c8fa 100644
--- a/src/LinAlg/SystemMatrix.C
+++ b/src/LinAlg/SystemMatrix.C
@@ -13,12 +13,15 @@
 
 #include "SystemMatrix.h"
 #include "DenseMatrix.h"
+#ifdef HAS_ISTL
+#include "ISTLMatrix.h"
+#endif
 #include "SPRMatrix.h"
 #include "SparseMatrix.h"
 #ifdef HAS_PETSC
-#include "LinSolParams.h"
 #include "PETScMatrix.h"
 #endif
+#include "LinSolParams.h"
 
 
 SystemVector* SystemVector::create (const ProcessAdm& adm, Type vectorType)
@@ -26,6 +29,9 @@ SystemVector* SystemVector::create (const ProcessAdm& adm, Type vectorType)
   switch (vectorType)
     {
     case STD   : return new StdVector();
+#ifdef HAS_ISTL
+    case ISTL  : return new ISTLVector(adm);
+#endif
 #ifdef HAS_PETSC
     case PETSC : return new PETScVector(adm);
 #endif
@@ -73,6 +79,10 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
   if (matrixType == PETSC)
     return new PETScMatrix(padm,spar,ltype);
 #endif
+#ifdef HAS_ISTL
+  if (matrixType == ISTL)
+    return new ISTLMatrix(padm,spar,ltype);
+#endif
 
   return SystemMatrix::create(padm,matrixType,ltype);
 }
@@ -88,7 +98,16 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
               << std::endl;
     exit(1);
   }
-#else
+#endif
+#ifndef HAS_ISTL
+  if (matrixType == ISTL) {
+    std::cerr <<"SystemMatrix::create: ISTL not compiled in, bailing out..."
+              << std::endl;
+    exit(1);
+  }
+#endif
+
+#if defined(HAS_PETSC) || defined(HAS_ISTL)
   // Use default settings when no parameters are provided by user
   static LinSolParams defaultPar;
 #endif
@@ -99,6 +118,9 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
     case SPR   : return new SPRMatrix();
     case SPARSE: return new SparseMatrix(SparseMatrix::SUPERLU,num_thread_SLU);
     case SAMG  : return new SparseMatrix(SparseMatrix::S_A_M_G);
+#ifdef HAS_ISTL
+    case ISTL  : return new ISTLMatrix(padm,defaultPar,ltype);
+#endif
 #ifdef HAS_PETSC
     case PETSC :      return new PETScMatrix(padm,defaultPar,ltype);
 #endif
diff --git a/src/LinAlg/SystemMatrix.h b/src/LinAlg/SystemMatrix.h
index 4926198d..f085df5e 100644
--- a/src/LinAlg/SystemMatrix.h
+++ b/src/LinAlg/SystemMatrix.h
@@ -30,7 +30,7 @@ class SystemVector
 {
 public:
   //! \brief The available system vector formats.
-  enum Type { STD = 0, PETSC = 1 };
+  enum Type { STD = 0, PETSC = 1, ISTL = 2 };
 
   //! \brief Static method creating a vector of the given type.
   static SystemVector* create(const ProcessAdm& padm, Type vectorType = STD);
@@ -203,7 +203,7 @@ class SystemMatrix
 public:
   //! \brief The available system matrix formats.
   enum Type { DENSE = 0, SPR = 1, SPARSE = 2, SAMG = 3,
-              PETSC = 4 };
+              PETSC = 4, ISTL = 5 };
 
   //! \brief Static method creating a matrix of the given type.
   static SystemMatrix* create(const ProcessAdm& padm, Type matrixType,
diff --git a/src/LinAlg/Test/MPI/TestISTLMatrix.C b/src/LinAlg/Test/MPI/TestISTLMatrix.C
new file mode 100644
index 00000000..46bffab1
--- /dev/null
+++ b/src/LinAlg/Test/MPI/TestISTLMatrix.C
@@ -0,0 +1,110 @@
+//==============================================================================
+//!
+//! \file TestISTLMatrix.C
+//!
+//! \date Apr 27 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Unit tests for parallel ISTL matrices
+//!
+//==============================================================================
+
+#include "AlgEqSystem.h"
+#include "ISTLMatrix.h"
+#include "ISTLSolParams.h"
+#include "SIM2D.h"
+#include "ASMs2D.h"
+#include "IntegrandBase.h"
+#include "SAM.h"
+#include "IFEM.h"
+
+#include "gtest/gtest.h"
+
+#include <fstream>
+
+
+typedef std::vector<int> IntVec;
+
+static IntVec readIntVector(const std::string& file)
+{
+  std::vector<int> result;
+  std::ifstream f(file);
+  size_t size;
+  f >> size;
+  result.resize(size);
+  for (size_t j=0;j<size;++j)
+    f >> result[j];
+
+  return result;
+}
+
+
+class DummyIntegrandForDummies : public IntegrandBase {
+public:
+  DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
+};
+
+
+class InspectMatrixSIM : public SIM2D {
+public:
+  InspectMatrixSIM(unsigned char n1 = 2, bool check = false) :
+    SIM2D(n1, check) { myProblem = new DummyIntegrandForDummies; }
+
+  InspectMatrixSIM(const std::vector<unsigned char>& n, bool check = false) :
+    SIM2D(n, check) { myProblem = new DummyIntegrandForDummies; }
+
+  SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
+  SystemVector* getVector() { return myEqSys->getVector(0); }
+};
+
+
+TEST(TestISTLMatrix, AssembleMPI)
+{
+  InspectMatrixSIM sim(1);
+  sim.read("src/LinAlg/Test/refdata/petsc_test.xinp");
+  sim.opt.solver = SystemMatrix::ISTL;
+  sim.preprocess();
+  sim.initSystem(SystemMatrix::ISTL);
+
+  Matrix stencil(4,4);
+  stencil(1,1) = stencil(2,2) = stencil(3,3) = stencil(4,4) = 1.0;
+
+  for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
+    sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
+
+  sim.getMatrix()->beginAssembly();
+  sim.getMatrix()->endAssembly();
+
+  // now inspect the matrix
+  const ProcessAdm& adm = sim.getProcessAdm();
+  ISTL::Mat& mat = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
+  ISTL::Vec b(mat.N()), b2(mat.N());
+
+  try {
+    Dune::OwnerOverlapCopyCommunication<int,int> comm(*adm.getCommunicator());
+    comm.indexSet().beginResize();
+    typedef Dune::ParallelLocalIndex<Dune::OwnerOverlapCopyAttributeSet::AttributeSet> LI;
+    for (size_t i = 0; i < adm.dd.getMLGEQ().size(); ++i) {
+      int gid = adm.dd.getGlobalEq(i+1);
+      comm.indexSet().add(gid-1, LI(i, gid >= adm.dd.getMinEq() ?
+                                       Dune::OwnerOverlapCopyAttributeSet::owner :
+                                       Dune::OwnerOverlapCopyAttributeSet::overlap));
+    }
+    comm.indexSet().endResize();
+    comm.remoteIndices().setIncludeSelf(true);
+    comm.remoteIndices().template rebuild<false>();
+
+    ISTL::ParMatrixAdapter op(mat, comm);
+
+    b = 1.0;
+    op.apply(b, b2);
+  } catch (Dune::ISTLError e) {
+    std::cerr << e << std::endl;
+    ASSERT_TRUE(false);
+  }
+
+  IntVec v = readIntVector("src/LinAlg/Test/refdata/petsc_matrix_diagonal.ref");
+  for (size_t i = 1; i <= adm.dd.getMLGEQ().size(); ++i)
+    ASSERT_FLOAT_EQ(v[adm.dd.getGlobalEq(i)-1], b2[i-1]);
+}
diff --git a/src/LinAlg/Test/TestISTLMatrix.C b/src/LinAlg/Test/TestISTLMatrix.C
new file mode 100644
index 00000000..879ed923
--- /dev/null
+++ b/src/LinAlg/Test/TestISTLMatrix.C
@@ -0,0 +1,205 @@
+//==============================================================================
+//!
+//! \file TestISTLMatrix.C
+//!
+//! \date Apr 27 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Unit tests for ISTL matrices
+//!
+//==============================================================================
+
+#include "AlgEqSystem.h"
+#include "ISTLMatrix.h"
+#include "SIM2D.h"
+#include "ASMs2D.h"
+#include "IntegrandBase.h"
+#include "SAM.h"
+#include <dune/istl/io.hh>
+
+#include "gtest/gtest.h"
+
+#include <fstream>
+
+
+typedef std::vector<int> IntVec;
+
+static IntVec readIntVector(const std::string& file)
+{
+  std::vector<int> result;
+  std::ifstream f(file);
+  size_t size;
+  f >> size;
+  result.resize(size);
+  for (size_t j=0;j<size;++j)
+    f >> result[j];
+
+  return result;
+}
+
+
+class DummyIntegrandForDummies : public IntegrandBase {
+public:
+  DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
+};
+
+
+class InspectMatrixSIM : public SIM2D {
+public:
+  InspectMatrixSIM(unsigned char n1 = 2, bool check = false) :
+    SIM2D(n1, check) { myProblem = new DummyIntegrandForDummies; }
+
+  InspectMatrixSIM(const std::vector<unsigned char>& n, bool check = false) :
+    SIM2D(n, check) { myProblem = new DummyIntegrandForDummies; }
+
+  SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
+};
+
+
+class InspectBlockPreconditioner : public ISTL::BlockPreconditioner {
+public:
+  InspectBlockPreconditioner(const ISTL::Mat& A,
+                             const DomainDecomposition& dd) :
+    ISTL::BlockPreconditioner(A, dd, "upper") {}
+
+  void extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
+                    const std::set<int>& eqs_row,
+                    const std::set<int>& eqs_col,
+                    const std::map<int,int>& eqs_col_g2l)
+  {
+    ISTL::BlockPreconditioner::extractBlock(B, A, eqs_row, eqs_col, eqs_col_g2l);
+  }
+};
+
+
+TEST(TestISTLMatrix, Assemble)
+{
+  InspectMatrixSIM sim(1);
+  sim.read("src/LinAlg/Test/refdata/petsc_test.xinp");
+  sim.opt.solver = SystemMatrix::ISTL;
+  sim.preprocess();
+  sim.initSystem(SystemMatrix::ISTL);
+
+  Matrix stencil(4,4);
+  stencil(1,1) = stencil(2,2) = stencil(3,3) = stencil(4,4) = 1.0;
+
+  for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
+    sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
+
+  sim.getMatrix()->beginAssembly();
+  sim.getMatrix()->endAssembly();
+
+  // now inspect the matrix
+  ISTL::Mat& mat = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
+  IntVec v = readIntVector("src/LinAlg/Test/refdata/petsc_matrix_diagonal.ref");
+
+  ASSERT_EQ(mat.N(), v.size());
+
+  for (size_t i = 0; i < v.size(); ++i)
+    ASSERT_FLOAT_EQ(double(v[i]), mat[i][i]);
+}
+
+
+TEST(TestISTLMatrix, AssembleBasisBlocks)
+{
+  InspectMatrixSIM sim({1,1});
+  sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
+  sim.opt.solver = SystemMatrix::ISTL;
+  sim.preprocess();
+  sim.initSystem(SystemMatrix::ISTL);
+
+  Matrix stencil(13,13);
+  for (size_t i = 1; i<= 9; ++i)
+    stencil(i,i) = 1.0;
+  for (size_t i = 10; i<= 13; ++i)
+    stencil(i,i) = 2.0;
+
+  for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
+    sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
+
+  sim.getMatrix()->beginAssembly();
+  sim.getMatrix()->endAssembly();
+
+  const ProcessAdm& adm = sim.getProcessAdm();
+  ISTL::Mat& A = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
+  InspectBlockPreconditioner block(A, adm.dd);
+
+  // now inspect the matrix blocks
+  std::vector<ISTL::Mat> mat(5);
+  block.extractBlock(mat[0], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
+  block.extractBlock(mat[1], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
+  block.extractBlock(mat[2], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
+  block.extractBlock(mat[3], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
+  mat[4] = block.getBlock(1); // S = D
+
+  for (size_t b = 0; b < mat.size(); ++b) {
+    if (b == 0 || b >= 3) { // diagonal blocks
+      std::stringstream str;
+      str << "src/LinAlg/Test/refdata/petsc_matrix_diagonal_basis_block" << std::min(b/2 + 1, 2ul) << ".ref";
+      IntVec v = readIntVector(str.str());
+      ASSERT_EQ(v.size(), mat[b].N());
+      for (size_t i = 0; i < v.size(); ++i)
+        ASSERT_FLOAT_EQ(double(v[i]), mat[b][i][i]);
+    }
+
+    // check that no values outside the diagonal are != 0
+    for (auto r = mat[b].begin(); r != mat[b].end(); ++r) {
+      for (auto it = r->begin(); it != r->end(); ++it)
+        if (r.index() != it.index())
+          ASSERT_FLOAT_EQ(*it, 0.0);
+    }
+  }
+}
+
+
+TEST(TestISTLMatrix, AssembleComponentBlocks)
+{
+  InspectMatrixSIM sim(2);
+  sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_components.xinp");
+  sim.opt.solver = SystemMatrix::ISTL;
+  sim.preprocess();
+  sim.initSystem(SystemMatrix::ISTL);
+
+  Matrix stencil(4*2,4*2);
+  for (size_t i = 1; i<= 4; ++i) {
+    stencil(2*i-1,2*i-1) = 1.0;
+    stencil(2*i,2*i) = 2.0;
+  }
+
+  for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
+    sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
+
+  sim.getMatrix()->beginAssembly();
+  sim.getMatrix()->endAssembly();
+
+  const ProcessAdm& adm = sim.getProcessAdm();
+  ISTL::Mat& A = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
+  InspectBlockPreconditioner block(A, adm.dd);
+
+  // now inspect the matrix blocks
+  std::vector<ISTL::Mat> mat(5);
+  block.extractBlock(mat[0], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
+  block.extractBlock(mat[1], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
+  block.extractBlock(mat[2], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
+  block.extractBlock(mat[3], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
+  mat[4] = block.getBlock(1); // S = D
+
+  for (size_t b = 0; b < mat.size(); ++b) {
+    if (b == 0 || b >= 3) { // diagonal blocks
+      std::stringstream str;
+      str << "src/LinAlg/Test/refdata/petsc_matrix_diagonal_components_block" << std::min(b/2 + 1, 2ul) << ".ref";
+      IntVec v = readIntVector(str.str());
+      ASSERT_EQ(v.size(), mat[b].N());
+      for (size_t i = 0; i < v.size(); ++i)
+        ASSERT_FLOAT_EQ(double(v[i]), mat[b][i][i]);
+    }
+
+    // check that no values outside the diagonal are != 0
+    for (auto r = mat[b].begin(); r != mat[b].end(); ++r) {
+      for (auto it = r->begin(); it != r->end(); ++it)
+        if (r.index() != it.index())
+          ASSERT_FLOAT_EQ(*it, 0.0);
+    }
+  }
+}
diff --git a/src/LinAlg/Test/TestISTLPETScMatrix.C b/src/LinAlg/Test/TestISTLPETScMatrix.C
new file mode 100644
index 00000000..e301b61d
--- /dev/null
+++ b/src/LinAlg/Test/TestISTLPETScMatrix.C
@@ -0,0 +1,82 @@
+//==============================================================================
+//!
+//! \file TestISTLPETScMatrix.C
+//!
+//! \date Apr 27 2016
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Unit tests comparing ISTL and PETSc matrices
+//!
+//==============================================================================
+
+#include "AlgEqSystem.h"
+#include "ISTLMatrix.h"
+#include "PETScMatrix.h"
+#include "SIM2D.h"
+#include "ASMs2D.h"
+#include "IntegrandBase.h"
+#include "SAM.h"
+#include <dune/istl/io.hh>
+
+#include "gtest/gtest.h"
+
+#include <fstream>
+
+
+class DummyIntegrandForDummies : public IntegrandBase {
+public:
+  DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
+};
+
+
+class InspectMatrixSIM : public SIM2D {
+public:
+  InspectMatrixSIM() :
+    SIM2D({1,1}, false) { myProblem = new DummyIntegrandForDummies; }
+
+  SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
+};
+
+
+
+TEST(TestISTLPETScMatrix, SchurComplement)
+{
+  Matrix stencil(13,13);
+  for (size_t i = 1; i<= 13; ++i)
+    for (size_t j = 1; j <= 13; ++j)
+      stencil(i,j) = 1.0;
+
+  std::array<InspectMatrixSIM,2> sim;
+  for (size_t i = 0; i < 2; ++i) {
+    sim[i].read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
+    sim[i].opt.solver = i == 0 ? SystemMatrix::PETSC : SystemMatrix::ISTL;
+    sim[i].preprocess();
+    sim[i].initSystem(i == 0 ? SystemMatrix::PETSC : SystemMatrix::ISTL);
+
+    for (int iel = 1; iel <= sim[i].getSAM()->getNoElms(); ++iel)
+      sim[i].getMatrix()->assemble(stencil, *sim[i].getSAM(), iel);
+
+    sim[i].getMatrix()->beginAssembly();
+    sim[i].getMatrix()->endAssembly();
+  }
+
+  const ProcessAdm& adm = sim[1].getProcessAdm();
+  ISTL::Mat& A = static_cast<ISTLMatrix*>(sim[1].getMatrix())->getMatrix();
+  ISTL::BlockPreconditioner block(A, adm.dd, "upper");
+
+  ISTL::Mat& S = block.getBlock(1);
+  PETScSolParams params(LinSolParams(), adm);
+  params.setupSchurComplement(static_cast<PETScMatrix*>(sim[0].getMatrix())->getBlockMatrices());
+
+  // check that matrices are the same
+  for (size_t r = 0; r < S.N(); ++r) {
+    const PetscInt* cols;
+    PetscInt ncols;
+    const PetscScalar* vals;
+    MatGetRow(params.getSchurComplement(), r, &ncols, &cols, &vals);
+    for (PetscInt i = 0; i < ncols; ++i)
+      ASSERT_FLOAT_EQ(vals[i], S[r][cols[i]]);
+    MatRestoreRow(params.getSchurComplement(), r, &ncols, &cols, &vals);
+  }
+}
diff --git a/src/SIM/SIMoptions.C b/src/SIM/SIMoptions.C
index bd6ad97b..166d43b2 100644
--- a/src/SIM/SIMoptions.C
+++ b/src/SIM/SIMoptions.C
@@ -68,6 +68,8 @@ void SIMoptions::setLinearSolver (const std::string& eqsolver)
     solver = SystemMatrix::SAMG;
   else if (eqsolver == "petsc")
     solver = SystemMatrix::PETSC;
+  else if (eqsolver == "istl")
+    solver = SystemMatrix::ISTL;
 }
 
 
@@ -223,6 +225,8 @@ bool SIMoptions::parseOldOptions (int argc, char** argv, int& i)
     solver = SystemMatrix::SAMG;
   else if (!strcmp(argv[i],"-petsc"))
     solver = SystemMatrix::PETSC;
+  else if (!strcmp(argv[i],"-istl"))
+    solver = SystemMatrix::ISTL;
   else if (!strncmp(argv[i],"-lag",4))
     discretization = ASM::Lagrange;
   else if (!strncmp(argv[i],"-spec",5))