diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index d9b03ac93fb2f1c4baeff98c90daac17547ff7cf..98297d759675fdb7fdbf8be1e638f19cba4cc562 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -17,7 +17,7 @@ class CRSMatrix
   using MutableDataType = std::remove_const_t<DataType>;
 
  private:
-  using HostMatrix = Kokkos::StaticCrsGraph<IndexType, Kokkos::HostSpace>;
+  using HostMatrix = Kokkos::StaticCrsGraph<const IndexType, Kokkos::HostSpace>;
 
   HostMatrix m_host_matrix;
   Array<const DataType> m_values;
@@ -37,9 +37,15 @@ class CRSMatrix
   }
 
   auto
-  hostMatrix() const
+  rowIndices() const
   {
-    return m_host_matrix;
+    return encapsulate(m_host_matrix.row_map);
+  }
+
+  auto
+  row(size_t i) const
+  {
+    return m_host_matrix.rowConst(i);
   }
 
   template <typename DataType2>
@@ -70,8 +76,17 @@ class CRSMatrix
 
   CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
   {
-    m_host_matrix = Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", M.graphVector());
-    m_values      = M.valueArray();
+    {
+      auto host_matrix =
+        Kokkos::create_staticcrsgraph<Kokkos::StaticCrsGraph<IndexType, Kokkos::HostSpace>>("connectivity_matrix",
+                                                                                            M.graphVector());
+
+      // This is a bit crappy but it is the price to pay to avoid
+      m_host_matrix.entries           = host_matrix.entries;
+      m_host_matrix.row_map           = host_matrix.row_map;
+      m_host_matrix.row_block_offsets = host_matrix.row_block_offsets;
+    }
+    m_values = M.valueArray();
   }
   ~CRSMatrix() = default;
 };
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 6a6beaa05bd3ac46b1a81f6e8645721003fe5577..1965a9966cafdf181b366feb1e41909077edd8ac 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -24,9 +24,11 @@ struct LinearSolver::Internals
       return false;
 #endif
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("Linear system library (" + ::name(library) + ") was not set!");
     }
+      // LCOV_EXCL_STOP
     }
   }
 
@@ -34,7 +36,9 @@ struct LinearSolver::Internals
   checkHasLibrary(const LSLibrary library)
   {
     if (not hasLibrary(library)) {
+      // LCOV_EXCL_START
       throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!");
+      // LCOV_EXCL_STOP
     }
   }
 
@@ -64,9 +68,11 @@ struct LinearSolver::Internals
     case LSMethod::choleski: {
       break;
     }
+      // LCOV_EXCL_START
     default: {
       throw NormalError(name(method) + " is not a builtin linear solver!");
     }
+      // LCOV_EXCL_STOP
     }
   }
 
@@ -94,9 +100,11 @@ struct LinearSolver::Internals
     case LSPrecond::incomplete_LU: {
       break;
     }
+      // LCOV_EXCL_START
     default: {
       throw NormalError(name(precond) + " is not a PETSc preconditioner!");
     }
+      // LCOV_EXCL_STOP
     }
   }
 
@@ -114,9 +122,11 @@ struct LinearSolver::Internals
       checkPETScPrecond(options.precond());
       break;
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("undefined options compatibility for this library (" + ::name(options.library()) + ")!");
     }
+      // LCOV_EXCL_STOP
     }
   }
 
@@ -127,7 +137,9 @@ struct LinearSolver::Internals
                           const LinearSolverOptions& options)
   {
     if (options.precond() != LSPrecond::none) {
+      // LCOV_EXCL_START
       throw UnexpectedError("builtin linear solver do not allow any preconditioner!");
+      // LCOV_EXCL_STOP
     }
     switch (options.method()) {
     case LSMethod::cg: {
@@ -138,14 +150,15 @@ struct LinearSolver::Internals
       BiCGStab{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
       break;
     }
+      // LCOV_EXCL_START
     default: {
       throw NotImplementedError("undefined builtin method: " + name(options.method()));
     }
+      // LCOV_EXCL_STOP
     }
   }
 
 #ifdef PUGS_HAS_PETSC
-
   static int
   petscMonitor(KSP, int i, double residu, void*)
   {
@@ -173,17 +186,17 @@ struct LinearSolver::Internals
     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];
+    const auto A_row_indices = A.rowIndices();
+    Array<PetscInt> row_indices{A_row_indices.size()};
+    for (size_t i = 0; i < row_indices.size(); ++i) {
+      row_indices[i] = A_row_indices[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);
+      const auto row_i = A.row(i);
       for (size_t j = 0; j < row_i.length; ++j) {
         column_indices[l++] = row_i.colidx(j);
       }
@@ -238,9 +251,11 @@ struct LinearSolver::Internals
       direct_solver = true;
       break;
     }
+      // LCOV_EXCL_START
     default: {
       throw UnexpectedError("unexpected method: " + name(options.method()));
     }
+      // LCOV_EXCL_STOP
     }
 
     if (not direct_solver) {
@@ -265,9 +280,11 @@ struct LinearSolver::Internals
         PCSetType(pc, PCNONE);
         break;
       }
+        // LCOV_EXCL_START
       default: {
         throw UnexpectedError("unexpected preconditioner: " + name(options.precond()));
       }
+        // LCOV_EXCL_STOP
       }
     }
     if (options.verbose()) {
@@ -278,7 +295,9 @@ struct LinearSolver::Internals
     MatDestroy(&petscMat);
   }
 
-#else
+#else   // PUGS_HAS_PETSC
+
+  // LCOV_EXCL_START
   static void
   petscSolveLocalSystem(const CRSMatrix<double, size_t>&,
                         Vector<double>&,
@@ -288,6 +307,7 @@ struct LinearSolver::Internals
     checkHasLibrary(LSLibrary::petsc);
     throw UnexpectedError("unexpected situation should not reach this point!");
   }
+  // LCOV_EXCL_STOP
 
 #endif   // PUGS_HAS_PETSC
 };
@@ -298,6 +318,12 @@ LinearSolver::hasLibrary(LSLibrary library) const
   return Internals::hasLibrary(library);
 }
 
+void
+LinearSolver::checkOptions(const LinearSolverOptions& options) const
+{
+  Internals::checkOptions(options);
+}
+
 void
 LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b)
 {
@@ -306,13 +332,18 @@ LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double
     Internals::builtinSolveLocalSystem(A, x, b, m_options);
     break;
   }
+    // LCOV_EXCL_START
   case LSLibrary::petsc: {
+    // not covered since if PETSc is not linked this point is
+    // unreachable: LinearSolver throws an exception at construction
+    // in this case.
     Internals::petscSolveLocalSystem(A, x, b, m_options);
     break;
   }
   default: {
     throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems for sparse matrices");
   }
+    // LCOV_EXCL_STOP
   }
 }
 
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index b72b09b0f6a3255f914c96ee3df7717923931294..3ed5660bd9ce51c1e5ba5de3fe00c595aa96f1e7 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -19,6 +19,7 @@ class LinearSolver
 
  public:
   bool hasLibrary(LSLibrary library) const;
+  void checkOptions(const LinearSolverOptions& options) const;
 
   void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b);
 
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index d5c14713bb6b91235b796a56f40b9a2f5a605805..f63a17576d619bf692f36a00e5c0c4e373effdbf 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -22,12 +22,6 @@ class SparseMatrixDescriptor
     std::map<IndexType, DataType> m_id_value_map;
 
    public:
-    const auto&
-    columnIdToValueMap() const
-    {
-      return m_id_value_map;
-    }
-
     IndexType
     numberOfValues() const noexcept
     {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a08e1a689cf3041dd62b90c81630381581b83f38..aaf4bfbad904a6466b4442b0d99f35a41840ddb3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -73,6 +73,8 @@ add_executable (unit_tests
   test_IncDecExpressionProcessor.cpp
   test_INodeProcessor.cpp
   test_ItemType.cpp
+  test_LinearSolver.cpp
+  test_LinearSolverOptions.cpp
   test_ListAffectationProcessor.cpp
   test_MathModule.cpp
   test_NameProcessor.cpp
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 45b733f94d6877a5c6792ed7439965a00e2aa428..87953af0c1d98a6b099b8d99606d8f4eaa686a40 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -3,6 +3,7 @@
 
 #include <Kokkos_Core.hpp>
 
+#include <algebra/PETScWrapper.hpp>
 #include <mesh/DiamondDualConnectivityManager.hpp>
 #include <mesh/DiamondDualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -21,6 +22,8 @@ main(int argc, char* argv[])
   parallel::Messenger::create(argc, argv);
   Kokkos::initialize({4, -1, -1, true});
 
+  PETScWrapper::initialize(argc, argv);
+
   const std::string output_base_name{"mpi_test_rank_"};
 
   std::filesystem::path parallel_output(std::string{PUGS_BINARY_DIR});
@@ -86,6 +89,8 @@ main(int argc, char* argv[])
     }
   }
 
+  PETScWrapper::finalize();
+
   Kokkos::finalize();
   parallel::Messenger::destroy();
 
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index 2ce9febb0a31b087d2a7651e4442245b9d6c51ce..75ab83063b9c19121e598ec68a5260ce95b02ae8 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -45,7 +45,7 @@ TEST_CASE("CG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    CG{A, x, b, 1e-12, 10};
+    CG{A, x, b, 1e-12, 10, true};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
   }
@@ -62,7 +62,7 @@ TEST_CASE("CG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    CG{A, x, b, 1e-12, 10, false};
+    CG{A, x, b, 1e-12, 10};
     REQUIRE(std::sqrt((x, x)) == 0);
   }
 
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index 7a7a943d14c943f5a315544fc6eab947ffc28463..6c8584796b8c71991941870c723cc23617c0c2ac 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -94,7 +94,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[4] == -2);
   }
 
-  SECTION("matrix vector product (complet)")
+  SECTION("matrix vector product (complete)")
   {
     SparseMatrixDescriptor<int, uint8_t> S{4};
     S(0, 0) = 1;
@@ -129,6 +129,60 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[3] == 150);
   }
 
+  SECTION("check values")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{4};
+    S(3, 0) = 13;
+    S(0, 0) = 1;
+    S(0, 1) = 2;
+    S(1, 1) = 6;
+    S(1, 2) = 7;
+    S(2, 2) = 11;
+    S(3, 2) = 15;
+    S(2, 0) = 9;
+    S(3, 3) = 16;
+    S(2, 3) = 12;
+    S(0, 3) = 4;
+    S(1, 0) = 5;
+    S(2, 1) = 10;
+
+    CRSMatrix<int, uint8_t> A{S};
+
+    auto values = A.values();
+    REQUIRE(values.size() == 13);
+    REQUIRE(values[0] == 1);
+    REQUIRE(values[1] == 2);
+    REQUIRE(values[2] == 4);
+    REQUIRE(values[3] == 5);
+    REQUIRE(values[4] == 6);
+    REQUIRE(values[5] == 7);
+    REQUIRE(values[6] == 9);
+    REQUIRE(values[7] == 10);
+    REQUIRE(values[8] == 11);
+    REQUIRE(values[9] == 12);
+    REQUIRE(values[10] == 13);
+    REQUIRE(values[11] == 15);
+    REQUIRE(values[12] == 16);
+
+    auto row_indices = A.rowIndices();
+
+    REQUIRE(row_indices.size() == 5);
+
+    REQUIRE(A.row(0).colidx(0) == 0);
+    REQUIRE(A.row(0).colidx(1) == 1);
+    REQUIRE(A.row(0).colidx(2) == 3);
+    REQUIRE(A.row(1).colidx(0) == 0);
+    REQUIRE(A.row(1).colidx(1) == 1);
+    REQUIRE(A.row(1).colidx(2) == 2);
+    REQUIRE(A.row(2).colidx(0) == 0);
+    REQUIRE(A.row(2).colidx(1) == 1);
+    REQUIRE(A.row(2).colidx(2) == 2);
+    REQUIRE(A.row(2).colidx(3) == 3);
+    REQUIRE(A.row(3).colidx(0) == 0);
+    REQUIRE(A.row(3).colidx(1) == 2);
+    REQUIRE(A.row(3).colidx(2) == 3);
+  }
+
 #ifndef NDEBUG
   SECTION("incompatible runtime matrix/vector product")
   {
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f56bb59ed1d6b7df2e81d8e4076f0ee23bb27ec6
--- /dev/null
+++ b/tests/test_LinearSolver.cpp
@@ -0,0 +1,527 @@
+#include <catch2/catch.hpp>
+
+#include <utils/pugs_config.hpp>
+
+#include <algebra/LinearSolver.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LinearSolver", "[algebra]")
+{
+  SECTION("check has library")
+  {
+    LinearSolver linear_solver;
+
+    REQUIRE(linear_solver.hasLibrary(LSLibrary::builtin) == true);
+
+#ifdef PUGS_HAS_PETSC
+    REQUIRE(linear_solver.hasLibrary(LSLibrary::petsc) == true);
+#else    // PUGS_HAS_PETSC
+    REQUIRE(linear_solver.hasLibrary(LSLibrary::petsc) == false);
+#endif   // PUGS_HAS_PETSC
+  }
+
+  SECTION("check linear solver building")
+  {
+    LinearSolverOptions options;
+
+    SECTION("builtin")
+    {
+      SECTION("builtin methods")
+      {
+        options.library() = LSLibrary::builtin;
+        options.precond() = LSPrecond::none;
+
+        options.method() = LSMethod::cg;
+        REQUIRE_NOTHROW(LinearSolver{options});
+
+        options.method() = LSMethod::bicgstab;
+        REQUIRE_NOTHROW(LinearSolver{options});
+
+        options.method() = LSMethod::bicgstab2;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not a builtin linear solver!");
+
+        options.method() = LSMethod::lu;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: LU is not a builtin linear solver!");
+
+        options.method() = LSMethod::choleski;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Choleski is not a builtin linear solver!");
+
+        options.method() = LSMethod::gmres;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: GMRES is not a builtin linear solver!");
+      }
+
+      SECTION("builtin precond")
+      {
+        options.library() = LSLibrary::builtin;
+        options.method()  = LSMethod::cg;
+
+        options.precond() = LSPrecond::none;
+        REQUIRE_NOTHROW(LinearSolver{options});
+
+        options.precond() = LSPrecond::diagonal;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: diagonal is not a builtin preconditioner!");
+
+        options.precond() = LSPrecond::incomplete_LU;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ILU is not a builtin preconditioner!");
+
+        options.precond() = LSPrecond::incomplete_choleski;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ICholeski is not a builtin preconditioner!");
+
+        options.precond() = LSPrecond::amg;
+        REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not a builtin preconditioner!");
+      }
+    }
+
+    SECTION("PETSc")
+    {
+      LinearSolverOptions always_valid;
+      always_valid.library() = LSLibrary::builtin;
+      always_valid.method()  = LSMethod::cg;
+      always_valid.precond() = LSPrecond::none;
+
+      LinearSolver linear_solver{always_valid};
+
+      SECTION("PETSc methods")
+      {
+        options.library() = LSLibrary::petsc;
+        options.precond() = LSPrecond::none;
+
+        options.method() = LSMethod::cg;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::bicgstab;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::bicgstab2;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::lu;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::choleski;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::gmres;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+      }
+
+      SECTION("builtin precond")
+      {
+        options.library() = LSLibrary::petsc;
+        options.method()  = LSMethod::cg;
+
+        options.precond() = LSPrecond::none;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::diagonal;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::incomplete_LU;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::incomplete_choleski;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::amg;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+      }
+    }
+
+#ifndef PUGS_HAS_PETSC
+    SECTION("not linked PETSc")
+    {
+      options.library() = LSLibrary::petsc;
+      options.method()  = LSMethod::cg;
+      options.precond() = LSPrecond::none;
+
+      REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+    }
+#endif   // PUGS_HAS_PETSC
+  }
+
+  SECTION("check linear solvers")
+  {
+    SECTION("symmetric system")
+    {
+      SparseMatrixDescriptor S{5};
+      S(0, 0) = 2;
+      S(0, 1) = -1;
+
+      S(1, 0) = -1;
+      S(1, 1) = 2;
+      S(1, 2) = -1;
+
+      S(2, 1) = -1;
+      S(2, 2) = 2;
+      S(2, 3) = -1;
+
+      S(3, 2) = -1;
+      S(3, 3) = 2;
+      S(3, 4) = -1;
+
+      S(4, 3) = -1;
+      S(4, 4) = 2;
+
+      CRSMatrix A{S};
+
+      Vector<const double> x_exact = [] {
+        Vector<double> y{5};
+        y[0] = 1;
+        y[1] = 3;
+        y[2] = 2;
+        y[3] = 4;
+        y[4] = 5;
+        return y;
+      }();
+
+      Vector<double> b = A * x_exact;
+
+      SECTION("builtin")
+      {
+        SECTION("CG no preconditioner")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::builtin;
+          options.method()  = LSMethod::cg;
+          options.precond() = LSPrecond::none;
+
+          Vector<double> x{5};
+          x = 0;
+
+          LinearSolver solver{options};
+
+          solver.solveLocalSystem(A, x, b);
+          Vector error = x - x_exact;
+          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+        }
+      }
+
+      SECTION("PETSc")
+      {
+#ifdef PUGS_HAS_PETSC
+
+        SECTION("CG")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::cg;
+          options.precond() = LSPrecond::none;
+          options.verbose() = true;
+
+          SECTION("CG no preconditioner")
+          {
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("CG Diagonal")
+          {
+            options.precond() = LSPrecond::diagonal;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("CG ICholeski")
+          {
+            options.precond() = LSPrecond::incomplete_choleski;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("CG AMG")
+          {
+            options.precond() = LSPrecond::amg;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+        }
+
+        SECTION("Choleski")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::choleski;
+          options.precond() = LSPrecond::none;
+
+          Vector<double> x{5};
+          x = 0;
+
+          LinearSolver solver{options};
+
+          solver.solveLocalSystem(A, x, b);
+          Vector error = x - x_exact;
+          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+        }
+
+#else    // PUGS_HAS_PETSC
+        SECTION("PETSc not linked")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::cg;
+          options.precond() = LSPrecond::none;
+
+          REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+        }
+#endif   // PUGS_HAS_PETSC
+      }
+    }
+
+    SECTION("none symmetric system")
+    {
+      SparseMatrixDescriptor S{5};
+      S(0, 0) = 2;
+      S(0, 1) = -1;
+
+      S(1, 0) = -0.2;
+      S(1, 1) = 2;
+      S(1, 2) = -1;
+
+      S(2, 1) = -1;
+      S(2, 2) = 4;
+      S(2, 3) = -2;
+
+      S(3, 2) = -1;
+      S(3, 3) = 2;
+      S(3, 4) = -0.1;
+
+      S(4, 3) = 1;
+      S(4, 4) = 3;
+
+      CRSMatrix A{S};
+
+      Vector<const double> x_exact = [] {
+        Vector<double> y{5};
+        y[0] = 1;
+        y[1] = 3;
+        y[2] = 2;
+        y[3] = 4;
+        y[4] = 5;
+        return y;
+      }();
+
+      Vector<double> b = A * x_exact;
+
+      SECTION("builtin")
+      {
+        SECTION("BICGStab no preconditioner")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::builtin;
+          options.method()  = LSMethod::bicgstab;
+          options.precond() = LSPrecond::none;
+
+          Vector<double> x{5};
+          x = 0;
+
+          LinearSolver solver{options};
+
+          solver.solveLocalSystem(A, x, b);
+          Vector error = x - x_exact;
+          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+        }
+      }
+
+      SECTION("PETSc")
+      {
+#ifdef PUGS_HAS_PETSC
+
+        SECTION("BICGStab")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::bicgstab;
+          options.precond() = LSPrecond::none;
+          options.verbose() = true;
+
+          SECTION("BICGStab no preconditioner")
+          {
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("BICGStab Diagonal")
+          {
+            options.precond() = LSPrecond::diagonal;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("BICGStab ILU")
+          {
+            options.precond() = LSPrecond::incomplete_LU;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+        }
+
+        SECTION("BICGStab2")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::bicgstab2;
+          options.precond() = LSPrecond::none;
+
+          SECTION("BICGStab2 no preconditioner")
+          {
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("BICGStab2 Diagonal")
+          {
+            options.precond() = LSPrecond::diagonal;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+        }
+
+        SECTION("GMRES")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::gmres;
+          options.precond() = LSPrecond::none;
+
+          SECTION("GMRES no preconditioner")
+          {
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("GMRES Diagonal")
+          {
+            options.precond() = LSPrecond::diagonal;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+
+          SECTION("GMRES ILU")
+          {
+            options.precond() = LSPrecond::incomplete_LU;
+
+            Vector<double> x{5};
+            x = 0;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          }
+        }
+
+        SECTION("LU")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::lu;
+          options.precond() = LSPrecond::none;
+
+          Vector<double> x{5};
+          x = 0;
+
+          LinearSolver solver{options};
+
+          solver.solveLocalSystem(A, x, b);
+          Vector error = x - x_exact;
+          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+        }
+
+#else    // PUGS_HAS_PETSC
+        SECTION("PETSc not linked")
+        {
+          LinearSolverOptions options;
+          options.library() = LSLibrary::petsc;
+          options.method()  = LSMethod::cg;
+          options.precond() = LSPrecond::none;
+
+          REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+        }
+#endif   // PUGS_HAS_PETSC
+      }
+    }
+  }
+}
diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..494c36760d7b48c6e40c569423a100fed669d0d6
--- /dev/null
+++ b/tests/test_LinearSolverOptions.cpp
@@ -0,0 +1,177 @@
+#include <catch2/catch.hpp>
+
+#include <algebra/LinearSolverOptions.hpp>
+
+#include <sstream>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LinearSolverOptions", "[algebra]")
+{
+  SECTION("print options")
+  {
+    SECTION("first set")
+    {
+      LinearSolverOptions options;
+      options.verbose()          = true;
+      options.library()          = LSLibrary::builtin;
+      options.method()           = LSMethod::cg;
+      options.precond()          = LSPrecond::incomplete_choleski;
+      options.epsilon()          = 1E-3;
+      options.maximumIteration() = 100;
+
+      std::stringstream os;
+      os << '\n' << options;
+
+      std::stringstream expected_output;
+      expected_output << R"(
+  library: builtin
+  method : CG
+  precond: ICholeski
+  epsilon: )" << 1E-3 << R"(
+  maxiter: 100
+  verbose: true
+)";
+
+      REQUIRE(os.str() == expected_output.str());
+    }
+
+    SECTION("second set")
+    {
+      LinearSolverOptions options;
+
+      options.verbose()          = false;
+      options.library()          = LSLibrary::petsc;
+      options.method()           = LSMethod::gmres;
+      options.precond()          = LSPrecond::incomplete_LU;
+      options.epsilon()          = 1E-6;
+      options.maximumIteration() = 200;
+
+      std::stringstream os;
+      os << '\n' << options;
+
+      std::stringstream expected_output;
+      expected_output << R"(
+  library: PETSc
+  method : GMRES
+  precond: ILU
+  epsilon: )" << 1E-6 << R"(
+  maxiter: 200
+  verbose: false
+)";
+
+      REQUIRE(os.str() == expected_output.str());
+    }
+  }
+
+  SECTION("library name")
+  {
+    REQUIRE(name(LSLibrary::builtin) == "builtin");
+    REQUIRE(name(LSLibrary::petsc) == "PETSc");
+    REQUIRE_THROWS_WITH(name(LSLibrary::LS__end), "unexpected error: Linear system library name is not defined!");
+  }
+
+  SECTION("method name")
+  {
+    REQUIRE(name(LSMethod::cg) == "CG");
+    REQUIRE(name(LSMethod::bicgstab) == "BICGStab");
+    REQUIRE(name(LSMethod::bicgstab2) == "BICGStab2");
+    REQUIRE(name(LSMethod::gmres) == "GMRES");
+    REQUIRE(name(LSMethod::lu) == "LU");
+    REQUIRE(name(LSMethod::choleski) == "Choleski");
+    REQUIRE_THROWS_WITH(name(LSMethod::LS__end), "unexpected error: Linear system method name is not defined!");
+  }
+
+  SECTION("precond name")
+  {
+    REQUIRE(name(LSPrecond::none) == "none");
+    REQUIRE(name(LSPrecond::diagonal) == "diagonal");
+    REQUIRE(name(LSPrecond::incomplete_choleski) == "ICholeski");
+    REQUIRE(name(LSPrecond::incomplete_LU) == "ILU");
+    REQUIRE(name(LSPrecond::amg) == "AMG");
+    REQUIRE_THROWS_WITH(name(LSPrecond::LS__end),
+                        "unexpected error: Linear system preconditioner name is not defined!");
+  }
+
+  SECTION("library from name")
+  {
+    REQUIRE(LSLibrary::builtin == getLSEnumFromName<LSLibrary>("builtin"));
+    REQUIRE(LSLibrary::petsc == getLSEnumFromName<LSLibrary>("PETSc"));
+
+    REQUIRE_THROWS_WITH(getLSEnumFromName<LSLibrary>("__invalid_lib"),
+                        "error: could not find '__invalid_lib' associate type!");
+  }
+
+  SECTION("method from name")
+  {
+    REQUIRE(LSMethod::cg == getLSEnumFromName<LSMethod>("CG"));
+    REQUIRE(LSMethod::bicgstab == getLSEnumFromName<LSMethod>("BICGStab"));
+    REQUIRE(LSMethod::bicgstab2 == getLSEnumFromName<LSMethod>("BICGStab2"));
+    REQUIRE(LSMethod::lu == getLSEnumFromName<LSMethod>("LU"));
+    REQUIRE(LSMethod::choleski == getLSEnumFromName<LSMethod>("Choleski"));
+    REQUIRE(LSMethod::gmres == getLSEnumFromName<LSMethod>("GMRES"));
+
+    REQUIRE_THROWS_WITH(getLSEnumFromName<LSMethod>("__invalid_method"),
+                        "error: could not find '__invalid_method' associate type!");
+  }
+
+  SECTION("precond from name")
+  {
+    REQUIRE(LSPrecond::none == getLSEnumFromName<LSPrecond>("none"));
+    REQUIRE(LSPrecond::diagonal == getLSEnumFromName<LSPrecond>("diagonal"));
+    REQUIRE(LSPrecond::incomplete_choleski == getLSEnumFromName<LSPrecond>("ICholeski"));
+    REQUIRE(LSPrecond::incomplete_LU == getLSEnumFromName<LSPrecond>("ILU"));
+
+    REQUIRE_THROWS_WITH(getLSEnumFromName<LSPrecond>("__invalid_precond"),
+                        "error: could not find '__invalid_precond' associate type!");
+  }
+
+  SECTION("library list")
+  {
+    std::stringstream os;
+    os << '\n';
+    printLSEnumListNames<LSLibrary>(os);
+
+    const std::string library_list = R"(
+  - builtin
+  - PETSc
+)";
+
+    REQUIRE(os.str() == library_list);
+  }
+
+  SECTION("method list")
+  {
+    std::stringstream os;
+    os << '\n';
+    printLSEnumListNames<LSMethod>(os);
+
+    const std::string library_list = R"(
+  - CG
+  - BICGStab
+  - BICGStab2
+  - GMRES
+  - LU
+  - Choleski
+)";
+
+    REQUIRE(os.str() == library_list);
+  }
+
+  SECTION("precond list")
+  {
+    std::stringstream os;
+    os << '\n';
+    printLSEnumListNames<LSPrecond>(os);
+
+    const std::string library_list = R"(
+  - none
+  - diagonal
+  - ICholeski
+  - ILU
+  - AMG
+)";
+
+    REQUIRE(os.str() == library_list);
+  }
+}
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index 7ab5a66fb05301c3e9e41e4140a898bda96cac76..74f695682992ad951c3b2c33e3c8c0f6b53b657a 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -3,6 +3,7 @@
 
 #include <Kokkos_Core.hpp>
 
+#include <algebra/PETScWrapper.hpp>
 #include <mesh/DiamondDualConnectivityManager.hpp>
 #include <mesh/DiamondDualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -17,6 +18,8 @@ main(int argc, char* argv[])
   parallel::Messenger::create(argc, argv);
   Kokkos::initialize({4, -1, -1, true});
 
+  PETScWrapper::initialize(argc, argv);
+
   Catch::Session session;
   int result = session.applyCommandLine(argc, argv);
 
@@ -46,6 +49,8 @@ main(int argc, char* argv[])
     }
   }
 
+  PETScWrapper::finalize();
+
   Kokkos::finalize();
   parallel::Messenger::destroy();