diff --git a/CMakeLists.txt b/CMakeLists.txt
index 93e29a96ff0090554ba3021094a2b61e3a660dc8..dc0f693980d6771d005c961cfae8f023ae81d174 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,6 +2,7 @@ cmake_minimum_required (VERSION 3.10)
 
 # CMake utils
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/packages/cmake-modules")
 
 # Forbids in-source builds
 include(CheckNotInSources)
@@ -113,13 +114,12 @@ if (NOT PUGS_ENABLE_MPI MATCHES "^(AUTO|ON|OFF)$")
   message(FATAL_ERROR "PUGS_ENABLE_MPI='${PUGS_ENABLE_MPI}'. Must be set to one of AUTO, ON or OFF")
 endif()
 
-# checks for MPI
+# Check for MPI
 if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$")
   set(MPI_DETERMINE_LIBRARY_VERSION TRUE)
   find_package(MPI)
 endif()
 
-
 #------------------------------------------------------
 # Search for ParMETIS
 
@@ -137,6 +137,24 @@ if(${MPI_FOUND})
   endif()
 endif()
 
+#------------------------------------------------------
+# Check for PETSc
+# defaults use PETSc
+set(PUGS_ENABLE_PETSC AUTO CACHE STRING
+  "Choose one of: AUTO ON OFF")
+
+if (PUGS_ENABLE_PETSC MATCHES "^(AUTO|ON)$")
+  find_package(PETSc)
+  set(PUGS_HAS_PETSC ${PETSC_FOUND})
+else()
+  unset(PUGS_HAS_PETSC)
+endif()
+
+if (${PETSC_FOUND})
+  include_directories(SYSTEM ${PETSC_INCLUDES})
+endif()
+
+
 # -----------------------------------------------------
 
 if (${MPI_FOUND})
@@ -414,15 +432,18 @@ add_executable(
 target_link_libraries(
   pugs
   PugsMesh
+  PugsAlgebra
   PugsUtils
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
   PugsLanguageAlgorithms
   PugsMesh
+  PugsAlgebra
   PugsUtils
   PugsLanguageUtils
   kokkos
+  ${PETSC_LIBRARIES}
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
   ${KOKKOS_CXX_FLAGS}
@@ -450,37 +471,51 @@ endif()
 
 message("")
 message("====== pugs build options ======")
-message(STATUS "version: ${PUGS_VERSION}")
-message(STATUS "build type: ${CMAKE_BUILD_TYPE}")
-message(STATUS "compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
-message(STATUS "kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}")
+message(" version: ${PUGS_VERSION}")
+message(" build type: ${CMAKE_BUILD_TYPE}")
+message(" compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
+message(" kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}")
 
 if (MPI_FOUND)
-  message(STATUS "MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}")
+  message(" MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}")
 else()
   if(NOT PARMETIS_LIBRARIES)
-    message(STATUS "MPI: deactivated: ParMETIS cannot be found!")
+    message(" MPI: deactivated: ParMETIS cannot be found!")
+  else()
+    if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$")
+      message(" MPI: not found!")
+    else()
+      message(" MPI: explicitly deactivated!")
+    endif()
+  endif()
+endif()
+
+if (PETSC_FOUND)
+  message(" PETSc: ${PETSC_VERSION}")
+else()
+  if (PUGS_ENABLE_PETSC MATCHES "^(AUTO|ON)$")
+    message(" PETSc: not found!")
   else()
-    message(STATUS "MPI: not found!")
+      message(" PETSc: explicitly deactivated!")
   endif()
 endif()
 
 if(CLANG_FORMAT)
-  message(STATUS "clang-format: ${CLANG_FORMAT}")
+  message(" clang-format: ${CLANG_FORMAT}")
 else()
-  message(STATUS "clang-format: not found!")
+  message(" clang-format: not found!")
 endif()
 
 if(CLAZY_STANDALONE)
-  message(STATUS "clazy-standalone: ${CLAZY_STANDALONE}")
+  message(" clazy-standalone: ${CLAZY_STANDALONE}")
 else()
-  message(STATUS "clazy-standalone: no found!")
+  message(" clazy-standalone: no found!")
 endif()
 
 if (DOXYGEN_FOUND)
-  message(STATUS "doxygen: ${DOXYGEN_EXECUTABLE}")
+  message(" doxygen: ${DOXYGEN_EXECUTABLE}")
 else()
-  message(STATUS "doxygen: no found!")
+  message(" doxygen: no found!")
 endif()
 
 message("================================")
diff --git a/packages/cmake-modules/.gitrepo b/packages/cmake-modules/.gitrepo
new file mode 100644
index 0000000000000000000000000000000000000000..40d20d1839ad63c531af16567f906d2042f3e2f7
--- /dev/null
+++ b/packages/cmake-modules/.gitrepo
@@ -0,0 +1,12 @@
+; DO NOT EDIT (unless you know what you are doing)
+;
+; This subdirectory is a git "subrepo", and this file is maintained by the
+; git-subrepo command. See https://github.com/git-commands/git-subrepo#readme
+;
+[subrepo]
+	remote = https://github.com/jedbrown/cmake-modules.git
+	branch = master
+	commit = 91f96174a8b3f65e19519fa592b1571391c0e3d0
+	parent = 14615ef666ef3a26dcccd36f0a357b493e6cfbe5
+	method = merge
+	cmdver = 0.4.1
diff --git a/packages/cmake-modules/CorrectWindowsPaths.cmake b/packages/cmake-modules/CorrectWindowsPaths.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..09bcdd67dcd04fd001d2b7acbd904b5014ebe42b
--- /dev/null
+++ b/packages/cmake-modules/CorrectWindowsPaths.cmake
@@ -0,0 +1,14 @@
+# CorrectWindowsPaths - this module defines one macro
+#
+# CONVERT_CYGWIN_PATH( PATH )
+#  This uses the command cygpath (provided by cygwin) to convert
+#  unix-style paths into paths useable by cmake on windows
+
+macro (CONVERT_CYGWIN_PATH _path)
+  if (WIN32)
+    EXECUTE_PROCESS(COMMAND cygpath.exe -m ${${_path}}
+      OUTPUT_VARIABLE ${_path})
+    string (STRIP ${${_path}} ${_path})
+  endif (WIN32)
+endmacro (CONVERT_CYGWIN_PATH)
+
diff --git a/packages/cmake-modules/FindFFTW.cmake b/packages/cmake-modules/FindFFTW.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..00c3401c11d7b32871290d4a4d2e9e9717411962
--- /dev/null
+++ b/packages/cmake-modules/FindFFTW.cmake
@@ -0,0 +1,22 @@
+# - Find FFTW
+# Find the native FFTW includes and library
+#
+#  FFTW_INCLUDES    - where to find fftw3.h
+#  FFTW_LIBRARIES   - List of libraries when using FFTW.
+#  FFTW_FOUND       - True if FFTW found.
+
+if (FFTW_INCLUDES)
+  # Already in cache, be silent
+  set (FFTW_FIND_QUIETLY TRUE)
+endif (FFTW_INCLUDES)
+
+find_path (FFTW_INCLUDES fftw3.h)
+
+find_library (FFTW_LIBRARIES NAMES fftw3)
+
+# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (FFTW DEFAULT_MSG FFTW_LIBRARIES FFTW_INCLUDES)
+
+mark_as_advanced (FFTW_LIBRARIES FFTW_INCLUDES)
diff --git a/packages/cmake-modules/FindGSL.cmake b/packages/cmake-modules/FindGSL.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..64a578fa6904af76e113c7b5d2b1e2ae74e630d7
--- /dev/null
+++ b/packages/cmake-modules/FindGSL.cmake
@@ -0,0 +1,32 @@
+# - Find GSL
+# Find the native GSL includes and library
+#
+#  GSL_INCLUDES    - where to find gsl/gsl_*.h, etc.
+#  GSL_LIBRARIES   - List of libraries when using GSL.
+#  GSL_FOUND       - True if GSL found.
+
+
+if (GSL_INCLUDES)
+  # Already in cache, be silent
+  set (GSL_FIND_QUIETLY TRUE)
+endif (GSL_INCLUDES)
+
+find_path (GSL_INCLUDES gsl/gsl_math.h)
+
+find_library (GSL_LIB NAMES gsl)
+
+set (GSL_CBLAS_LIB "" CACHE FILEPATH "If your program fails to link
+(usually because GSL is not automatically linking a CBLAS and no other
+component of your project provides a CBLAS) then you may need to point
+this variable to a valid CBLAS.  Usually GSL is distributed with
+libgslcblas.{a,so} (next to GSL_LIB) which you may use if an optimized
+CBLAS is unavailable.")
+
+set (GSL_LIBRARIES "${GSL_LIB}" "${GSL_CBLAS_LIB}")
+
+# handle the QUIETLY and REQUIRED arguments and set GSL_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (GSL DEFAULT_MSG GSL_LIBRARIES GSL_INCLUDES)
+
+mark_as_advanced (GSL_LIB GSL_CBLAS_LIB GSL_INCLUDES)
diff --git a/packages/cmake-modules/FindGit.cmake b/packages/cmake-modules/FindGit.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..05c9bf0be638ad2049c8fbf6a260dd2f43270e8f
--- /dev/null
+++ b/packages/cmake-modules/FindGit.cmake
@@ -0,0 +1,34 @@
+SET(Git_FOUND FALSE)
+
+FIND_PROGRAM(Git_EXECUTABLE git
+  DOC "git command line client")
+MARK_AS_ADVANCED(Git_EXECUTABLE)
+
+IF(Git_EXECUTABLE)
+  SET(Git_FOUND TRUE)
+  MACRO(Git_WC_INFO dir prefix)
+    EXECUTE_PROCESS(COMMAND ${Git_EXECUTABLE} rev-list -n 1 HEAD
+       WORKING_DIRECTORY ${dir}
+        ERROR_VARIABLE Git_error
+       OUTPUT_VARIABLE ${prefix}_WC_REVISION_HASH
+        OUTPUT_STRIP_TRAILING_WHITESPACE)
+    if(NOT ${Git_error} EQUAL 0)
+      MESSAGE(SEND_ERROR "Command \"${Git_EXECUTBALE} rev-list -n 1 HEAD\" in directory ${dir} failed with output:\n${Git_error}")
+    ELSE(NOT ${Git_error} EQUAL 0)
+      EXECUTE_PROCESS(COMMAND ${Git_EXECUTABLE} name-rev ${${prefix}_WC_REVISION_HASH}
+         WORKING_DIRECTORY ${dir}
+         OUTPUT_VARIABLE ${prefix}_WC_REVISION_NAME
+          OUTPUT_STRIP_TRAILING_WHITESPACE)
+    ENDIF(NOT ${Git_error} EQUAL 0)
+  ENDMACRO(Git_WC_INFO)
+ENDIF(Git_EXECUTABLE)
+
+IF(NOT Git_FOUND)
+  IF(NOT Git_FIND_QUIETLY)
+    MESSAGE(STATUS "Git was not found")
+  ELSE(NOT Git_FIND_QUIETLY)
+    if(Git_FIND_REQUIRED)
+      MESSAGE(FATAL_ERROR "Git was not found")
+    ENDIF(Git_FIND_REQUIRED)
+  ENDIF(NOT Git_FIND_QUIETLY)
+ENDIF(NOT Git_FOUND)
diff --git a/packages/cmake-modules/FindITAPS.cmake b/packages/cmake-modules/FindITAPS.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..91ba0ca2139ace6c04bf4df4b76e8845c2e44d90
--- /dev/null
+++ b/packages/cmake-modules/FindITAPS.cmake
@@ -0,0 +1,187 @@
+# - Try to find ITAPS
+#
+# This will define
+#
+#  ITAPS_FOUND          - Requested components were found
+#  ITAPS_INCLUDES       - Includes for all available components
+#  ITAPS_LIBRARIES      - Libraries for all available components
+#
+#  ITAPS_MESH_FOUND     - System has iMesh
+#  ITAPS_MESH_INCLUDES  - The iMesh include directory
+#  ITAPS_MESH_LIBRARIES - Link these to use iMesh
+#
+#  ITAPS_GEOM_FOUND     - System has iGeom
+#  ITAPS_GEOM_INCLUDES  - The iGeom include directory
+#  ITAPS_GEOM_LIBRARIES - Link these to use iGeom
+#
+#  ITAPS_REL_FOUND      - System has iRel
+#  ITAPS_REL_INCLUDES   - The iRel include directory
+#  ITAPS_REL_LIBRARIES  - Link these to use iRel
+#
+# Setting this changes the behavior of the search
+#  ITAPS_MESH_DEFS_FILE - path to iMesh-Defs.inc
+#  ITAPS_GEOM_DEFS_FILE - path to iGeom-Defs.inc
+#  ITAPS_REL_DEFS_FILE  - path to iRel-Defs.inc
+#
+# If any of these variables are in your environment, they will be used as hints
+#  IMESH_DIR - directory in which iMesh resides
+#  IGEOM_DIR - directory in which iGeom resides
+#  IREL_DIR  - directory in which iRel resides
+#
+# Redistribution and use is allowed according to the terms of the BSD license.
+# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
+#
+
+include (FindPackageMultipass)
+include (ResolveCompilerPaths)
+include (CheckCSourceRuns)
+include (FindPackageHandleStandardArgs)
+
+find_program (MAKE_EXECUTABLE NAMES make gmake)
+
+macro (ITAPS_PREPARE_COMPONENT component name)
+  find_file (ITAPS_${component}_DEFS_FILE ${name}-Defs.inc
+    HINTS ENV I${component}_DIR
+    PATH_SUFFIXES lib64 lib)
+  # If ITAPS_XXX_DEFS_FILE has changed, the library will be found again
+  find_package_multipass (ITAPS_${component} itaps_${component}_config_current
+    STATES DEFS_FILE
+    DEPENDENTS INCLUDES LIBRARIES EXECUTABLE_RUNS)
+endmacro ()
+
+macro (ITAPS_GET_VARIABLE makefile name var)
+  set (${var} "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+  execute_process (COMMAND ${MAKE_EXECUTABLE} -f ${${makefile}} show VARIABLE=${name}
+    OUTPUT_VARIABLE ${var}
+    RESULT_VARIABLE itaps_return)
+endmacro (ITAPS_GET_VARIABLE)
+
+macro (ITAPS_TEST_RUNS component name includes libraries program runs)
+  # message (STATUS "Starting run test: ${includes} ${libraries} ${runs}")
+  multipass_c_source_runs ("${includes}" "${libraries}" "${program}" ${runs})
+  if (NOT ITAPS_${component}_EXECUTABLE_RUNS)
+    set (ITAPS_${component}_EXECUTABLE_RUNS "${${runs}}" CACHE BOOL
+      "Can the system successfully run an ${name} executable?  This variable can be manually set to \"YES\" to force CMake to accept a given configuration, but this will almost always result in a broken build." FORCE)
+  endif ()
+endmacro (ITAPS_TEST_RUNS)
+
+macro (ITAPS_REQUIRED_LIBS component name includes libraries_all program libraries_required)
+  # message (STATUS "trying program: ${program}")
+  resolve_libraries (_all_libraries "${libraries_all}")
+  list (GET _all_libraries 0 _first_library)
+  itaps_test_runs (${component} ${name} "${includes}" "${_first_library};${itaps_rel_libs}" "${program}" ${name}_works_minimal)
+  if (${name}_works_minimal)
+    set (${libraries_required} "${_first_library}")
+    message (STATUS "${name} executable works when only linking to the interface lib, this probably means you have shared libs.")
+  else ()
+    itaps_test_runs (${component} ${name} "${includes}" "${_all_libraries};${itaps_rel_libs}" "${itaps_mesh_program}" ${name}_works_extra)
+    if (${name}_works_extra)
+      set (${libraries_required} "${_all_libraries}")
+      message (STATUS "${name} executable requires linking to extra libs, this probably means it's statically linked.")
+    else ()
+      message (STATUS "${name} could not be used, maybe the install is broken.")
+    endif ()
+  endif ()
+endmacro ()
+
+macro (ITAPS_HANDLE_COMPONENT component name program)
+  itaps_prepare_component ("${component}" "${name}")
+  if (ITAPS_${component}_DEFS_FILE AND NOT itaps_${component}_config_current)
+    # A temporary makefile to probe this ITAPS components's configuration
+    set (itaps_config_makefile "${PROJECT_BINARY_DIR}/Makefile.${name}")
+    file (WRITE ${itaps_config_makefile}
+      "## This file was autogenerated by FindITAPS.cmake
+include ${ITAPS_${component}_DEFS_FILE}
+show :
+\t-@echo -n \${\${VARIABLE}}")
+    itaps_get_variable (itaps_config_makefile I${component}_INCLUDEDIR   itaps_includedir)
+    itaps_get_variable (itaps_config_makefile I${component}_LIBS         itaps_libs)
+    file (REMOVE ${itaps_config_makefile})
+    find_path (itaps_include_tmp ${name}.h HINTS ${itaps_includedir} NO_DEFAULT_PATH)
+    set (ITAPS_${component}_INCLUDES "${itaps_include_tmp}" CACHE STRING "Include directories for ${name}")
+    set (itaps_include_tmp "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+    itaps_required_libs ("${component}" "${name}" "${ITAPS_${component}_INCLUDES}" "${itaps_libs}" "${program}" itaps_${component}_required_libraries)
+    set (ITAPS_${component}_LIBRARIES "${itaps_${component}_required_libraries}" CACHE STRING "Libraries for ${name}")
+    mark_as_advanced (ITAPS_${component}_EXECUTABLE_RUNS ITAPS_${component}_LIBRARIES)
+  endif()
+  set (ITAPS_${component}_FOUND "${ITAPS_${component}_EXECUTABLE_RUNS}")
+endmacro()
+
+itaps_handle_component (MESH iMesh "
+/* iMesh test program */
+#include <iMesh.h>
+#define CHK(err) if (err) return 1
+int main(int argc,char *argv[]) {
+  int err;
+  iMesh_Instance m;
+  iMesh_newMesh(\"\",&m,&err,0);CHK(err);
+  iMesh_dtor(m,&err);CHK(err);
+  return 0;
+}
+")
+find_path (imesh_include_tmp iMeshP.h HINTS ${ITAPS_MESH_INCLUDES} NO_DEFAULT_PATH)
+if (imesh_include_tmp)
+  set (ITAPS_MESH_HAS_PARALLEL "YES")
+else ()
+  set (ITAPS_MESH_HAS_PARALLEL "NO")
+endif ()
+set (imesh_include_tmp "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+
+set (itaps_rel_libs)       # Extra libraries which should only be set when linking with iRel
+
+itaps_handle_component (GEOM iGeom "
+/* iGeom test program */
+#include <iGeom.h>
+#define CHK(err) if (err) return 1
+int main() {
+  int ierr;
+  iGeom_Instance g;
+  iGeom_newGeom(\"\",&g,&ierr,0);CHK(ierr);
+  iGeom_dtor(g,&ierr);CHK(ierr);
+  return 0;
+}
+")
+
+if (ITAPS_MESH_FOUND AND ITAPS_GEOM_FOUND) # iRel only makes sense if iMesh and iGeom are found
+  set (itaps_rel_libs "${ITAPS_MESH_LIBRARIES}" "${ITAPS_GEOM_LIBRARIES}")
+  itaps_handle_component (REL iRel "
+/* iRel test program */
+#include <iRel.h>
+#define CHK(err) if (err) return 1
+int main() {
+  int ierr;
+  iRel_Instance rel;
+  iRel_create(\"\",&rel,&ierr,0);CHK(ierr);
+  iRel_destroy(rel,&ierr);CHK(ierr);
+  return 0;
+}
+")
+endif ()
+
+set (ITAPS_INCLUDES)
+set (ITAPS_LIBRARIES)
+foreach (component REL GEOM MESH)
+  if (ITAPS_${component}_INCLUDES)
+    list (APPEND ITAPS_INCLUDES "${ITAPS_${component}_INCLUDES}")
+  endif ()
+  if (ITAPS_${component}_LIBRARIES)
+    list (APPEND ITAPS_LIBRARIES "${ITAPS_${component}_LIBRARIES}")
+  endif ()
+  message (STATUS "ITAPS_${component}: ${ITAPS_${component}_INCLUDES} ${ITAPS_${component}_LIBRARIES}")
+endforeach()
+list (REMOVE_DUPLICATES ITAPS_INCLUDES)
+list (REMOVE_DUPLICATES ITAPS_LIBRARIES)
+
+set (ITAPS_FOUND_REQUIRED_COMPONENTS YES)
+if (ITAPS_FIND_REQUIRED)
+  foreach (component ${ITAPS_FIND_COMPONENTS})
+    if (NOT ITAPS_${component}_FOUND)
+      set (ITAPS_FOUND_REQUIRED_COMPONENTS NOTFOUND)
+    endif()
+  endforeach()
+endif()
+
+message (STATUS "ITAPS: ${ITAPS_INCLUDES}  ${ITAPS_LIBRARIES}")
+
+find_package_handle_standard_args (ITAPS "ITAPS not found, check environment variables I{MESH,GEOM,REL}_DIR"
+  ITAPS_INCLUDES ITAPS_LIBRARIES ITAPS_FOUND_REQUIRED_COMPONENTS)
diff --git a/packages/cmake-modules/FindNetCDF.cmake b/packages/cmake-modules/FindNetCDF.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..4db4b0d384c7c2d2bfb57b4847bbb9d5c0650364
--- /dev/null
+++ b/packages/cmake-modules/FindNetCDF.cmake
@@ -0,0 +1,71 @@
+# - Find NetCDF
+# Find the native NetCDF includes and library
+#
+#  NETCDF_INCLUDES    - where to find netcdf.h, etc
+#  NETCDF_LIBRARIES   - Link these libraries when using NetCDF
+#  NETCDF_FOUND       - True if NetCDF found including required interfaces (see below)
+#
+# Your package can require certain interfaces to be FOUND by setting these
+#
+#  NETCDF_CXX         - require the C++ interface and link the C++ library
+#  NETCDF_F77         - require the F77 interface and link the fortran library
+#  NETCDF_F90         - require the F90 interface and link the fortran library
+#
+# The following are not for general use and are included in
+# NETCDF_LIBRARIES if the corresponding option above is set.
+#
+#  NETCDF_LIBRARIES_C    - Just the C interface
+#  NETCDF_LIBRARIES_CXX  - C++ interface, if available
+#  NETCDF_LIBRARIES_F77  - Fortran 77 interface, if available
+#  NETCDF_LIBRARIES_F90  - Fortran 90 interface, if available
+#
+# Normal usage would be:
+#  set (NETCDF_F90 "YES")
+#  find_package (NetCDF REQUIRED)
+#  target_link_libraries (uses_f90_interface ${NETCDF_LIBRARIES})
+#  target_link_libraries (only_uses_c_interface ${NETCDF_LIBRARIES_C})
+
+if (NETCDF_INCLUDES AND NETCDF_LIBRARIES)
+  # Already in cache, be silent
+  set (NETCDF_FIND_QUIETLY TRUE)
+endif (NETCDF_INCLUDES AND NETCDF_LIBRARIES)
+
+find_path (NETCDF_INCLUDES netcdf.h
+  HINTS NETCDF_DIR ENV NETCDF_DIR)
+
+find_library (NETCDF_LIBRARIES_C       NAMES netcdf)
+mark_as_advanced(NETCDF_LIBRARIES_C)
+
+set (NetCDF_has_interfaces "YES") # will be set to NO if we're missing any interfaces
+set (NetCDF_libs "${NETCDF_LIBRARIES_C}")
+
+get_filename_component (NetCDF_lib_dirs "${NETCDF_LIBRARIES_C}" PATH)
+
+macro (NetCDF_check_interface lang header libs)
+  if (NETCDF_${lang})
+    find_path (NETCDF_INCLUDES_${lang} NAMES ${header}
+      HINTS "${NETCDF_INCLUDES}" NO_DEFAULT_PATH)
+    find_library (NETCDF_LIBRARIES_${lang} NAMES ${libs}
+      HINTS "${NetCDF_lib_dirs}" NO_DEFAULT_PATH)
+    mark_as_advanced (NETCDF_INCLUDES_${lang} NETCDF_LIBRARIES_${lang})
+    if (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang})
+      list (INSERT NetCDF_libs 0 ${NETCDF_LIBRARIES_${lang}}) # prepend so that -lnetcdf is last
+    else (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang})
+      set (NetCDF_has_interfaces "NO")
+      message (STATUS "Failed to find NetCDF interface for ${lang}")
+    endif (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang})
+  endif (NETCDF_${lang})
+endmacro (NetCDF_check_interface)
+
+NetCDF_check_interface (CXX netcdfcpp.h netcdf_c++)
+NetCDF_check_interface (F77 netcdf.inc  netcdff)
+NetCDF_check_interface (F90 netcdf.mod  netcdff)
+
+set (NETCDF_LIBRARIES "${NetCDF_libs}" CACHE STRING "All NetCDF libraries required for interface level")
+
+# handle the QUIETLY and REQUIRED arguments and set NETCDF_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (NetCDF DEFAULT_MSG NETCDF_LIBRARIES NETCDF_INCLUDES NetCDF_has_interfaces)
+
+mark_as_advanced (NETCDF_LIBRARIES NETCDF_INCLUDES)
diff --git a/packages/cmake-modules/FindPETSc.cmake b/packages/cmake-modules/FindPETSc.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..41e02d9798c47f32ae57e7ecff2f544eb5256c1d
--- /dev/null
+++ b/packages/cmake-modules/FindPETSc.cmake
@@ -0,0 +1,346 @@
+# - Try to find PETSc
+# Once done this will define
+#
+#  PETSC_FOUND        - system has PETSc
+#  PETSC_INCLUDES     - the PETSc include directories
+#  PETSC_LIBRARIES    - Link these to use PETSc
+#  PETSC_COMPILER     - Compiler used by PETSc, helpful to find a compatible MPI
+#  PETSC_DEFINITIONS  - Compiler switches for using PETSc
+#  PETSC_MPIEXEC      - Executable for running MPI programs
+#  PETSC_VERSION      - Version string (MAJOR.MINOR.SUBMINOR)
+#
+#  Usage:
+#  find_package(PETSc COMPONENTS CXX)  - required if build --with-clanguage=C++ --with-c-support=0
+#  find_package(PETSc COMPONENTS C)    - standard behavior of checking build using a C compiler
+#  find_package(PETSc)                 - same as above
+#
+# Setting these changes the behavior of the search
+#  PETSC_DIR - directory in which PETSc resides
+#  PETSC_ARCH - build architecture
+#
+# Redistribution and use is allowed according to the terms of the BSD license.
+# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
+#
+
+cmake_policy(VERSION 3.3)
+
+set(PETSC_VALID_COMPONENTS
+  C
+  CXX)
+
+if(NOT PETSc_FIND_COMPONENTS)
+  get_property (_enabled_langs GLOBAL PROPERTY ENABLED_LANGUAGES)
+  if ("C" IN_LIST _enabled_langs)
+    set(PETSC_LANGUAGE_BINDINGS "C")
+  else ()
+    set(PETSC_LANGUAGE_BINDINGS "CXX")
+  endif ()
+else()
+  # Right now, this is designed for compatability with the --with-clanguage option, so
+  # only allow one item in the components list.
+  list(LENGTH ${PETSc_FIND_COMPONENTS} components_length)
+  if(${components_length} GREATER 1)
+    message(FATAL_ERROR "Only one component for PETSc is allowed to be specified")
+  endif()
+  # This is a stub for allowing multiple components should that time ever come. Perhaps
+  # to also test Fortran bindings?
+  foreach(component ${PETSc_FIND_COMPONENTS})
+    list(FIND PETSC_VALID_COMPONENTS ${component} component_location)
+    if(${component_location} EQUAL -1)
+      message(FATAL_ERROR "\"${component}\" is not a valid PETSc component.")
+    else()
+      list(APPEND PETSC_LANGUAGE_BINDINGS ${component})
+    endif()
+  endforeach()
+endif()
+
+function (petsc_get_version)
+  if (EXISTS "${PETSC_DIR}/include/petscversion.h")
+    file (STRINGS "${PETSC_DIR}/include/petscversion.h" vstrings REGEX "#define PETSC_VERSION_(RELEASE|MAJOR|MINOR|SUBMINOR|PATCH) ")
+    foreach (line ${vstrings})
+      string (REGEX REPLACE " +" ";" fields ${line}) # break line into three fields (the first is always "#define")
+      list (GET fields 1 var)
+      list (GET fields 2 val)
+      set (${var} ${val} PARENT_SCOPE)
+      set (${var} ${val})         # Also in local scope so we have access below
+    endforeach ()
+    if (PETSC_VERSION_RELEASE)
+      if ($(PETSC_VERSION_PATCH) GREATER 0)
+        set (PETSC_VERSION "${PETSC_VERSION_MAJOR}.${PETSC_VERSION_MINOR}.${PETSC_VERSION_SUBMINOR}p${PETSC_VERSION_PATCH}" CACHE INTERNAL "PETSc version")
+      else ()
+        set (PETSC_VERSION "${PETSC_VERSION_MAJOR}.${PETSC_VERSION_MINOR}.${PETSC_VERSION_SUBMINOR}" CACHE INTERNAL "PETSc version")
+      endif ()
+    else ()
+      # make dev version compare higher than any patch level of a released version
+      set (PETSC_VERSION "${PETSC_VERSION_MAJOR}.${PETSC_VERSION_MINOR}.${PETSC_VERSION_SUBMINOR}.99" CACHE INTERNAL "PETSc version")
+    endif ()
+  else ()
+    message (SEND_ERROR "PETSC_DIR can not be used, ${PETSC_DIR}/include/petscversion.h does not exist")
+  endif ()
+endfunction ()
+
+# Debian uses versioned paths e.g /usr/lib/petscdir/3.5/
+file (GLOB DEB_PATHS "/usr/lib/petscdir/*")
+
+find_path (PETSC_DIR include/petsc.h
+  HINTS ENV PETSC_DIR
+  PATHS
+  /usr/lib/petsc
+  # Debian paths
+  ${DEB_PATHS}
+  # Arch Linux path
+  /opt/petsc/linux-c-opt
+  # MacPorts path
+  /opt/local/lib/petsc
+  $ENV{HOME}/petsc
+  DOC "PETSc Directory")
+
+find_program (MAKE_EXECUTABLE NAMES make gmake)
+
+if (PETSC_DIR AND NOT PETSC_ARCH)
+  set (_petsc_arches
+    $ENV{PETSC_ARCH}                   # If set, use environment variable first
+    linux-gnu-c-debug linux-gnu-c-opt  # Debian defaults
+    x86_64-unknown-linux-gnu i386-unknown-linux-gnu)
+  set (petscconf "NOTFOUND" CACHE FILEPATH "Cleared" FORCE)
+  foreach (arch ${_petsc_arches})
+    if (NOT PETSC_ARCH)
+      find_path (petscconf petscconf.h
+        HINTS ${PETSC_DIR}
+        PATH_SUFFIXES ${arch}/include bmake/${arch}
+        NO_DEFAULT_PATH)
+      if (petscconf)
+        set (PETSC_ARCH "${arch}" CACHE STRING "PETSc build architecture")
+      endif (petscconf)
+    endif (NOT PETSC_ARCH)
+  endforeach (arch)
+  set (petscconf "NOTFOUND" CACHE INTERNAL "Scratch variable" FORCE)
+endif (PETSC_DIR AND NOT PETSC_ARCH)
+
+set (petsc_slaves LIBRARIES_SYS LIBRARIES_VEC LIBRARIES_MAT LIBRARIES_DM LIBRARIES_KSP LIBRARIES_SNES LIBRARIES_TS
+  INCLUDE_DIR INCLUDE_CONF)
+include (FindPackageMultipass)
+find_package_multipass (PETSc petsc_config_current
+  STATES DIR ARCH
+  DEPENDENTS INCLUDES LIBRARIES COMPILER MPIEXEC ${petsc_slaves})
+
+# Determine whether the PETSc layout is old-style (through 2.3.3) or
+# new-style (>= 3.0.0)
+if (EXISTS "${PETSC_DIR}/${PETSC_ARCH}/lib/petsc/conf/petscvariables") # > 3.5
+  set (petsc_conf_rules "${PETSC_DIR}/lib/petsc/conf/rules")
+  set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables")
+elseif (EXISTS "${PETSC_DIR}/${PETSC_ARCH}/include/petscconf.h")   # > 2.3.3
+  set (petsc_conf_rules "${PETSC_DIR}/conf/rules")
+  set (petsc_conf_variables "${PETSC_DIR}/conf/variables")
+elseif (EXISTS "${PETSC_DIR}/bmake/${PETSC_ARCH}/petscconf.h") # <= 2.3.3
+  set (petsc_conf_rules "${PETSC_DIR}/bmake/common/rules")
+  set (petsc_conf_variables "${PETSC_DIR}/bmake/common/variables")
+elseif (PETSC_DIR)
+  message (SEND_ERROR "The pair PETSC_DIR=${PETSC_DIR} PETSC_ARCH=${PETSC_ARCH} do not specify a valid PETSc installation")
+endif ()
+
+if (petsc_conf_rules AND petsc_conf_variables AND NOT petsc_config_current)
+  petsc_get_version()
+
+  # Put variables into environment since they are needed to get
+  # configuration (petscvariables) in the PETSc makefile
+  set (ENV{PETSC_DIR} "${PETSC_DIR}")
+  set (ENV{PETSC_ARCH} "${PETSC_ARCH}")
+
+  # A temporary makefile to probe the PETSc configuration
+  set (petsc_config_makefile "${PROJECT_BINARY_DIR}/Makefile.petsc")
+  file (WRITE "${petsc_config_makefile}"
+"## This file was autogenerated by FindPETSc.cmake
+# PETSC_DIR  = ${PETSC_DIR}
+# PETSC_ARCH = ${PETSC_ARCH}
+include ${petsc_conf_rules}
+include ${petsc_conf_variables}
+show :
+\t-@echo -n \${\${VARIABLE}}
+")
+
+  macro (PETSC_GET_VARIABLE name var)
+    set (${var} "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+    execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} show VARIABLE=${name}
+      OUTPUT_VARIABLE ${var}
+      RESULT_VARIABLE petsc_return)
+  endmacro (PETSC_GET_VARIABLE)
+  petsc_get_variable (PETSC_LIB_DIR            petsc_lib_dir)
+  petsc_get_variable (PETSC_EXTERNAL_LIB_BASIC petsc_libs_external)
+  petsc_get_variable (PETSC_CCPPFLAGS          petsc_cpp_line)
+  petsc_get_variable (PETSC_INCLUDE            petsc_include)
+  petsc_get_variable (PCC                      petsc_cc)
+  petsc_get_variable (PCC_FLAGS                petsc_cc_flags)
+  petsc_get_variable (MPIEXEC                  petsc_mpiexec)
+  # We are done with the temporary Makefile, calling PETSC_GET_VARIABLE after this point is invalid!
+  file (REMOVE ${petsc_config_makefile})
+
+  include (ResolveCompilerPaths)
+  # Extract include paths and libraries from compile command line
+  resolve_includes (petsc_includes_all "${petsc_cpp_line}")
+
+  #on windows we need to make sure we're linking against the right
+  #runtime library
+  if (WIN32)
+    if (petsc_cc_flags MATCHES "-MT")
+      set(using_md False)
+      foreach(flag_var
+          CMAKE_C_FLAGS CMAKE_C_FLAGS_DEBUG CMAKE_C_FLAGS_RELEASE
+          CMAKE_C_FLAGS_MINSIZEREL CMAKE_C_FLAGS_RELWITHDEBINFO
+          CMAKE_CXX_FLAGS CMAKE_CXX_FLAGS_DEBUG CMAKE_CXX_FLAGS_RELEASE
+          CMAKE_CXX_FLAGS_MINSIZEREL CMAKE_CXX_FLAGS_RELWITHDEBINFO)
+        if(${flag_var} MATCHES "/MD")
+          set(using_md True)
+        endif(${flag_var} MATCHES "/MD")
+      endforeach(flag_var)
+      if(${using_md} MATCHES "True")
+        message(WARNING "PETSc was built with /MT, but /MD is currently set.
+ See http://www.cmake.org/Wiki/CMake_FAQ#How_can_I_build_my_MSVC_application_with_a_static_runtime.3F")
+      endif(${using_md} MATCHES "True")
+    endif (petsc_cc_flags MATCHES "-MT")
+  endif (WIN32)
+
+  include (CorrectWindowsPaths)
+  convert_cygwin_path(petsc_lib_dir)
+  message (STATUS "petsc_lib_dir ${petsc_lib_dir}")
+
+  macro (PETSC_FIND_LIBRARY suffix name)
+    set (PETSC_LIBRARY_${suffix} "NOTFOUND" CACHE INTERNAL "Cleared" FORCE) # Clear any stale value, if we got here, we need to find it again
+    if (WIN32)
+      set (libname lib${name}) #windows expects "libfoo", linux expects "foo"
+    else (WIN32)
+      set (libname ${name})
+    endif (WIN32)
+    find_library (PETSC_LIBRARY_${suffix} NAMES ${libname} HINTS ${petsc_lib_dir} NO_DEFAULT_PATH)
+    set (PETSC_LIBRARIES_${suffix} "${PETSC_LIBRARY_${suffix}}")
+    mark_as_advanced (PETSC_LIBRARY_${suffix})
+  endmacro (PETSC_FIND_LIBRARY suffix name)
+
+  # Look for petscvec first, if it doesn't exist, we must be using single-library
+  petsc_find_library (VEC petscvec)
+  if (PETSC_LIBRARY_VEC)
+    petsc_find_library (SYS  "petscsys;petsc") # libpetscsys is called libpetsc prior to 3.1 (when single-library was introduced)
+    petsc_find_library (MAT  petscmat)
+    petsc_find_library (DM   petscdm)
+    petsc_find_library (KSP  petscksp)
+    petsc_find_library (SNES petscsnes)
+    petsc_find_library (TS   petscts)
+    macro (PETSC_JOIN libs deps)
+      list (APPEND PETSC_LIBRARIES_${libs} ${PETSC_LIBRARIES_${deps}})
+    endmacro (PETSC_JOIN libs deps)
+    petsc_join (VEC  SYS)
+    petsc_join (MAT  VEC)
+    petsc_join (DM   MAT)
+    petsc_join (KSP  DM)
+    petsc_join (SNES KSP)
+    petsc_join (TS   SNES)
+    petsc_join (ALL  TS)
+  else ()
+    set (PETSC_LIBRARY_VEC "NOTFOUND" CACHE INTERNAL "Cleared" FORCE) # There is no libpetscvec
+    petsc_find_library (SINGLE petsc)
+    # Debian 9/Ubuntu 16.04 uses _real and _complex extensions when using libraries in /usr/lib/petsc.
+    if (NOT PETSC_LIBRARY_SINGLE)
+      petsc_find_library (SINGLE petsc_real)
+    endif()
+    if (NOT PETSC_LIBRARY_SINGLE)
+      petsc_find_library (SINGLE petsc_complex)
+    endif()
+    foreach (pkg SYS VEC MAT DM KSP SNES TS ALL)
+      set (PETSC_LIBRARIES_${pkg} "${PETSC_LIBRARY_SINGLE}")
+    endforeach ()
+  endif ()
+  if (PETSC_LIBRARY_TS)
+    message (STATUS "Recognized PETSc install with separate libraries for each package")
+  else ()
+    message (STATUS "Recognized PETSc install with single library for all packages")
+  endif ()
+
+  include(Check${PETSC_LANGUAGE_BINDINGS}SourceRuns)
+  macro (PETSC_TEST_RUNS includes libraries runs)
+    if (PETSC_VERSION VERSION_GREATER 3.1)
+      set (_PETSC_TSDestroy "TSDestroy(&ts)")
+    else ()
+      set (_PETSC_TSDestroy "TSDestroy(ts)")
+    endif ()
+
+    set(_PETSC_TEST_SOURCE "
+static const char help[] = \"PETSc test program.\";
+#include <petscts.h>
+int main(int argc,char *argv[]) {
+  PetscErrorCode ierr;
+  TS ts;
+
+  ierr = PetscInitialize(&argc,&argv,0,help);CHKERRQ(ierr);
+  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+  ierr = ${_PETSC_TSDestroy};CHKERRQ(ierr);
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}
+")
+    multipass_source_runs ("${includes}" "${libraries}" "${_PETSC_TEST_SOURCE}" ${runs} "${PETSC_LANGUAGE_BINDINGS}")
+    if (${${runs}})
+      set (PETSC_EXECUTABLE_RUNS "YES" CACHE BOOL
+        "Can the system successfully run a PETSc executable?  This variable can be manually set to \"YES\" to force CMake to accept a given PETSc configuration, but this will almost always result in a broken build.  If you change PETSC_DIR, PETSC_ARCH, or PETSC_CURRENT you would have to reset this variable." FORCE)
+    endif (${${runs}})
+  endmacro (PETSC_TEST_RUNS)
+
+
+  find_path (PETSC_INCLUDE_DIR petscts.h HINTS "${PETSC_DIR}" PATH_SUFFIXES include NO_DEFAULT_PATH)
+  find_path (PETSC_INCLUDE_CONF petscconf.h HINTS "${PETSC_DIR}" PATH_SUFFIXES "${PETSC_ARCH}/include" "bmake/${PETSC_ARCH}" NO_DEFAULT_PATH)
+  mark_as_advanced (PETSC_INCLUDE_DIR PETSC_INCLUDE_CONF)
+  set (petsc_includes_minimal ${PETSC_INCLUDE_CONF} ${PETSC_INCLUDE_DIR})
+
+  petsc_test_runs ("${petsc_includes_minimal}" "${PETSC_LIBRARIES_TS}" petsc_works_minimal)
+  if (petsc_works_minimal)
+    message (STATUS "Minimal PETSc includes and libraries work.  This probably means we are building with shared libs.")
+    set (petsc_includes_needed "${petsc_includes_minimal}")
+  else (petsc_works_minimal)     # Minimal includes fail, see if just adding full includes fixes it
+    petsc_test_runs ("${petsc_includes_all}" "${PETSC_LIBRARIES_TS}" petsc_works_allincludes)
+    if (petsc_works_allincludes) # It does, we just need all the includes (
+      message (STATUS "PETSc requires extra include paths, but links correctly with only interface libraries.  This is an unexpected configuration (but it seems to work fine).")
+      set (petsc_includes_needed ${petsc_includes_all})
+    else (petsc_works_allincludes) # We are going to need to link the external libs explicitly
+      resolve_libraries (petsc_libraries_external "${petsc_libs_external}")
+      foreach (pkg SYS VEC MAT DM KSP SNES TS ALL)
+        list (APPEND PETSC_LIBRARIES_${pkg}  ${petsc_libraries_external})
+      endforeach (pkg)
+      petsc_test_runs ("${petsc_includes_minimal}" "${PETSC_LIBRARIES_TS}" petsc_works_alllibraries)
+      if (petsc_works_alllibraries)
+         message (STATUS "PETSc only need minimal includes, but requires explicit linking to all dependencies.  This is expected when PETSc is built with static libraries.")
+        set (petsc_includes_needed ${petsc_includes_minimal})
+      else (petsc_works_alllibraries)
+        # It looks like we really need everything, should have listened to Matt
+        set (petsc_includes_needed ${petsc_includes_all})
+        petsc_test_runs ("${petsc_includes_all}" "${PETSC_LIBRARIES_TS}" petsc_works_all)
+        if (petsc_works_all) # We fail anyways
+          message (STATUS "PETSc requires extra include paths and explicit linking to all dependencies.  This probably means you have static libraries and something unexpected in PETSc headers.")
+        else (petsc_works_all) # We fail anyways
+          message (STATUS "PETSc could not be used, maybe the install is broken.")
+        endif (petsc_works_all)
+      endif (petsc_works_alllibraries)
+    endif (petsc_works_allincludes)
+  endif (petsc_works_minimal)
+
+  # We do an out-of-source build so __FILE__ will be an absolute path, hence __INSDIR__ is superfluous
+  if (${PETSC_VERSION} VERSION_LESS 3.1)
+    set (PETSC_DEFINITIONS "-D__SDIR__=\"\"" CACHE STRING "PETSc definitions" FORCE)
+  else ()
+    set (PETSC_DEFINITIONS "-D__INSDIR__=" CACHE STRING "PETSc definitions" FORCE)
+  endif ()
+  # Sometimes this can be used to assist FindMPI.cmake
+  set (PETSC_MPIEXEC ${petsc_mpiexec} CACHE FILEPATH "Executable for running PETSc MPI programs" FORCE)
+  set (PETSC_INCLUDES ${petsc_includes_needed} CACHE STRING "PETSc include path" FORCE)
+  set (PETSC_LIBRARIES ${PETSC_LIBRARIES_ALL} CACHE STRING "PETSc libraries" FORCE)
+  set (PETSC_COMPILER ${petsc_cc} CACHE FILEPATH "PETSc compiler" FORCE)
+  # Note that we have forced values for all these choices.  If you
+  # change these, you are telling the system to trust you that they
+  # work.  It is likely that you will end up with a broken build.
+  mark_as_advanced (PETSC_INCLUDES PETSC_LIBRARIES PETSC_COMPILER PETSC_DEFINITIONS PETSC_MPIEXEC PETSC_EXECUTABLE_RUNS)
+endif ()
+
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (PETSc
+  REQUIRED_VARS PETSC_INCLUDES PETSC_LIBRARIES PETSC_EXECUTABLE_RUNS
+  VERSION_VAR PETSC_VERSION
+  FAIL_MESSAGE "PETSc could not be found.  Be sure to set PETSC_DIR and PETSC_ARCH.")
diff --git a/packages/cmake-modules/FindPackageMultipass.cmake b/packages/cmake-modules/FindPackageMultipass.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..fbf06a7f0fc3aa20a0387f091eac4f74e7ffdab2
--- /dev/null
+++ b/packages/cmake-modules/FindPackageMultipass.cmake
@@ -0,0 +1,106 @@
+# PackageMultipass - this module defines two macros
+#
+# FIND_PACKAGE_MULTIPASS (Name CURRENT
+#  STATES VAR0 VAR1 ...
+#  DEPENDENTS DEP0 DEP1 ...)
+#
+#  This function creates a cache entry <UPPERCASED-Name>_CURRENT which
+#  the user can set to "NO" to trigger a reconfiguration of the package.
+#  The first time this function is called, the values of
+#  <UPPERCASED-Name>_VAR0, ... are saved.  If <UPPERCASED-Name>_CURRENT
+#  is false or if any STATE has changed since the last time
+#  FIND_PACKAGE_MULTIPASS() was called, then CURRENT will be set to "NO",
+#  otherwise CURRENT will be "YES".  IF not CURRENT, then
+#  <UPPERCASED-Name>_DEP0, ... will be FORCED to NOTFOUND.
+#  Example:
+#    find_path (FOO_DIR include/foo.h)
+#    FIND_PACKAGE_MULTIPASS (Foo foo_current
+#      STATES DIR
+#      DEPENDENTS INCLUDES LIBRARIES)
+#    if (NOT foo_current)
+#      # Make temporary files, run programs, etc, to determine FOO_INCLUDES and FOO_LIBRARIES
+#    endif (NOT foo_current)
+#
+# MULTIPASS_SOURCE_RUNS (Name INCLUDES LIBRARIES SOURCE RUNS LANGUAGE)
+#  Always runs the given test, use this when you need to re-run tests
+#  because parent variables have made old cache entries stale. The LANGUAGE
+#  variable is either C or CXX indicating which compiler the test should
+#  use.
+# MULTIPASS_C_SOURCE_RUNS (Name INCLUDES LIBRARIES SOURCE RUNS)
+#  DEPRECATED! This is only included for backwards compatability. Use
+#  the more general MULTIPASS_SOURCE_RUNS instead.
+#  Always runs the given test, use this when you need to re-run tests
+#  because parent variables have made old cache entries stale.
+
+macro (FIND_PACKAGE_MULTIPASS _name _current)
+  string (TOUPPER ${_name} _NAME)
+  set (_args ${ARGV})
+  list (REMOVE_AT _args 0 1)
+
+  set (_states_current "YES")
+  list (GET _args 0 _cmd)
+  if (_cmd STREQUAL "STATES")
+    list (REMOVE_AT _args 0)
+    list (GET _args 0 _state)
+    while (_state AND NOT _state STREQUAL "DEPENDENTS")
+      # The name of the stored value for the given state
+      set (_stored_var PACKAGE_MULTIPASS_${_NAME}_${_state})
+      if (NOT "${${_stored_var}}" STREQUAL "${${_NAME}_${_state}}")
+        set (_states_current "NO")
+      endif (NOT "${${_stored_var}}" STREQUAL "${${_NAME}_${_state}}")
+      set (${_stored_var} "${${_NAME}_${_state}}" CACHE INTERNAL "Stored state for ${_name}." FORCE)
+      list (REMOVE_AT _args 0)
+      list (GET _args 0 _state)
+    endwhile (_state AND NOT _state STREQUAL "DEPENDENTS")
+  endif (_cmd STREQUAL "STATES")
+
+  set (_stored ${_NAME}_CURRENT)
+  if (NOT ${_stored})
+    set (${_stored} "YES" CACHE BOOL "Is the configuration for ${_name} current?  Set to \"NO\" to reconfigure." FORCE)
+    set (_states_current "NO")
+  endif (NOT ${_stored})
+
+  set (${_current} ${_states_current})
+  if (NOT ${_current} AND PACKAGE_MULTIPASS_${_name}_CALLED)
+    message (STATUS "Clearing ${_name} dependent variables")
+    # Clear all the dependent variables so that the module can reset them
+    list (GET _args 0 _cmd)
+    if (_cmd STREQUAL "DEPENDENTS")
+      list (REMOVE_AT _args 0)
+      foreach (dep ${_args})
+        set (${_NAME}_${dep} "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+      endforeach (dep)
+    endif (_cmd STREQUAL "DEPENDENTS")
+    set (${_NAME}_FOUND "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
+  endif ()
+  set (PACKAGE_MULTIPASS_${name}_CALLED YES CACHE INTERNAL "Private" FORCE)
+endmacro (FIND_PACKAGE_MULTIPASS)
+
+
+macro (MULTIPASS_SOURCE_RUNS includes libraries source runs language)
+  include (Check${language}SourceRuns)
+  # This is a ridiculous hack.  CHECK_${language}_SOURCE_* thinks that if the
+  # *name* of the return variable doesn't change, then the test does
+  # not need to be re-run.  We keep an internal count which we
+  # increment to guarantee that every test name is unique.  If we've
+  # gotten here, then the configuration has changed enough that the
+  # test *needs* to be rerun.
+  if (NOT MULTIPASS_TEST_COUNT)
+    set (MULTIPASS_TEST_COUNT 00)
+  endif (NOT MULTIPASS_TEST_COUNT)
+  math (EXPR _tmp "${MULTIPASS_TEST_COUNT} + 1") # Why can't I add to a cache variable?
+  set (MULTIPASS_TEST_COUNT ${_tmp} CACHE INTERNAL "Unique test ID")
+  set (testname MULTIPASS_TEST_${MULTIPASS_TEST_COUNT}_${runs})
+  set (CMAKE_REQUIRED_INCLUDES ${includes})
+  set (CMAKE_REQUIRED_LIBRARIES ${libraries})
+  if(${language} STREQUAL "C")
+    check_c_source_runs ("${source}" ${testname})
+  elseif(${language} STREQUAL "CXX")
+    check_cxx_source_runs ("${source}" ${testname})
+  endif()
+  set (${runs} "${${testname}}")
+endmacro (MULTIPASS_SOURCE_RUNS)
+
+macro (MULTIPASS_C_SOURCE_RUNS includes libraries source runs)
+  multipass_source_runs("${includes}" "${libraries}" "${source}" ${runs} "C")
+endmacro (MULTIPASS_C_SOURCE_RUNS)
diff --git a/packages/cmake-modules/LICENSE b/packages/cmake-modules/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..49bac9892af49f867817ca73496499cd6a6fbf89
--- /dev/null
+++ b/packages/cmake-modules/LICENSE
@@ -0,0 +1,23 @@
+Copyright $(git shortlog -s)
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice, this
+  list of conditions and the following disclaimer in the documentation and/or
+  other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/packages/cmake-modules/README b/packages/cmake-modules/README
new file mode 100644
index 0000000000000000000000000000000000000000..6f2183ca12ec2ea808f9f69b53e1f1b0a9bd1c97
--- /dev/null
+++ b/packages/cmake-modules/README
@@ -0,0 +1,7 @@
+A collection of CMake modules, these can mostly be used independently.
+
+The utilities for writing robust Find* modules might be useful until
+CMake takes static libraries and multiple active configurations
+seriously.
+
+  http://www.cmake.org/Wiki/CMake:Improving_Find*_Modules
diff --git a/packages/cmake-modules/ResolveCompilerPaths.cmake b/packages/cmake-modules/ResolveCompilerPaths.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..54787fa38ffa50136414e6c788c50fb3c63746b8
--- /dev/null
+++ b/packages/cmake-modules/ResolveCompilerPaths.cmake
@@ -0,0 +1,105 @@
+# ResolveCompilerPaths - this module defines two macros
+#
+# RESOLVE_LIBRARIES (XXX_LIBRARIES LINK_LINE)
+#  This macro is intended to be used by FindXXX.cmake modules.
+#  It parses a compiler link line and resolves all libraries
+#  (-lfoo) using the library path contexts (-L/path) in scope.
+#  The result in XXX_LIBRARIES is the list of fully resolved libs.
+#  Example:
+#
+#    RESOLVE_LIBRARIES (FOO_LIBRARIES "-L/A -la -L/B -lb -lc -ld")
+#
+#  will be resolved to
+#
+#    FOO_LIBRARIES:STRING="/A/liba.so;/B/libb.so;/A/libc.so;/usr/lib/libd.so"
+#
+#  if the filesystem looks like
+#
+#    /A:       liba.so         libc.so
+#    /B:       liba.so libb.so
+#    /usr/lib: liba.so libb.so libc.so libd.so
+#
+#  and /usr/lib is a system directory.
+#
+#  Note: If RESOLVE_LIBRARIES() resolves a link line differently from
+#  the native linker, there is a bug in this macro (please report it).
+#
+# RESOLVE_INCLUDES (XXX_INCLUDES INCLUDE_LINE)
+#  This macro is intended to be used by FindXXX.cmake modules.
+#  It parses a compile line and resolves all includes
+#  (-I/path/to/include) to a list of directories.  Other flags are ignored.
+#  Example:
+#
+#    RESOLVE_INCLUDES (FOO_INCLUDES "-I/A -DBAR='\"irrelevant -I/string here\"' -I/B")
+#
+#  will be resolved to
+#
+#    FOO_INCLUDES:STRING="/A;/B"
+#
+#  assuming both directories exist.
+#  Note: as currently implemented, the -I/string will be picked up mistakenly (cry, cry)
+include (CorrectWindowsPaths)
+
+macro (RESOLVE_LIBRARIES LIBS LINK_LINE)
+  string (REGEX MATCHALL "((-L|-l|-Wl)([^\" ]+|\"[^\"]+\")|[^\" ]+\\.(a|so|dll|lib))" _all_tokens "${LINK_LINE}")
+  set (_libs_found "")
+  set (_directory_list "")
+  foreach (token ${_all_tokens})
+    if (token MATCHES "-L([^\" ]+|\"[^\"]+\")")
+      # If it's a library path, add it to the list
+      string (REGEX REPLACE "^-L" "" token ${token})
+      string (REGEX REPLACE "//" "/" token ${token})
+      convert_cygwin_path(token)
+      list (APPEND _directory_list ${token})
+    elseif (token MATCHES "^(-l([^\" ]+|\"[^\"]+\")|[^\" ]+\\.(a|so|dll|lib))")
+      # It's a library, resolve the path by looking in the list and then (by default) in system directories
+      if (WIN32) #windows expects "libfoo", linux expects "foo"
+        string (REGEX REPLACE "^-l" "lib" token ${token})
+      else (WIN32)
+        string (REGEX REPLACE "^-l" "" token ${token})
+      endif (WIN32)
+      set (_root "")
+      if (token MATCHES "^/")	# We have an absolute path
+        #separate into a path and a library name:
+        string (REGEX MATCH "[^/]*\\.(a|so|dll|lib)$" libname ${token})
+        string (REGEX MATCH ".*[^${libname}$]" libpath ${token})
+        convert_cygwin_path(libpath)
+        set (_directory_list ${_directory_list} ${libpath})
+        set (token ${libname})
+      endif (token MATCHES "^/")
+      set (_lib "NOTFOUND" CACHE FILEPATH "Cleared" FORCE)
+      find_library (_lib ${token} HINTS ${_directory_list} ${_root})
+      if (_lib)
+	string (REPLACE "//" "/" _lib ${_lib})
+        list (APPEND _libs_found ${_lib})
+      else (_lib)
+        message (STATUS "Unable to find library ${token}")
+      endif (_lib)
+    endif (token MATCHES "-L([^\" ]+|\"[^\"]+\")")
+  endforeach (token)
+  set (_lib "NOTFOUND" CACHE INTERNAL "Scratch variable" FORCE)
+  # only the LAST occurence of each library is required since there should be no circular dependencies
+  if (_libs_found)
+    list (REVERSE _libs_found)
+    list (REMOVE_DUPLICATES _libs_found)
+    list (REVERSE _libs_found)
+  endif (_libs_found)
+  set (${LIBS} "${_libs_found}")
+endmacro (RESOLVE_LIBRARIES)
+
+macro (RESOLVE_INCLUDES INCS COMPILE_LINE)
+  string (REGEX MATCHALL "-I([^\" ]+|\"[^\"]+\")" _all_tokens "${COMPILE_LINE}")
+  set (_incs_found "")
+  foreach (token ${_all_tokens})
+    string (REGEX REPLACE "^-I" "" token ${token})
+    string (REGEX REPLACE "//" "/" token ${token})
+    convert_cygwin_path(token)
+    if (EXISTS ${token})
+      list (APPEND _incs_found ${token})
+    else (EXISTS ${token})
+      message (STATUS "Include directory ${token} does not exist")
+    endif (EXISTS ${token})
+  endforeach (token)
+  list (REMOVE_DUPLICATES _incs_found)
+  set (${INCS} "${_incs_found}")
+endmacro (RESOLVE_INCLUDES)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0d6e8cbb996d4f6c18c9098bb9af7051c11b0dbf..1b21f79a7545c442be5448792132acbdf840947a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -10,7 +10,7 @@ add_subdirectory(utils)
 add_subdirectory(language)
 
 # Pugs algebra
-#add_subdirectory(algebra)
+add_subdirectory(algebra)
 
 # Pugs mesh
 add_subdirectory(mesh)
diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index bc822668a7f63c4d5d7cc4c45f4d1b1d70084fc0..ff34b89e148d844215cf7fffc21c2ad78bca5ec9 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -7,16 +7,20 @@
 
 #include <rang.hpp>
 
-template <bool verbose = true>
 struct BiCGStab
 {
   template <typename VectorType, typename MatrixType>
-  BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
+  BiCGStab(const MatrixType& A,
+           VectorType& x,
+           const VectorType& b,
+           const double epsilon,
+           const size_t maximum_iteration,
+           const bool verbose)
   {
-    if constexpr (verbose) {
+    if (verbose) {
       std::cout << "- bi-conjugate gradient stabilized\n";
       std::cout << "  epsilon = " << epsilon << '\n';
-      std::cout << "  maximum number of iterations: " << max_iter << '\n';
+      std::cout << "  maximum number of iterations: " << maximum_iteration << '\n';
     }
 
     VectorType r_k_1{b.size()};
@@ -38,13 +42,13 @@ struct BiCGStab
 
       VectorType r_k{x.size()};
 
-      if constexpr (verbose) {
+      if (verbose) {
         std::cout << "   initial residu: " << resid0 << '\n';
       }
-      for (size_t i = 1; i <= max_iter; ++i) {
-        if constexpr (verbose) {
-          std::cout << "  - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
-                    << "\tabsolute: " << residu << '\n';
+      for (size_t i = 1; i <= maximum_iteration; ++i) {
+        if (verbose) {
+          std::cout << "  - iteration: " << std::setw(6) << i << " residu: " << std::scientific << residu / resid0
+                    << " absolute: " << std::scientific << residu << '\n';
         }
 
         Ap_k = A * p_k;
@@ -77,8 +81,8 @@ struct BiCGStab
                   << '\n';
         ;
         std::cout << "  - epsilon:          " << epsilon << '\n';
-        std::cout << "  - relative residu : " << residu / resid0 << '\n';
-        std::cout << "  - absolute residu : " << residu << '\n';
+        std::cout << "  - relative residu : " << std::scientific << residu / resid0 << '\n';
+        std::cout << "  - absolute residu : " << std::scientific << residu << '\n';
       }
     }
   }
diff --git a/src/algebra/PCG.hpp b/src/algebra/CG.hpp
similarity index 60%
rename from src/algebra/PCG.hpp
rename to src/algebra/CG.hpp
index 36a4635b313d2bba13d97e1b7a4ea1882e2a66a2..538c77571dac4d551d7476ea79bf855e8026fb53 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/CG.hpp
@@ -6,30 +6,30 @@
 
 #include <rang.hpp>
 
-template <bool verbose = true>
-struct PCG
+struct CG
 {
-  template <typename VectorType, typename MatrixType, typename PreconditionerType>
-  PCG(const VectorType& f,
-      const MatrixType& A,
-      [[maybe_unused]] const PreconditionerType& C,
-      VectorType& x,
-      const size_t maxiter,
-      const double epsilon = 1e-6)
+  template <typename VectorType, typename MatrixType>
+  CG(const MatrixType& A,
+     VectorType& x,
+     const VectorType& f,
+     const double epsilon,
+     const size_t maximum_iteration,
+     const bool verbose = false)
   {
-    if constexpr (verbose) {
+    if (verbose) {
       std::cout << "- conjugate gradient\n";
       std::cout << "  epsilon = " << epsilon << '\n';
-      std::cout << "  maximum number of iterations: " << maxiter << '\n';
+      std::cout << "  maximum number of iterations: " << maximum_iteration << '\n';
     }
 
     VectorType h{f.size()};
     VectorType b = copy(f);
 
-    if constexpr (verbose) {
+    if (verbose) {
       h = A * x;
       h -= f;
-      std::cout << "- initial *real* residu :   " << (h, h) << '\n';
+      std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu :   " << (h, h)
+                << '\n';
     }
 
     VectorType g{b.size()};
@@ -40,7 +40,7 @@ struct PCG
 
     double relativeEpsilon = epsilon;
 
-    for (size_t i = 1; i <= maxiter; ++i) {
+    for (size_t i = 1; i <= maximum_iteration; ++i) {
       if (i == 1) {
         h = A * x;
 
@@ -74,13 +74,13 @@ struct PCG
       if ((i == 1) && (gcg != 0)) {
         relativeEpsilon = epsilon * gcg;
         gcg0            = gcg;
-        if constexpr (verbose) {
+        if (verbose) {
           std::cout << "  initial residu: " << gcg << '\n';
         }
       }
-      if constexpr (verbose) {
-        std::cout << "  - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
-        std::cout << "\tabsolute: " << gcg << '\n';
+      if (verbose) {
+        std::cout << "  - iteration " << std::setw(6) << i << std::scientific << " residu: " << gcg / gcg0;
+        std::cout << " absolute: " << std::scientific << gcg << '\n';
       }
 
       if (gcg < relativeEpsilon) {
@@ -96,8 +96,8 @@ struct PCG
     if (gcg > relativeEpsilon) {
       std::cout << "  conjugate gradient: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset << '\n';
       std::cout << "  - epsilon:          " << epsilon << '\n';
-      std::cout << "  - relative residu : " << gcg / gcg0 << '\n';
-      std::cout << "  - absolute residu : " << gcg << '\n';
+      std::cout << "  - relative residu : " << std::scientific << gcg / gcg0 << '\n';
+      std::cout << "  - absolute residu : " << std::scientific << gcg << '\n';
     }
   }
 };
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d1f01ec310858be0e775514f934ce0a7c0e9da28
--- /dev/null
+++ b/src/algebra/CMakeLists.txt
@@ -0,0 +1,7 @@
+# ------------------- Source files --------------------
+
+add_library(
+  PugsAlgebra
+  LinearSolver.cpp
+  LinearSolverOptions.cpp
+  PETScWrapper.cpp)
diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 97a0a8a8ed44d660f9524eac269abee0324cd136..d9b03ac93fb2f1c4baeff98c90daac17547ff7cf 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -30,8 +30,21 @@ class CRSMatrix
     return m_host_matrix.numRows();
   }
 
+  auto
+  values() const
+  {
+    return m_values;
+  }
+
+  auto
+  hostMatrix() const
+  {
+    return m_host_matrix;
+  }
+
   template <typename DataType2>
-  Vector<MutableDataType> operator*(const Vector<DataType2>& x) const
+  Vector<MutableDataType>
+  operator*(const Vector<DataType2>& x) const
   {
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
                   "Cannot multiply matrix and vector of different type");
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a275cd5816d6ce04715fd1d160008b4dd8508051
--- /dev/null
+++ b/src/algebra/LinearSolver.cpp
@@ -0,0 +1,320 @@
+#include <algebra/LinearSolver.hpp>
+#include <utils/pugs_config.hpp>
+
+#include <algebra/BiCGStab.hpp>
+#include <algebra/CG.hpp>
+
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
+struct LinearSolver::Internals
+{
+  static bool
+  hasLibrary(const LSLibrary library)
+  {
+    switch (library) {
+    case LSLibrary::builtin: {
+      return true;
+    }
+    case LSLibrary::petsc: {
+#ifdef PUGS_HAS_PETSC
+      return true;
+#else
+      return false;
+#endif
+    }
+    default: {
+      throw UnexpectedError("Linear system library (" + ::name(library) + ") was not set!");
+    }
+    }
+  }
+
+  static void
+  checkHasLibrary(const LSLibrary library)
+  {
+    if (not hasLibrary(library)) {
+      throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!");
+    }
+  }
+
+  static void
+  checkBuiltinMethod(const LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg:
+    case LSMethod::bicgstab: {
+      break;
+    }
+    default: {
+      throw NormalError(name(method) + " is not a builtin linear solver!");
+    }
+    }
+  }
+
+  static void
+  checkPETScMethod(const LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg:
+    case LSMethod::bicgstab:
+    case LSMethod::bicgstab2:
+    case LSMethod::gmres:
+    case LSMethod::lu:
+    case LSMethod::choleski: {
+      break;
+    }
+    default: {
+      throw NormalError(name(method) + " is not a builtin linear solver!");
+    }
+    }
+  }
+
+  static void
+  checkBuiltinPrecond(const LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none: {
+      break;
+    }
+    default: {
+      throw NormalError(name(precond) + " is not a builtin preconditioner!");
+    }
+    }
+  }
+
+  static void
+  checkPETScPrecond(const LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none:
+    case LSPrecond::diagonal:
+    case LSPrecond::incomplete_choleski:
+    case LSPrecond::incomplete_LU: {
+      break;
+    }
+    default: {
+      throw NormalError(name(precond) + " is not a PETSc preconditioner!");
+    }
+    }
+  }
+
+  static void
+  checkOptions(const LinearSolverOptions& options)
+  {
+    switch (options.library()) {
+    case LSLibrary::builtin: {
+      checkBuiltinMethod(options.method());
+      checkBuiltinPrecond(options.precond());
+      break;
+    }
+    case LSLibrary::petsc: {
+      checkPETScMethod(options.method());
+      checkPETScPrecond(options.precond());
+      break;
+    }
+    default: {
+      throw UnexpectedError("undefined options compatibility for this library (" + ::name(options.library()) + ")!");
+    }
+    }
+  }
+
+  static void
+  builtinSolveLocalSystem(const CRSMatrix<double, size_t>& A,
+                          Vector<double>& x,
+                          const Vector<double>& b,
+                          const LinearSolverOptions& options)
+  {
+    if (options.precond() != LSPrecond::none) {
+      throw UnexpectedError("builtin linear solver do not allow any preconditioner!");
+    }
+    switch (options.method()) {
+    case LSMethod::cg: {
+      CG{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
+      break;
+    }
+    case LSMethod::bicgstab: {
+      BiCGStab{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
+      break;
+    }
+    default: {
+      throw NotImplementedError("undefined builtin method: " + name(options.method()));
+    }
+    }
+  }
+
+#ifdef PUGS_HAS_PETSC
+
+  static int
+  petscMonitor(KSP, int i, double residu, void*)
+  {
+    std::cout << "  - iteration: " << std::setw(6) << i << " residu: " << std::scientific << residu << '\n';
+    return 0;
+  }
+
+  static void
+  petscSolveLocalSystem(const CRSMatrix<double, size_t>& A,
+                        Vector<double>& x,
+                        const Vector<double>& b,
+                        const LinearSolverOptions& options)
+  {
+    Assert(x.size() == b.size() and x.size() == A.numberOfRows());
+
+    Vec petscB;
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB);
+    Vec petscX;
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX);
+
+    Mat petscMat;
+    MatCreate(PETSC_COMM_WORLD, &petscMat);
+    MatSetSizes(petscMat, PETSC_DECIDE, PETSC_DECIDE, x.size(), x.size());
+
+    MatSetType(petscMat, MATAIJ);
+
+    Array<PetscScalar> values = copy(A.values());
+    auto hm                   = A.hostMatrix();
+
+    Array<PetscInt> row_indices{hm.row_map.size()};
+    for (size_t i = 0; i < hm.row_map.size(); ++i) {
+      row_indices[i] = hm.row_map[i];
+    }
+
+    Array<PetscInt> column_indices{values.size()};
+    size_t l = 0;
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      const auto row_i = hm.rowConst(i);
+      for (size_t j = 0; j < row_i.length; ++j) {
+        column_indices[l++] = row_i.colidx(j);
+      }
+    }
+
+    MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, x.size(), x.size(), &row_indices[0], &column_indices[0], &values[0],
+                              &petscMat);
+
+    MatAssemblyBegin(petscMat, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(petscMat, MAT_FINAL_ASSEMBLY);
+
+    KSP ksp;
+    KSPCreate(PETSC_COMM_WORLD, &ksp);
+    KSPSetTolerances(ksp, options.epsilon(), 1E-100, 1E5, options.maximumIteration());
+
+    KSPSetOperators(ksp, petscMat, petscMat);
+
+    PC pc;
+    KSPGetPC(ksp, &pc);
+
+    bool direct_solver = false;
+
+    switch (options.method()) {
+    case LSMethod::bicgstab: {
+      KSPSetType(ksp, KSPBCGS);
+      break;
+    }
+    case LSMethod::bicgstab2: {
+      KSPSetType(ksp, KSPBCGSL);
+      KSPBCGSLSetEll(ksp, 2);
+      break;
+    }
+    case LSMethod::cg: {
+      KSPSetType(ksp, KSPCG);
+      break;
+    }
+    case LSMethod::gmres: {
+      KSPSetType(ksp, KSPGMRES);
+
+      break;
+    }
+    case LSMethod::lu: {
+      KSPSetType(ksp, KSPPREONLY);
+      PCSetType(pc, PCLU);
+      PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
+      direct_solver = true;
+      break;
+    }
+    case LSMethod::choleski: {
+      KSPSetType(ksp, KSPPREONLY);
+      PCSetType(pc, PCCHOLESKY);
+      direct_solver = true;
+      break;
+    }
+    default: {
+      throw UnexpectedError("unexpected method: " + name(options.method()));
+    }
+    }
+
+    if (not direct_solver) {
+      switch (options.precond()) {
+      case LSPrecond::diagonal: {
+        PCSetType(pc, PCJACOBI);
+        break;
+      }
+      case LSPrecond::incomplete_LU: {
+        PCSetType(pc, PCILU);
+        break;
+      }
+      case LSPrecond::incomplete_choleski: {
+        PCSetType(pc, PCICC);
+        break;
+      }
+      case LSPrecond::none: {
+        PCSetType(pc, PCNONE);
+        break;
+      }
+      default: {
+        throw UnexpectedError("unexpected preconditioner: " + name(options.precond()));
+      }
+      }
+    }
+    if (options.verbose()) {
+      KSPMonitorSet(ksp, petscMonitor, 0, 0);
+    }
+
+    KSPSolve(ksp, petscB, petscX);
+    MatDestroy(&petscMat);
+  }
+
+#else
+  static void
+  petscSolveLocalSystem(const CRSMatrix<double, size_t>&,
+                        Vector<double>&,
+                        const Vector<double>&,
+                        const LinearSolverOptions&)
+  {
+    checkHasLibrary(LSLibrary::petsc);
+    throw UnexpectedError("unexpected situation should not reach this point!");
+  }
+
+#endif   // PUGS_HAS_PETSC
+};
+
+bool
+LinearSolver::hasLibrary(LSLibrary library) const
+{
+  return Internals::hasLibrary(library);
+}
+
+void
+LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b)
+{
+  switch (m_options.library()) {
+  case LSLibrary::builtin: {
+    Internals::builtinSolveLocalSystem(A, x, b, m_options);
+    break;
+  }
+  case LSLibrary::petsc: {
+    Internals::petscSolveLocalSystem(A, x, b, m_options);
+    break;
+  }
+  default: {
+    throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems for sparse matrices");
+  }
+  }
+}
+
+LinearSolver::LinearSolver(const LinearSolverOptions& options) : m_options{options}
+{
+  Internals::checkHasLibrary(m_options.library());
+  Internals::checkOptions(options);
+}
+
+LinearSolver::LinearSolver() : LinearSolver{LinearSolverOptions::default_options} {}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b72b09b0f6a3255f914c96ee3df7717923931294
--- /dev/null
+++ b/src/algebra/LinearSolver.hpp
@@ -0,0 +1,36 @@
+#ifndef LINEAR_SOLVER_HPP
+#define LINEAR_SOLVER_HPP
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LinearSolverOptions.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <utils/Exceptions.hpp>
+
+class LinearSolver
+{
+ private:
+  struct Internals;
+
+  const LinearSolverOptions m_options;
+
+  void _solveLocalDense(size_t N, const double* A, double* x, const double* b);
+
+ public:
+  bool hasLibrary(LSLibrary library) const;
+
+  void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b);
+
+  template <size_t N>
+  void
+  solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
+  {
+    this->_solveLocalDense(N, &A(0, 0), &x[0], &b[0]);
+  }
+
+  LinearSolver();
+  LinearSolver(const LinearSolverOptions& options);
+};
+
+#endif   // LINEAR_SOLVER_HPP
diff --git a/src/algebra/LinearSolverOptions.cpp b/src/algebra/LinearSolverOptions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b40c5fa92f3af73015afa133eaae99b91feb6fda
--- /dev/null
+++ b/src/algebra/LinearSolverOptions.cpp
@@ -0,0 +1,15 @@
+#include <algebra/LinearSolverOptions.hpp>
+
+#include <rang.hpp>
+
+std::ostream&
+operator<<(std::ostream& os, const LinearSolverOptions& options)
+{
+  os << "  library: " << rang::style::bold << name(options.library()) << rang::style::reset << '\n';
+  os << "  method : " << rang::style::bold << name(options.method()) << rang::style::reset << '\n';
+  os << "  precond: " << rang::style::bold << name(options.precond()) << rang::style::reset << '\n';
+  os << "  epsilon: " << rang::style::bold << options.epsilon() << rang::style::reset << '\n';
+  os << "  maxiter: " << rang::style::bold << options.maximumIteration() << rang::style::reset << '\n';
+  os << "  verbose: " << rang::style::bold << std::boolalpha << options.verbose() << rang::style::reset << '\n';
+  return os;
+}
diff --git a/src/algebra/LinearSolverOptions.hpp b/src/algebra/LinearSolverOptions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd7538b53df8740ee6c151e0d2e72a106fdc3c8c
--- /dev/null
+++ b/src/algebra/LinearSolverOptions.hpp
@@ -0,0 +1,233 @@
+#ifndef LINEAR_SOLVER_OPTIONS_HPP
+#define LINEAR_SOLVER_OPTIONS_HPP
+
+#include <utils/Exceptions.hpp>
+
+#include <iostream>
+
+enum class LSLibrary : int8_t
+{
+  LS__begin = 0,
+  //
+  builtin = LS__begin,
+  petsc,
+  //
+  LS__end
+};
+
+enum class LSMethod : int8_t
+{
+  LS__begin = 0,
+  //
+  cg = LS__begin,
+  bicgstab,
+  bicgstab2,
+  gmres,
+  lu,
+  choleski,
+  //
+  LS__end
+};
+
+enum class LSPrecond : int8_t
+{
+  LS__begin = 0,
+  //
+  none = LS__begin,
+  diagonal,
+  incomplete_choleski,
+  incomplete_LU,
+  //
+  LS__end
+};
+
+inline std::string
+name(const LSLibrary library)
+{
+  switch (library) {
+  case LSLibrary::builtin: {
+    return "builtin";
+  }
+  case LSLibrary::petsc: {
+    return "PETSc";
+  }
+  case LSLibrary::LS__end: {
+  }
+  }
+  throw UnexpectedError("Linear system library name is not defined!");
+}
+
+inline std::string
+name(const LSMethod method)
+{
+  switch (method) {
+  case LSMethod::cg: {
+    return "CG";
+  }
+  case LSMethod::bicgstab: {
+    return "BICGStab";
+  }
+  case LSMethod::bicgstab2: {
+    return "BICGStab2";
+  }
+  case LSMethod::gmres: {
+    return "GMRES";
+  }
+  case LSMethod::lu: {
+    return "LU";
+  }
+  case LSMethod::choleski: {
+    return "Choleski";
+  }
+  case LSMethod::LS__end: {
+  }
+  }
+  throw UnexpectedError("Linear system method name is not defined!");
+}
+
+inline std::string
+name(const LSPrecond precond)
+{
+  switch (precond) {
+  case LSPrecond::none: {
+    return "none";
+  }
+  case LSPrecond::diagonal: {
+    return "diagonal";
+  }
+  case LSPrecond::incomplete_choleski: {
+    return "ICholeski";
+  }
+  case LSPrecond::incomplete_LU: {
+    return "ILU";
+  }
+  case LSPrecond::LS__end: {
+  }
+  }
+  throw UnexpectedError("Linear system preconditioner name is not defined!");
+}
+
+template <typename LSEnumType>
+inline LSEnumType
+getLSEnumFromName(const std::string& enum_name)
+{
+  using BaseT = std::underlying_type_t<LSEnumType>;
+  for (BaseT enum_value = static_cast<BaseT>(LSEnumType::LS__begin);
+       enum_value < static_cast<BaseT>(LSEnumType::LS__end); ++enum_value) {
+    if (name(LSEnumType{enum_value}) == enum_name) {
+      return LSEnumType{enum_value};
+    }
+  }
+  throw NormalError(std::string{"could not find '"} + enum_name + "' associate type!");
+}
+
+template <typename LSEnumType>
+inline void
+printLSEnumListNames(std::ostream& os)
+{
+  using BaseT = std::underlying_type_t<LSEnumType>;
+  for (BaseT enum_value = static_cast<BaseT>(LSEnumType::LS__begin);
+       enum_value < static_cast<BaseT>(LSEnumType::LS__end); ++enum_value) {
+    os << "  - " << name(LSEnumType{enum_value}) << '\n';
+  }
+}
+
+class LinearSolverOptions
+{
+ private:
+  LSLibrary m_library = LSLibrary::builtin;
+  LSMethod m_method   = LSMethod::bicgstab;
+  LSPrecond m_precond = LSPrecond::none;
+
+  double m_epsilon           = 1E-6;
+  size_t m_maximum_iteration = 200;
+
+  bool m_verbose = false;
+
+ public:
+  static LinearSolverOptions default_options;
+
+  friend std::ostream& operator<<(std::ostream& os, const LinearSolverOptions& options);
+
+  LSLibrary&
+  library()
+  {
+    return m_library;
+  }
+
+  LSLibrary
+  library() const
+  {
+    return m_library;
+  }
+
+  LSMethod
+  method() const
+  {
+    return m_method;
+  }
+
+  LSMethod&
+  method()
+  {
+    return m_method;
+  }
+
+  LSPrecond
+  precond() const
+  {
+    return m_precond;
+  }
+
+  LSPrecond&
+  precond()
+  {
+    return m_precond;
+  }
+
+  double
+  epsilon() const
+  {
+    return m_epsilon;
+  }
+
+  double&
+  epsilon()
+  {
+    return m_epsilon;
+  }
+
+  size_t&
+  maximumIteration()
+  {
+    return m_maximum_iteration;
+  }
+
+  size_t
+  maximumIteration() const
+  {
+    return m_maximum_iteration;
+  }
+
+  bool&
+  verbose()
+  {
+    return m_verbose;
+  };
+
+  bool
+  verbose() const
+  {
+    return m_verbose;
+  };
+
+  LinearSolverOptions(const LinearSolverOptions&) = default;
+  LinearSolverOptions(LinearSolverOptions&&)      = default;
+
+  LinearSolverOptions()  = default;
+  ~LinearSolverOptions() = default;
+};
+
+inline LinearSolverOptions LinearSolverOptions::default_options;
+
+#endif   // LINEAR_SOLVER_OPTIONS_HPP
diff --git a/src/algebra/PETScWrapper.cpp b/src/algebra/PETScWrapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd11dc9772f11ad29ca4970c5b88656c623a3508
--- /dev/null
+++ b/src/algebra/PETScWrapper.cpp
@@ -0,0 +1,27 @@
+#include <algebra/PETScWrapper.hpp>
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
+namespace PETScWrapper
+{
+void
+initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
+{
+#ifdef PUGS_HAS_PETSC
+  PetscOptionsSetValue(NULL, "-no_signal_handler", "true");
+  PetscInitialize(&argc, &argv, 0, 0);
+#endif   // PUGS_HAS_PETSC
+}
+
+void
+finalize()
+{
+#ifdef PUGS_HAS_PETSC
+  PetscFinalize();
+#endif   // PUGS_HAS_PETSC
+}
+}   // namespace PETScWrapper
diff --git a/src/algebra/PETScWrapper.hpp b/src/algebra/PETScWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2f40b37f8c4562b7abeccffeb6e8cfb27688c7c4
--- /dev/null
+++ b/src/algebra/PETScWrapper.hpp
@@ -0,0 +1,10 @@
+#ifndef PETSC_WRAPPER_HPP
+#define PETSC_WRAPPER_HPP
+
+namespace PETScWrapper
+{
+void initialize(int& argc, char* argv[]);
+void finalize();
+}   // namespace PETScWrapper
+
+#endif   // PETSC_WRAPPER_HPP
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index 100000366d3bd0e5a072084864e82cdb4b73ba29..d5c14713bb6b91235b796a56f40b9a2f5a605805 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -22,6 +22,12 @@ class SparseMatrixDescriptor
     std::map<IndexType, DataType> m_id_value_map;
 
    public:
+    const auto&
+    columnIdToValueMap() const
+    {
+      return m_id_value_map;
+    }
+
     IndexType
     numberOfValues() const noexcept
     {
@@ -145,7 +151,12 @@ class SparseMatrixDescriptor
     return values;
   }
 
-  SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row} {}
+  SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row}
+  {
+    for (size_t i = 0; i < nb_row; ++i) {
+      m_row_array[i][i] = 0;
+    }
+  }
 
   ~SparseMatrixDescriptor() = default;
 };
diff --git a/src/language/ast/ASTNodeDataType.cpp b/src/language/ast/ASTNodeDataType.cpp
index 10ec58b451ead8402a328744cb977c29557dbbbf..ce2bbcc60182e0fe025d4836ee43fe0cb6d37a59 100644
--- a/src/language/ast/ASTNodeDataType.cpp
+++ b/src/language/ast/ASTNodeDataType.cpp
@@ -48,11 +48,15 @@ dataTypeName(const ASTNodeDataType& data_type)
   case ASTNodeDataType::list_t: {
     std::ostringstream data_type_name_list;
     const auto& data_type_list = data_type.contentTypeList();
-    data_type_name_list << dataTypeName(*data_type_list[0]);
-    for (size_t i = 1; i < data_type_list.size(); ++i) {
-      data_type_name_list << '*' << dataTypeName(*data_type_list[i]);
+    if (data_type_list.size() > 0) {
+      data_type_name_list << dataTypeName(*data_type_list[0]);
+      for (size_t i = 1; i < data_type_list.size(); ++i) {
+        data_type_name_list << '*' << dataTypeName(*data_type_list[i]);
+      }
+      name = "list(" + data_type_name_list.str() + ")";
+    } else {
+      name = "list(void)";
     }
-    name = "list(" + data_type_name_list.str() + ")";
     break;
   }
   case ASTNodeDataType::string_t:
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index f8cf2f0b445882a672d54e60dbc75858bbb246c7..0a996fd6764b4d2eed57c02e8682283f1e640901 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -2,6 +2,7 @@
 
 add_library(PugsLanguageModules
   BuiltinModule.cpp
+  LinearSolverModule.cpp
   MathModule.cpp
   MeshModule.cpp
   ModuleRepository.cpp
diff --git a/src/language/modules/LinearSolverModule.cpp b/src/language/modules/LinearSolverModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9e626ba2b0e4707df0495de2057a00d81ea4b118
--- /dev/null
+++ b/src/language/modules/LinearSolverModule.cpp
@@ -0,0 +1,85 @@
+#include <language/modules/LinearSolverModule.hpp>
+
+#include <algebra/LinearSolver.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <language/utils/TypeDescriptor.hpp>
+
+LinearSolverModule::LinearSolverModule()
+{
+  this->_addBuiltinFunction("setLSVerbosity", std::make_shared<BuiltinFunctionEmbedder<void(const bool&)>>(
+
+                                                [](const bool& verbose) -> void {
+                                                  LinearSolverOptions::default_options.verbose() = verbose;
+                                                }
+
+                                                ));
+
+  this->_addBuiltinFunction("setLSEpsilon", std::make_shared<BuiltinFunctionEmbedder<void(const double&)>>(
+
+                                              [](const double& epsilon) -> void {
+                                                LinearSolverOptions::default_options.epsilon() = epsilon;
+                                              }
+
+                                              ));
+
+  this->_addBuiltinFunction("setLSMaxIter", std::make_shared<BuiltinFunctionEmbedder<void(const uint64_t&)>>(
+
+                                              [](const uint64_t& max_iter) -> void {
+                                                LinearSolverOptions::default_options.maximumIteration() = max_iter;
+                                              }
+
+                                              ));
+
+  this->_addBuiltinFunction("setLSLibrary", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+                                              [](const std::string& library_name) -> void {
+                                                LinearSolverOptions::default_options.library() =
+                                                  getLSEnumFromName<LSLibrary>(library_name);
+                                              }
+
+                                              ));
+
+  this->_addBuiltinFunction("setLSMethod", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+                                             [](const std::string& method_name) -> void {
+                                               LinearSolverOptions::default_options.method() =
+                                                 getLSEnumFromName<LSMethod>(method_name);
+                                             }
+
+                                             ));
+
+  this->_addBuiltinFunction("setLSPrecond", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+                                              [](const std::string& precond_name) -> void {
+                                                LinearSolverOptions::default_options.precond() =
+                                                  getLSEnumFromName<LSPrecond>(precond_name);
+                                              }
+
+                                              ));
+
+  this->_addBuiltinFunction("printLSOptions", std::make_shared<BuiltinFunctionEmbedder<void()>>(
+
+                                                []() -> void {
+                                                  std::cout << rang::fgB::yellow << "Linear solver options"
+                                                            << rang::style::reset << '\n';
+                                                  std::cout << LinearSolverOptions::default_options;
+                                                }
+
+                                                ));
+
+  this->_addBuiltinFunction("printLSAvailable",
+                            std::make_shared<BuiltinFunctionEmbedder<void()>>(
+
+                              []() -> void {
+                                std::cout << rang::fgB::yellow << "Available linear solver options"
+                                          << rang::style::reset << '\n';
+                                std::cout << rang::fgB::blue << " libraries" << rang::style::reset << '\n';
+                                printLSEnumListNames<LSLibrary>(std::cout);
+                                std::cout << rang::fgB::blue << " methods" << rang::style::reset << '\n';
+                                printLSEnumListNames<LSMethod>(std::cout);
+                                std::cout << rang::fgB::blue << " preconditioners" << rang::style::reset << '\n';
+                                printLSEnumListNames<LSPrecond>(std::cout);
+                              }
+
+                              ));
+}
diff --git a/src/language/modules/LinearSolverModule.hpp b/src/language/modules/LinearSolverModule.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5dc7ae64efd3ed43f554e821693f761b0a66d412
--- /dev/null
+++ b/src/language/modules/LinearSolverModule.hpp
@@ -0,0 +1,19 @@
+#ifndef LINEAR_SOLVER_MODULE_HPP
+#define LINEAR_SOLVER_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+
+class LinearSolverModule : public BuiltinModule
+{
+ public:
+  std::string_view
+  name() const final
+  {
+    return "linear_solver";
+  }
+
+  LinearSolverModule();
+  ~LinearSolverModule() = default;
+};
+
+#endif   // LINEAR_SOLVER_MODULE_HPP
diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp
index 497ba083ac2d9d2ca0359c5d73e85e094b1fed0e..df403fa87f0c0d39778414a92ec7aa991b3d8d5b 100644
--- a/src/language/modules/ModuleRepository.cpp
+++ b/src/language/modules/ModuleRepository.cpp
@@ -1,6 +1,7 @@
 #include <language/modules/ModuleRepository.hpp>
 
 #include <language/ast/ASTNode.hpp>
+#include <language/modules/LinearSolverModule.hpp>
 #include <language/modules/MathModule.hpp>
 #include <language/modules/MeshModule.hpp>
 #include <language/modules/SchemeModule.hpp>
@@ -23,6 +24,7 @@ ModuleRepository::ModuleRepository()
   this->_subscribe(std::make_unique<MeshModule>());
   this->_subscribe(std::make_unique<VTKModule>());
   this->_subscribe(std::make_unique<SchemeModule>());
+  this->_subscribe(std::make_unique<LinearSolverModule>());
 }
 
 template <typename NameEmbedderMapT, typename EmbedderTableT>
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 13b8633c2664add44d91b11b0a6d55afe79d2955..7162aeccf69a7b5e6420e23ae7dca6cea6ced9b1 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -8,6 +8,10 @@
 #include <mpi.h>
 #endif   //  PUGS_HAS_MPI
 
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
 std::string
 BuildInfo::type()
 {
@@ -42,3 +46,14 @@ BuildInfo::mpiLibrary()
   return "none";
 #endif   // PUGS_HAS_MPI
 }
+
+std::string
+BuildInfo::petscLibrary()
+{
+#ifdef PUGS_HAS_PETSC
+  return std::to_string(PETSC_VERSION_MAJOR) + "." + std::to_string(PETSC_VERSION_MINOR) + "." +
+         std::to_string(PETSC_VERSION_SUBMINOR);
+#else
+  return "none";
+#endif   // PUGS_HAS_PETSC
+}
diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp
index 81d2b4960085fbeffff667fb92e8a4069848fb8d..67134a782f595ccf68dc89cd29a6f273c7ce6140 100644
--- a/src/utils/BuildInfo.hpp
+++ b/src/utils/BuildInfo.hpp
@@ -9,6 +9,7 @@ struct BuildInfo
   static std::string compiler();
   static std::string kokkosDevices();
   static std::string mpiLibrary();
+  static std::string petscLibrary();
 };
 
 #endif   // BUILD_INFO_HPP
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 7ea50b0ed95845d2cee2633477e33080453f0f30..8f388c5b492dcd8065ed77740d322556bdc76db3 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -9,6 +9,8 @@
 #include <utils/FPEManager.hpp>
 #include <utils/SignalManager.hpp>
 
+#include <algebra/PETScWrapper.hpp>
+
 #include <rang.hpp>
 
 #include <Kokkos_Core.hpp>
@@ -41,7 +43,8 @@ initialize(int& argc, char* argv[])
   std::cout << "type:     " << rang::style::bold << BuildInfo::type() << rang::style::reset << '\n';
   std::cout << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
   std::cout << "kokkos:   " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
-  std::cout << "mpi:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+  std::cout << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+  std::cout << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
   std::cout << "-------------------------------------------------------\n";
 
   std::string filename;
@@ -81,6 +84,8 @@ initialize(int& argc, char* argv[])
     SignalManager::init(enable_signals);
   }
 
+  PETScWrapper::initialize(argc, argv);
+
   Kokkos::initialize(argc, argv);
   std::cout << "-------------------- " << rang::fg::green << "exec info" << rang::fg::reset
             << " ------------------------" << '\n';
@@ -97,5 +102,6 @@ void
 finalize()
 {
   Kokkos::finalize();
+  PETScWrapper::finalize();
   parallel::Messenger::destroy();
 }
diff --git a/src/utils/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in
index 45878e64327797cd6549b18f4d6922f8cde8a29d..326884f97f401f0d70688fa4ad7259bdb64ab967 100644
--- a/src/utils/pugs_config.hpp.in
+++ b/src/utils/pugs_config.hpp.in
@@ -3,6 +3,7 @@
 
 #cmakedefine PUGS_HAS_FENV_H
 #cmakedefine PUGS_HAS_MPI
+#cmakedefine PUGS_HAS_PETSC
 
 #cmakedefine SYSTEM_IS_LINUX
 #cmakedefine SYSTEM_IS_DARWIN
@@ -10,4 +11,4 @@
 
 #define PUGS_BUILD_TYPE "@CMAKE_BUILD_TYPE@"
 
-#endif // PUGS_CONFIG_HPP
+#endif   // PUGS_CONFIG_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e5c35aab5062d7da7a1be90159d8440b2a6f8370..dbaf666c108fe6d96830c53bdb7c87aa50a912c2 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -48,7 +48,7 @@ add_executable (unit_tests
   test_BuiltinFunctionEmbedder.cpp
   test_BuiltinFunctionEmbedderTable.cpp
   test_BuiltinFunctionProcessor.cpp
-  test_MathModule.cpp
+  test_CG.cpp
   test_ContinueProcessor.cpp
   test_ConcatExpressionProcessor.cpp
   test_CRSMatrix.cpp
@@ -67,10 +67,10 @@ add_executable (unit_tests
   test_INodeProcessor.cpp
   test_ItemType.cpp
   test_ListAffectationProcessor.cpp
+  test_MathModule.cpp
   test_NameProcessor.cpp
   test_OStreamProcessor.cpp
   test_ParseError.cpp
-  test_PCG.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsAssert.cpp
   test_RevisionInfo.cpp
@@ -96,19 +96,23 @@ target_link_libraries (unit_tests
   PugsLanguageAlgorithms
   PugsLanguageUtils
   PugsMesh
+  PugsAlgebra
   PugsUtils
   kokkos
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+  ${PETSC_LIBRARIES}
   Catch2
   )
 
 target_link_libraries (mpi_unit_tests
   PugsUtils
+  PugsAlgebra
   PugsMesh
   kokkos
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+  ${PETSC_LIBRARIES}
   Catch2
   )
 
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index 99b94515ba7795903434db761d509b9052fcaf71..49a6955be4ae1da7c45e036977a3036be0ccf6d3 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    BiCGStab<false>{b, A, x, 10, 1e-12};
+    BiCGStab{A, x, b, 1e-12, 10, false};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
   }
@@ -88,7 +88,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    BiCGStab{b, A, x, 1, 1e-12};
+    BiCGStab{A, x, b, 1e-12, 1, true};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) > 1E-5 * std::sqrt((x, x)));
   }
diff --git a/tests/test_PCG.cpp b/tests/test_CG.cpp
similarity index 91%
rename from tests/test_PCG.cpp
rename to tests/test_CG.cpp
index d4564761d8c564082e4668360e5a2e2a8ecce008..2ce9febb0a31b087d2a7651e4442245b9d6c51ce 100644
--- a/tests/test_PCG.cpp
+++ b/tests/test_CG.cpp
@@ -1,11 +1,11 @@
 #include <catch2/catch.hpp>
 
+#include <algebra/CG.hpp>
 #include <algebra/CRSMatrix.hpp>
-#include <algebra/PCG.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
-TEST_CASE("PCG", "[algebra]")
+TEST_CASE("CG", "[algebra]")
 {
   SECTION("no preconditionner")
   {
@@ -45,7 +45,7 @@ TEST_CASE("PCG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    PCG{b, A, A, x, 10, 1e-12};
+    CG{A, x, b, 1e-12, 10};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
   }
@@ -62,7 +62,7 @@ TEST_CASE("PCG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    PCG<false>{b, A, A, x, 10, 1e-12};
+    CG{A, x, b, 1e-12, 10, false};
     REQUIRE(std::sqrt((x, x)) == 0);
   }
 
@@ -104,7 +104,7 @@ TEST_CASE("PCG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    PCG<false>{b, A, A, x, 1, 1e-12};
+    CG{A, x, b, 1e-12, 1, false};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
   }
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index cd66380b1d1beb6bcb9ba1a9c0d9c495947a3235..8616a41da0ff8a60857ef314c74e04becf6c651f 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -68,18 +68,18 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
     S(4, 1) = 1;
     S(4, 4) = -2;
 
-    REQUIRE(S.row(0).numberOfValues() == 1);
+    REQUIRE(S.row(0).numberOfValues() == 2);
     REQUIRE(S.row(1).numberOfValues() == 2);
     REQUIRE(S.row(2).numberOfValues() == 1);
-    REQUIRE(S.row(3).numberOfValues() == 1);
+    REQUIRE(S.row(3).numberOfValues() == 2);
     REQUIRE(S.row(4).numberOfValues() == 2);
 
     const auto& const_S = S;
 
-    REQUIRE(const_S.row(0).numberOfValues() == 1);
+    REQUIRE(const_S.row(0).numberOfValues() == 2);
     REQUIRE(const_S.row(1).numberOfValues() == 2);
     REQUIRE(const_S.row(2).numberOfValues() == 1);
-    REQUIRE(const_S.row(3).numberOfValues() == 1);
+    REQUIRE(const_S.row(3).numberOfValues() == 2);
     REQUIRE(const_S.row(4).numberOfValues() == 2);
 
 #ifndef NDEBUG
@@ -126,13 +126,14 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
     const auto graph = S.graphVector();
 
     REQUIRE(graph.size() == S.numberOfRows());
-    REQUIRE(graph[0].size() == 1);
+    REQUIRE(graph[0].size() == 2);
     REQUIRE(graph[1].size() == 2);
     REQUIRE(graph[2].size() == 1);
     REQUIRE(graph[3].size() == 2);
     REQUIRE(graph[4].size() == 2);
 
-    REQUIRE(graph[0][0] == 2);
+    REQUIRE(graph[0][0] == 0);
+    REQUIRE(graph[0][1] == 2);
     REQUIRE(graph[1][0] == 1);
     REQUIRE(graph[1][1] == 2);
     REQUIRE(graph[2][0] == 2);
@@ -157,16 +158,17 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
 
     const auto value_array = S.valueArray();
 
-    REQUIRE(value_array.size() == 8);
-
-    REQUIRE(value_array[0] == 5);
-    REQUIRE(value_array[1] == 1);
-    REQUIRE(value_array[2] == 11);
-    REQUIRE(value_array[3] == 4);
-    REQUIRE(value_array[4] == -3);
-    REQUIRE(value_array[5] == 5);
-    REQUIRE(value_array[6] == 1);
-    REQUIRE(value_array[7] == -2);
+    REQUIRE(value_array.size() == 9);
+
+    REQUIRE(value_array[0] == 0);
+    REQUIRE(value_array[1] == 5);
+    REQUIRE(value_array[2] == 1);
+    REQUIRE(value_array[3] == 11);
+    REQUIRE(value_array[4] == 4);
+    REQUIRE(value_array[5] == -3);
+    REQUIRE(value_array[6] == 5);
+    REQUIRE(value_array[7] == 1);
+    REQUIRE(value_array[8] == -2);
   }
 
   SECTION("output")
@@ -186,7 +188,7 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
     output << '\n' << S;
 
     std::string expected_output = R"(
-0 | 2:5
+0 | 0:0 2:5
 1 | 1:1 2:11
 2 | 2:4
 3 | 1:-3 3:5