diff --git a/libs/geom/include/psemek/geom/qr.hpp b/libs/geom/include/psemek/geom/qr.hpp
index 3ec82838..2c892793 100644
--- a/libs/geom/include/psemek/geom/qr.hpp
+++ b/libs/geom/include/psemek/geom/qr.hpp
@@ -9,9 +9,10 @@ namespace psemek::geom
 	struct qr_result
 	{
 		matrix<T, N, N> q; // orthogonal
-		matrix<T, N, N> r; // upper-triangular
+		matrix<T, N, N> r; // upper-triangular (for decomposition) or block upper-triangular (for eigenvalues)
 	};
 
+	// returns {q, r} such that m = q * r
 	template <typename T, std::size_t N>
 	qr_result<T, N> qr_decomposition(matrix<T, N, N> const & m)
 	{
@@ -59,4 +60,22 @@ namespace psemek::geom
 		return {q, r};
 	}
 
+	// returns {q, r} such that m = q * r * q^(-1)
+	// NB: q^(-1) = q^T
+	template <typename T, std::size_t N>
+	qr_result<T, N> qr_eigenvalues(matrix<T, N, N> const & m, std::size_t const iterations)
+	{
+		matrix<T, N, N> r = m;
+		matrix<T, N, N> q = matrix<T, N, N>::identity();
+
+		for (std::size_t iteration = 0; iteration < iterations; ++iteration)
+		{
+			auto qr = qr_decomposition(r);
+			r = qr.r * qr.q;
+			q = q * qr.q;
+		}
+
+		return {q, r};
+	}
+
 }
diff --git a/libs/geom/tests/matrix.cpp b/libs/geom/tests/matrix.cpp
index 2debf5f1..16fd5626 100644
--- a/libs/geom/tests/matrix.cpp
+++ b/libs/geom/tests/matrix.cpp
@@ -217,3 +217,70 @@ test_case(geom_matrix_qr)
 				expect_small(result.r[j][i], 1e-6);
 	}
 }
+
+test_case(geom_matrix_qr__eig_symmetric)
+{
+	random::generator rng;
+	random::normal_distribution<float> normal;
+
+	for (int iteration = 0; iteration < 64; ++iteration)
+	{
+		matrix<float, 8, 8> m;
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = 0; j < m.columns(); ++j)
+				m[i][j] = normal(rng);
+
+		m = transpose(m) * m;
+
+		auto result = qr_eigenvalues(m, 1024);
+
+		expect_small(frobenius_norm(m * result.q - result.q * result.r), 1e-3);
+		expect_small(frobenius_norm(matrix<float, 8, 8>::identity() - result.q * transpose(result.q)), 1e-4);
+
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = i + 1; j < m.columns(); ++j)
+				expect_small(result.r[j][i], 1e-3);
+
+		std::string test;
+
+		for (std::size_t i = 0; i < m.columns(); ++i)
+		{
+			auto v = column(result.q, i);
+			expect_small(length(m * v - result.r[i][i] * v), 1e-3);
+		}
+	}
+}
+
+test_case(geom_matrix_qr__eig_general)
+{
+	random::generator rng;
+	random::normal_distribution<float> normal;
+
+	for (int iteration = 0; iteration < 64; ++iteration)
+	{
+		matrix<float, 8, 8> m;
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = 0; j < m.columns(); ++j)
+				m[i][j] = normal(rng);
+
+		auto result = qr_eigenvalues(m, 4096);
+
+		expect_small(frobenius_norm(m * result.q - result.q * result.r), 1e-2);
+		expect_small(frobenius_norm(matrix<float, 8, 8>::identity() - result.q * transpose(result.q)), 1e-3);
+
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = i + 2; j < m.columns(); ++j)
+				expect_small(result.r[j][i], 1e-2);
+
+		for (std::size_t i = 0; i < m.columns(); ++i)
+		{
+			// Don't check complex eigenvalues
+			if (i > 0 && std::abs(result.r[i][i - 1]) > 1e-3f) continue;
+			if (i + 1 < m.rows() && std::abs(result.r[i + 1][i]) > 1e-3f) continue;
+
+			// TODO: properly check eigenvalues
+//			auto v = column(result.q, i);
+//			expect_small(length(m * v - result.r[i][i] * v), 1e-3);
+		}
+	}
+}