diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c8fa8e2c5e2b3d073dcaefdb586cce29035dcec8..2a5f5dfa3358649feda41dc4ea8e7318da0e8d53 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,4 +9,5 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # ---------------------------------------------------------------
 
-add_subdirectory(fortran)
+# add_subdirectory(fortran)
+add_subdirectory(c)
diff --git a/test/c/CMakeLists.txt b/test/c/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..be1af9e339fed50fd9967b9a0f2ba7466d75ee3c
--- /dev/null
+++ b/test/c/CMakeLists.txt
@@ -0,0 +1,30 @@
+# Fetch GoogleTest via FetchContent
+include(FetchContent)
+FetchContent_Declare(
+  googletest
+  URL https://github.com/google/googletest/archive/refs/tags/release-1.12.1.zip
+)
+set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+FetchContent_MakeAvailable(googletest)
+
+# Find Kokkos (or use your existing Kokkos installation)
+# find_package(Kokkos REQUIRED)
+
+set(SOURCES
+  main.cpp
+  test_tdma_solver.cpp
+)
+# Create the test executable from your test files, including main.cpp.
+add_executable(iconmath_test_c ${SOURCES})
+
+# Link the test executable with GoogleTest and Kokkos.
+target_link_libraries(iconmath_test_c
+  PUBLIC
+    iconmath-support
+  PRIVATE
+    gtest_main
+    Kokkos::kokkos
+)
+
+include(GoogleTest)
+gtest_discover_tests(iconmath_test_c)
diff --git a/test/c/main.cpp b/test/c/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2df720d1cecff67239f5cb2efe00b5f3a7f86f9c
--- /dev/null
+++ b/test/c/main.cpp
@@ -0,0 +1,14 @@
+#include <Kokkos_Core.hpp>
+#include <gtest/gtest.h>
+
+int main(int argc, char** argv) {
+  // Initialize Kokkos before any tests run.
+  Kokkos::initialize(argc, argv);
+  
+  ::testing::InitGoogleTest(&argc, argv);
+  int result = RUN_ALL_TESTS();
+  
+  // Finalize Kokkos after all tests have completed.
+  Kokkos::finalize();
+  return result;
+}
diff --git a/test/c/test_tdma_solver.cpp b/test/c/test_tdma_solver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7c3c3a8c6519fd3c40e44769df20fe07d5d10486
--- /dev/null
+++ b/test/c/test_tdma_solver.cpp
@@ -0,0 +1,77 @@
+#include <gtest/gtest.h>
+#include <vector>
+#include <algorithm>
+#include "mo_math_utilities.hpp"
+
+// Helper function to compute the 1D index for column-major storage.
+inline int idx(int i, int j, int nrows) {
+  return i + j * nrows;
+}
+
+// Test fixture for the TDMA solver tests.
+class TDMASolverTestFixture : public ::testing::Test {
+protected:
+  const int n = 10;             // Matrix dimension.
+  std::vector<double> a;        // Input matrix a.
+  std::vector<double> b;        // Input matrix b.
+  std::vector<double> c;        // Input matrix c.
+  std::vector<double> d;        // Input matrix d.
+  std::vector<double> x;        // Output matrix.
+
+  TDMASolverTestFixture()
+      : a(n * n), b(n * n), c(n * n), d(n * n), x(n * n, 0.0) {}
+
+  // SetUp is run before each test.
+  void SetUp() override {
+    // Fill arrays in column-major order.
+    for (int j = 0; j < n; j++) {
+      for (int i = 0; i < n; i++) {
+        double value = (i + 1) + (j + 1);
+        a[idx(i, j, n)] = 1.0 * value;
+        b[idx(i, j, n)] = 2.0 * value;
+        c[idx(i, j, n)] = 1.0 * value;
+        d[idx(i, j, n)] = 1.0 * value;
+      }
+    }
+    // Clear the output vector.
+    std::fill(x.begin(), x.end(), 0.0);
+  }
+};
+
+TEST_F(TDMASolverTestFixture, FullTest) {
+  // Call the solver over the full range:
+  tdma_solver_vec_double(a.data(), b.data(), c.data(), d.data(),
+                         0, n, 0, n, n, n, x.data());
+
+  // Compute the sum of all elements in the output matrix.
+  double sum = 0.0;
+  for (int j = 0; j < n; j++) {
+    for (int i = 0; i < n; i++) {
+      sum += x[idx(i, j, n)];
+    }
+  }
+
+  // Expected reference sum
+  double sum_ref = 27.2727272727272769;
+  double tol = 1e-13;
+  EXPECT_NEAR(sum, sum_ref, tol);
+}
+
+TEST_F(TDMASolverTestFixture, PartialTest) {
+  // Call the solver for a partial region:
+  // For C++: slev = 1, elev = n-1, startidx = 1, endidx = n-1.
+  tdma_solver_vec_double(a.data(), b.data(), c.data(), d.data(),
+                         1, n - 1, 1, n - 1, n, n, x.data());
+
+  // Compute the sum over a region
+  double sum = 0.0;
+  for (int j = 1; j < n - 1; j++) {
+    for (int i = 1; i < n - 1; i++) {
+      sum += x[idx(i, j, n)];
+    }
+  }
+
+  double sum_ref = 17.7777777777777679;
+  double tol = 1e-13;
+  EXPECT_NEAR(sum, sum_ref, tol);
+}