diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2747fa73af..c7e72fd664 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,6 +37,7 @@ topaz_add_library(
 		src/tz/core/error.cpp
 		src/tz/core/job.cpp
 		src/tz/core/vector.cpp
+		src/tz/core/matrix.cpp
 		src/tz/gpu/rhi_vulkan.cpp
 		src/tz/os/impl_win32.cpp
 )
diff --git a/include/tz/core/matrix.hpp b/include/tz/core/matrix.hpp
new file mode 100644
index 0000000000..ede91ff334
--- /dev/null
+++ b/include/tz/core/matrix.hpp
@@ -0,0 +1,79 @@
+#ifndef TOPAZ_CORE_MATRIX_HPP
+#define TOPAZ_CORE_MATRIX_HPP
+#include <array>
+#include <concepts>
+
+namespace tz
+{
+	template<typename T, int N>
+	requires std::integral<T> || std::floating_point<T>
+	struct matrix
+	{
+		static constexpr matrix<T, N> zero()
+		{
+			return matrix<T, N>::filled(T{0});
+		}
+
+		static constexpr matrix<T, N> iden()
+		{
+			auto ret = matrix<T, N>::zero();
+			for(std::size_t i = 0; i < N; i++)
+			{
+				ret.mat[i * N] = T{1};
+			}
+			return ret;
+		}
+
+		static constexpr matrix<T, N> filled(T t)
+		{
+			matrix<T, N> ret;
+			std::fill(ret.mat.begin(), ret.mat.end(), t);
+			return ret;
+		}
+
+		matrix<T, N>() = default;
+		matrix(const matrix<T, N>& cpy) = default;
+		matrix(matrix<T, N>&& move) = default;
+		matrix& operator=(const matrix<T, N>& rhs) = default;
+		matrix& operator=(matrix<T, N>&& rhs) = default;
+
+		const T& operator[](std::size_t idx) const;
+		T& operator[](std::size_t idx);
+
+		const T& operator()(std::size_t x, std::size_t y) const;
+		T& operator()(std::size_t x, std::size_t y);
+
+		matrix<T, N> inverse() const;
+		matrix<T, N> transpose() const;
+
+		// matrix-scalar
+		matrix<T, N>& operator+=(T scalar);
+		matrix<T, N>& operator-=(T scalar);
+		matrix<T, N>& operator*=(T scalar);
+		matrix<T, N>& operator/=(T scalar);
+
+		// matrix-matrix
+		matrix<T, N>& operator*=(const matrix<T, N>& rhs);
+		matrix<T, N> operator*(const matrix<T, N>& rhs){auto cpy = *this; return cpy *= rhs;}
+
+		bool operator==(const matrix<T, N>& rhs) const = default;
+
+		private:
+		std::array<T, N*N> mat;
+	};
+
+	using m2i = matrix<int, 2>;
+	using m3i = matrix<int, 3>;
+	using m4i = matrix<int, 4>;
+	using m2u = matrix<unsigned int, 2>;
+	using m3u = matrix<unsigned int, 3>;
+	using m4u = matrix<unsigned int, 4>;
+	using m2f = matrix<float, 2>;
+	using m3f = matrix<float, 3>;
+	using m4f = matrix<float, 4>;
+	using m2d = matrix<double, 2>;
+	using m3d = matrix<double, 3>;
+	using m4d = matrix<double, 4>;
+}
+
+#endif // TOPAZ_CORE_MATRIX_HPP
\ No newline at end of file
diff --git a/src/tz/core/matrix.cpp b/src/tz/core/matrix.cpp
new file mode 100644
index 0000000000..f2f93ac00f
--- /dev/null
+++ b/src/tz/core/matrix.cpp
@@ -0,0 +1,196 @@
+#include "tz/core/matrix.hpp"
+#include <limits>
+#include <utility>
+#include <cmath>
+
+namespace tz
+{
+#define MATRIX_TEMPLATE_IMPL \
+    template<typename T, int N> \
+    requires std::integral<T> || std::floating_point<T>
+
+	MATRIX_TEMPLATE_IMPL
+	const T& matrix<T, N>::operator[](std::size_t idx) const
+	{
+		return this->mat[idx];
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	T& matrix<T, N>::operator[](std::size_t idx)
+	{
+		return this->mat[idx];
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	const T& matrix<T, N>::operator()(std::size_t x, std::size_t y) const
+	{
+		return this->mat[x + (y * N)];
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	T& matrix<T, N>::operator()(std::size_t x, std::size_t y)
+	{
+		return this->mat[x + (y * N)];
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N> matrix<T, N>::inverse() const
+	{
+		matrix<T, N> inv = matrix<T, N>::iden();  // Start with identity matrix as the inverse
+		matrix<T, N> copy = *this; // Copy of the current matrix
+
+		for(std::size_t i = 0; i < N; i++)
+		{
+			// Find pivot row and swap
+			std::size_t pivot = i;
+			for(std::size_t j = i + 1; j < N; j++)
+			{
+				// abs madness.
+				// this is just: if(abs(mat[j * N + i] > abs(mat[pivot * N + i])))
+				bool cond;
+				if constexpr(std::is_floating_point_v<T>)
+				{
+					cond = std::fabs(copy.mat[j * N + i]) > std::fabs(copy.mat[pivot * N + i]);
+				}
+				else if constexpr(std::is_signed_v<T>)
+				{
+					cond = std::labs(copy.mat[j * N + i]) > std::labs(copy.mat[pivot * N + i]);
+				}
+				else
+				{
+					cond = std::abs(static_cast<std::int64_t>(copy.mat[j * N + i])) > std::abs(static_cast<std::int64_t>(copy.mat[pivot * N + i]));
+				}
+				if(cond)
+				{
+					pivot = j;
+				}
+			}
+
+			if(copy.mat[pivot * N + i] == T{0})
+			{
+				// matrix is singular so there is no inverse.
+				return matrix<T, N>::filled(T{0});
+			}
+
+			// Swap rows in both the copy and the inverse matrix
+			if(pivot != i)
+			{
+				for(std::size_t k = 0; k < N; k++)
+				{
+					std::swap(copy.mat[i * N + k], copy.mat[pivot * N + k]);
+					std::swap(inv.mat[i * N + k], inv.mat[pivot * N + k]);
+				}
+			}
+
+			// Normalize the pivot row
+			T pivot_val = copy.mat[i * N + i];
+			for(std::size_t k = 0; k < N; k++)
+			{
+				copy.mat[i * N + k] /= pivot_val;
+				inv.mat[i * N + k] /= pivot_val;
+			}
+
+			// Eliminate other rows
+			for(std::size_t j = 0; j < N; j++)
+			{
+				if(j != i)
+				{
+					T factor = copy.mat[j * N + i];
+					for(std::size_t k = 0; k < N; k++)
+					{
+						copy.mat[j * N + k] -= factor * copy.mat[i * N + k];
+						inv.mat[j * N + k] -= factor * inv.mat[i * N + k];
+					}
+				}
+			}
+		}
+		return inv;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N> matrix<T, N>::transpose() const
+	{
+		auto cpy = *this;
+		for(std::size_t i = 0; i < N; i++)
+		{
+			for(std::size_t j = 0; j < N; j++)
+			{
+				cpy.mat[j * N + i] = this->mat[i * N + j];
+			}
+		}
+		return cpy;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N>& matrix<T, N>::operator+=(T scalar)
+	{
+		for(std::size_t i = 0; i < N*N; i++)
+		{
+			this->mat[i] += scalar;
+		}
+		return *this;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N>& matrix<T, N>::operator-=(T scalar)
+	{
+		for(std::size_t i = 0; i < N*N; i++)
+		{
+			this->mat[i] -= scalar;
+		}
+		return *this;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N>& matrix<T, N>::operator*=(T scalar)
+	{
+		for(std::size_t i = 0; i < N*N; i++)
+		{
+			this->mat[i] *= scalar;
+		}
+		return *this;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N>& matrix<T, N>::operator/=(T scalar)
+	{
+		for(std::size_t i = 0; i < N*N; i++)
+		{
+			this->mat[i] /= scalar;
+		}
+		return *this;
+	}
+
+	MATRIX_TEMPLATE_IMPL
+	matrix<T, N>& matrix<T, N>::operator*=(const matrix<T, N>& rhs)
+	{
+		matrix<T, N> result = matrix<T, N>::zero(); // Start with zero matrix for the result
+
+		for(std::size_t i = 0; i < N; ++i) 
+		{
+			for(std::size_t j = 0; j < N; ++j) 
+			{
+				for(std::size_t k = 0; k < N; ++k) 
+				{
+					result.mat[i * N + j] += this->mat[i * N + k] * rhs.mat[k * N + j];
+				}
+			}
+		}
+
+		*this = result;
+		return *this;
+	}
+
+    template struct matrix<int, 2>;
+	template struct matrix<int, 3>;
+	template struct matrix<int, 4>;
+    template struct matrix<unsigned int, 2>;
+	template struct matrix<unsigned int, 3>;
+	template struct matrix<unsigned int, 4>;
+    template struct matrix<float, 2>;
+	template struct matrix<float, 3>;
+	template struct matrix<float, 4>;
+    template struct matrix<double, 2>;
+	template struct matrix<double, 3>;
+	template struct matrix<double, 4>;
+}
\ No newline at end of file
diff --git a/test/tz/CMakeLists.txt b/test/tz/CMakeLists.txt
index 84fba60431..1726ee6169 100644
--- a/test/tz/CMakeLists.txt
+++ b/test/tz/CMakeLists.txt
@@ -21,6 +21,12 @@ topaz_add_test(
     vector_test.cpp
 )
 
+topaz_add_test(
+  TARGET tz_matrix_test
+  SOURCES
+    matrix_test.cpp
+)
+
 topaz_add_test(
   TARGET tz_gpu_initialise_test
   SOURCES
diff --git a/test/tz/matrix_test.cpp b/test/tz/matrix_test.cpp
new file mode 100644
index 0000000000..a50e0c576f
--- /dev/null
+++ b/test/tz/matrix_test.cpp
@@ -0,0 +1,104 @@
+#include "tz/core/matrix.hpp"
+#include "tz/topaz.hpp"
+
+template<typename T, int N>
+void test_matrix_constructor()
+{
+    using matrix_t = tz::matrix<T, N>;
+
+    // Test zero matrix creation
+    auto zero_matrix = matrix_t::zero();
+    matrix_t expected_zero = matrix_t::filled(T{0});
+    tz_assert(zero_matrix == expected_zero, "zero() failed. Expected all zero matrix");
+
+    // Test identity matrix creation (only works for square matrices)
+    if constexpr (N == 2 || N == 3 || N == 4) // We only test identity for N >= 2
+    {
+        auto iden_matrix = matrix_t::iden();
+        matrix_t expected_identity = matrix_t::zero();
+        for (std::size_t i = 0; i < N; ++i)
+            expected_identity[i * N] = T{1};  // Diagonal elements = 1
+        tz_assert(iden_matrix == expected_identity, "iden() failed. Expected identity matrix");
+    }
+
+    // Test filled matrix creation
+    auto filled_matrix = matrix_t::filled(static_cast<T>(1.5f));
+    matrix_t expected_filled = matrix_t::filled(static_cast<T>(1.5f));
+    tz_assert(filled_matrix == expected_filled, "filled() failed. Expected all 1.5 matrix");
+}
+
+void test_matrix_accessors()
+{
+    tz::matrix<int, 3> mat;
+    mat[0] = 1;
+    mat[1] = 2;
+    mat[2] = 3;
+    mat[3] = 4;
+    mat[4] = 5;
+    mat[5] = 6;
+    mat[6] = 7;
+    mat[7] = 8;
+    mat[8] = 9;
+
+    // Test accessors (1D indexing)
+    tz_assert(mat[0] == 1, "matrix accessor failed. Expected {}, got {}", 1, mat[0]);
+    tz_assert(mat[4] == 5, "matrix accessor failed. Expected {}, got {}", 5, mat[4]);
+
+    // Test accessors (2D indexing)
+    tz_assert(mat(0, 0) == 1, "matrix accessor (2D) failed. Expected {}, got {}", 1, mat(0, 0));
+    tz_assert(mat(2, 2) == 9, "matrix accessor (2D) failed. Expected {}, got {}", 9, mat(2, 2));
+
+    mat(0, 0) = 10;
+    tz_assert(mat(0, 0) == 10, "matrix accessor (2D) failed after modification. Expected {}, got {}", 10, mat(0, 0));
+}
+
+void test_matrix_operations()
+{
+    // Test matrix addition and subtraction
+    tz::matrix<float, 3> mat1 = tz::matrix<float, 3>::filled(1.0f);
+    tz::matrix<float, 3> mat2 = tz::matrix<float, 3>::filled(2.0f);
+
+    // Test scalar multiplication
+    mat1 *= 2.0f;
+    tz_assert((mat1 == tz::matrix<float, 3>::filled(2.0f)), "Scalar multiplication failed. Expected all 2.0 matrix");
+
+    // Test matrix multiplication
+    tz::matrix<float, 3> mat3 = tz::matrix<float, 3>::iden();
+    mat3 *= mat2;
+    tz_assert((mat3 == tz::matrix<float, 3>::filled(2.0f)), "Matrix multiplication failed. Expected all 2.0 matrix");
+}
+
+void test_matrix_transpose()
+{
+    tz::matrix<int, 2> mat;
+    mat[0] = 1; mat[1] = 2;
+    mat[2] = 3; mat[3] = 4;
+
+    auto transposed = mat.transpose();
+    tz::matrix<int, 2> expected;
+    expected[0] = 1; expected[1] = 3;
+    expected[2] = 2; expected[3] = 4;
+
+    tz_assert(transposed == expected, "transpose() failed. Expected transposed matrix");
+}
+
+int main()
+{
+    test_matrix_constructor<int, 2>();
+    test_matrix_constructor<int, 3>();
+    test_matrix_constructor<int, 4>();
+
+    test_matrix_constructor<float, 2>();
+    test_matrix_constructor<float, 3>();
+    test_matrix_constructor<float, 4>();
+
+    test_matrix_constructor<double, 2>();
+    test_matrix_constructor<double, 3>();
+    test_matrix_constructor<double, 4>();
+
+    test_matrix_accessors();
+    test_matrix_operations();
+    test_matrix_transpose();
+    
+    return 0;
+}
\ No newline at end of file