diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e112df5a6cff9d7a7a83696017e16dd79ddfb136..06967be54d53c7ff8e88d6ec38d675381cd8ad07 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -86,6 +86,7 @@ add_executable (unit_tests
   test_NameProcessor.cpp
   test_OStreamProcessor.cpp
   test_ParseError.cpp
+  test_PETScUtils.cpp
   test_PugsAssert.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
diff --git a/tests/test_PETScUtils.cpp b/tests/test_PETScUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e0daa8519241fde2f7fd1749b40398977b839de
--- /dev/null
+++ b/tests/test_PETScUtils.cpp
@@ -0,0 +1,117 @@
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <algebra/PETScUtils.hpp>
+
+#include <algebra/CRSMatrixDescriptor.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PETScUtils", "[algebra]")
+{
+  SECTION("PETScAijMatrixEmbedder")
+  {
+    SECTION("from TinyMatrix")
+    {
+      TinyMatrix<3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value = 0;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from DenseMatrix")
+    {
+      DenseMatrix<double> A(3, 3);
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      REQUIRE(petsc_A.numberOfRows() == 3);
+      REQUIRE(petsc_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from DenseMatrix [non-square]")
+    {
+      DenseMatrix<double> A(4, 3);
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      REQUIRE(petsc_A.numberOfRows() == 4);
+      REQUIRE(petsc_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value = 0;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+  }
+
+  SECTION("from CRSMatrix")
+  {
+    Array<int> non_zeros(4);
+    non_zeros[0] = 1;
+    non_zeros[1] = 2;
+    non_zeros[2] = 3;
+    non_zeros[3] = 2;
+
+    CRSMatrixDescriptor<double, int> A(4, 3, non_zeros);
+
+    A(0, 0) = 1;
+    A(1, 0) = 2;
+    A(1, 2) = 3;
+    A(2, 0) = 4;
+    A(2, 1) = 5;
+    A(2, 2) = 6;
+    A(3, 1) = 7;
+    A(3, 2) = 8;
+
+    PETScAijMatrixEmbedder petsc_A{A.getCRSMatrix()};
+
+    auto get_value = [&](int i, int j) {
+      double value;
+      MatGetValue(petsc_A, i, j, &value);
+      return value;
+    };
+
+    REQUIRE(get_value(0, 0) == 1);
+    REQUIRE(get_value(1, 0) == 2);
+    REQUIRE(get_value(1, 2) == 3);
+    REQUIRE(get_value(2, 0) == 4);
+    REQUIRE(get_value(2, 1) == 5);
+    REQUIRE(get_value(2, 2) == 6);
+    REQUIRE(get_value(3, 1) == 7);
+    REQUIRE(get_value(3, 2) == 8);
+  }
+}
+
+#endif   // PUGS_HAS_PETSC