diff --git a/libs/gfx/include/psemek/gfx/gltf_accessor_iterator.hpp b/libs/gfx/include/psemek/gfx/gltf_accessor_iterator.hpp
index 16af4e84..890b69eb 100644
--- a/libs/gfx/include/psemek/gfx/gltf_accessor_iterator.hpp
+++ b/libs/gfx/include/psemek/gfx/gltf_accessor_iterator.hpp
@@ -93,14 +93,12 @@ namespace psemek::gfx
 
 			static auto pointer(math::matrix<T, R, C> & value)
 			{
-				return &value[0][0];
+				return value.values().begin();
 			}
 
 			static void finalize(math::matrix<T, R, C> & value)
 			{
-				math::matrix<T, C, R> temp;
-				std::copy(value.coords, value.coords + R * C, temp.coords);
-				value = math::transpose(temp);
+				value = math::transpose(value);
 			}
 		};
 
diff --git a/libs/gfx/source/program.cpp b/libs/gfx/source/program.cpp
index 05668990..a092048a 100644
--- a/libs/gfx/source/program.cpp
+++ b/libs/gfx/source/program.cpp
@@ -153,47 +153,47 @@ namespace psemek::gfx
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 2> const & m)
 	{
-		gl::UniformMatrix2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 3> const & m)
 	{
-		gl::UniformMatrix3x2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3x2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 4> const & m)
 	{
-		gl::UniformMatrix4x2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4x2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 2> const & m)
 	{
-		gl::UniformMatrix2x3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2x3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 3> const & m)
 	{
-		gl::UniformMatrix3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 4> const & m)
 	{
-		gl::UniformMatrix4x3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4x3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 2> const & m)
 	{
-		gl::UniformMatrix2x4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2x4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 3> const & m)
 	{
-		gl::UniformMatrix3x4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3x4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 4> const & m)
 	{
-		gl::UniformMatrix4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::interval<int> const & i)
diff --git a/libs/math/include/psemek/math/detail/array.hpp b/libs/math/include/psemek/math/detail/array.hpp
index 91456e95..a86f9f13 100644
--- a/libs/math/include/psemek/math/detail/array.hpp
+++ b/libs/math/include/psemek/math/detail/array.hpp
@@ -11,7 +11,7 @@ namespace psemek::math::detail
 		: util::exception
 	{
 		empty_array_exception(util::stacktrace stacktrace = {})
-			: util::exception("Indexing an empty array", std::move(stacktrace))
+			: util::exception("Indexing into a zero-dimensional array", std::move(stacktrace))
 		{}
 	};
 
@@ -29,7 +29,7 @@ namespace psemek::math::detail
 			T const & operator[](std::size_t) const { throw empty_array_exception{}; }
 			T & operator[](std::size_t) { throw empty_array_exception{}; }
 
-			type operator + (std::size_t) const { return *this; }
+			// type operator + (std::size_t) const { return *this; }
 		};
 	};
 
diff --git a/libs/math/include/psemek/math/detail/array_2d.hpp b/libs/math/include/psemek/math/detail/array_2d.hpp
new file mode 100644
index 00000000..bd257e33
--- /dev/null
+++ b/libs/math/include/psemek/math/detail/array_2d.hpp
@@ -0,0 +1,218 @@
+#pragma once
+
+#include <psemek/math/detail/array.hpp>
+
+namespace psemek::math::detail
+{
+
+	template <typename T, std::size_t R, std::size_t C>
+	struct array_2d
+	{
+		using type = T[R][C];
+	};
+
+	template <typename T, std::size_t R>
+	struct array_2d<T, R, 0>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = R;
+			static constexpr std::size_t columns = 0;
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, dynamic, 0>
+	{
+		struct type
+		{
+			std::size_t rows;
+			static constexpr std::size_t columns = 0;
+
+			type(std::size_t rows)
+				: rows(rows)
+			{}
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T, std::size_t C>
+	struct array_2d<T, 0, C>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			static constexpr std::size_t columns = C;
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, 0, dynamic>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			std::size_t columns;
+
+			type(std::size_t columns)
+				: columns(columns)
+			{}
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, 0, 0>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			static constexpr std::size_t columns = 0;
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T, std::size_t R>
+	struct array_2d<T, R, dynamic>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = R;
+			std::size_t columns;
+			std::unique_ptr<T[]> data;
+
+			type()
+				: columns(0)
+			{}
+
+			type(std::size_t columns)
+				: columns(columns)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			type([[maybe_unused]] std::size_t rows, std::size_t columns)
+				: type(columns)
+			{
+				assert(rows == R);
+			}
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.columns = columns;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+	template <typename T, std::size_t C>
+	struct array_2d<T, dynamic, C>
+	{
+		struct type
+		{
+			std::size_t rows;
+			static constexpr std::size_t columns = C;
+			std::unique_ptr<T[]> data;
+
+			type()
+				: rows(0)
+			{}
+
+			type(std::size_t rows)
+				: rows(rows)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			type(std::size_t rows, [[maybe_unused]] std::size_t columns)
+				: type(rows)
+			{
+				assert(columns == C);
+			}
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.rows = rows;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, dynamic, dynamic>
+	{
+		struct type
+		{
+			std::size_t rows;
+			std::size_t columns;
+			std::unique_ptr<T[]> data;
+
+			type()
+				: rows(0)
+				, columns(0)
+			{}
+
+			type(std::size_t rows, std::size_t columns)
+				: rows(rows)
+				, columns(columns)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.rows = rows;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+}
diff --git a/libs/math/include/psemek/math/matrix.hpp b/libs/math/include/psemek/math/matrix.hpp
index d3f81f6c..ecc5cfa3 100644
--- a/libs/math/include/psemek/math/matrix.hpp
+++ b/libs/math/include/psemek/math/matrix.hpp
@@ -1,8 +1,9 @@
 #pragma once
 
-#include <psemek/math/detail/array.hpp>
+#include <psemek/math/detail/array_2d.hpp>
 #include <psemek/math/vector.hpp>
 #include <psemek/math/math.hpp>
+#include <psemek/util/range.hpp>
 
 #include <iostream>
 #include <iomanip>
@@ -18,26 +19,64 @@ namespace psemek::math
 		static constexpr std::size_t static_rows = R;
 		static constexpr std::size_t static_columns = C;
 
-		typename detail::array<T, R * C>::type coords;
+		typename detail::array_2d<T, R, C>::type coords;
 
 		std::size_t rows() const
 		{
-			return R;
+			if constexpr (R == dynamic)
+			{
+				return coords.rows;
+			}
+			else
+			{
+				return R;
+			}
 		}
 
 		std::size_t columns() const
 		{
-			return C;
+			if constexpr (C == dynamic)
+			{
+				return coords.columns;
+			}
+			else
+			{
+				return C;
+			}
 		}
 
 		auto operator[](std::size_t i)
 		{
-			return coords + C * i;
+			return coords[i];
 		}
 
 		auto operator[](std::size_t i) const
 		{
-			return coords + C * i;
+			return coords[i];
+		}
+
+		util::range<T *> values()
+		{
+			return {&coords[0][0], &coords[0][0] + rows() * columns()};
+		}
+
+		util::range<T const *> values() const
+		{
+			return {&coords[0][0], &coords[0][0] + rows() * columns()};
+		}
+
+		matrix copy() const
+		{
+			if constexpr (R == dynamic || C == dynamic)
+			{
+				matrix r;
+				r.coords = coords.copy();
+				return r;
+			}
+			else
+			{
+				return *this;
+			}
 		}
 
 		matrix & operator *= (T const & s);
@@ -56,8 +95,8 @@ namespace psemek::math
 	matrix<T, R, C> matrix<T, R, C>::zero()
 	{
 		matrix<T, R, C> m;
-		for (std::size_t i = 0; i < R * C; ++i)
-			m.coords[i] = T(0);
+		for (auto & v : m.values())
+			v = T(0);
 		return m;
 	}
 
@@ -89,20 +128,24 @@ namespace psemek::math
 	matrix<T1, R, C> cast(matrix<T, R, C> const & m)
 	{
 		matrix<T1, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = static_cast<T1>(m.coords[i]);
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = static_cast<T1>(m[i][j]);
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	std::strong_ordering operator <=> (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		for (std::size_t i = 0; i < R * C; ++i)
+		for (std::size_t i = 0; i < R; ++i)
 		{
-			if (m1.coords[i] < m2.coords[i])
-				return std::strong_ordering::less;
-			else if (m1.coords[i] > m2.coords[i])
-				return std::strong_ordering::greater;
+			for (std::size_t j = 0; j < C; ++j)
+			{
+				if (m1[i][j] < m2[i][j])
+					return std::strong_ordering::less;
+				else if (m1[i][j] > m2[i][j])
+					return std::strong_ordering::greater;
+			}
 		}
 		return std::strong_ordering::equal;
 	}
@@ -147,8 +190,9 @@ namespace psemek::math
 	matrix<T, R, C> operator * (matrix<T, R, C> const & m, T const & s)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m.coords[i] * s;
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = m[i][j] * s;
 		return r;
 	}
 
@@ -156,8 +200,9 @@ namespace psemek::math
 	matrix<T, R, C> operator * (T const & s, matrix<T, R, C> const & m)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = s * m.coords[i];
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = m[i][j] * s;
 		return r;
 	}
 
@@ -165,8 +210,9 @@ namespace psemek::math
 	matrix<T, R, C> operator / (matrix<T, R, C> const & m, T const & s)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m.coords[i] / s;
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = m[i][j] / s;
 		return r;
 	}
 
@@ -188,8 +234,9 @@ namespace psemek::math
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = -m.coords[i];
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = -m[i][j];
 		return r;
 	}
 
@@ -197,8 +244,9 @@ namespace psemek::math
 	matrix<T, R, C> operator + (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m1.coords[i] + m2.coords[i];
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = m1[i][j] + m2[i][j];
 		return r;
 	}
 
@@ -206,8 +254,9 @@ namespace psemek::math
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
 		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m1.coords[i] - m2.coords[i];
+		for (std::size_t i = 0; i < R; ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				r[i][j] = m1[i][j] - m2[i][j];
 		return r;
 	}
 
@@ -283,11 +332,11 @@ namespace psemek::math
 	{
 		vector<T, N> m[] = {v, rows...};
 
-		matrix<T, sizeof...(Rows) + 1, N> result;
-		for (std::size_t i = 0; i < result.rows(); ++i)
-			for (std::size_t j = 0; j < result.columns(); ++j)
-				result[i][j] = m[i][j];
-		return result;
+		matrix<T, sizeof...(Rows) + 1, N> r;
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
+				r[i][j] = m[i][j];
+		return r;
 	}
 
 	template <typename T, std::size_t N, typename ... Columns>
@@ -295,11 +344,11 @@ namespace psemek::math
 	{
 		vector<T, N> m[] = {v, columns...};
 
-		matrix<T, N, sizeof...(Columns) + 1> result;
-		for (std::size_t i = 0; i < result.rows(); ++i)
-			for (std::size_t j = 0; j < result.columns(); ++j)
-				result[i][j] = m[j][i];
-		return result;
+		matrix<T, N, sizeof...(Columns) + 1> r;
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
+				r[i][j] = m[j][i];
+		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
diff --git a/libs/math/include/psemek/math/vector.hpp b/libs/math/include/psemek/math/vector.hpp
index e53a5ffb..0be0e1e7 100644
--- a/libs/math/include/psemek/math/vector.hpp
+++ b/libs/math/include/psemek/math/vector.hpp
@@ -51,6 +51,20 @@ namespace psemek::math
 			}
 		}
 
+		vector copy() const
+		{
+			if constexpr (N == dynamic)
+			{
+				vector result;
+				result.coords = coords.copy();
+				return result;
+			}
+			else
+			{
+				return *this;
+			}
+		}
+
 		T & operator[](std::size_t i)
 		{
 			return coords[i];