Skip to content

Commit

Permalink
[core] added matrix library
Browse files Browse the repository at this point in the history
  • Loading branch information
harrand committed Oct 20, 2024
1 parent afec52b commit 5151f12
Show file tree
Hide file tree
Showing 5 changed files with 386 additions and 0 deletions.
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand Down
79 changes: 79 additions & 0 deletions include/tz/core/matrix.hpp
Original file line number Diff line number Diff line change
@@ -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
196 changes: 196 additions & 0 deletions src/tz/core/matrix.cpp
Original file line number Diff line number Diff line change
@@ -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>;
}
6 changes: 6 additions & 0 deletions test/tz/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 5151f12

Please sign in to comment.