codereckons/rotgen
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Implement Eigen arithmetic functions stub

Browse source
This commit is contained in:
Joel Falcou 2025-05-20 15:54:34 +02:00
commit 03591ac323
6 changed files with 501 additions and 12 deletions
Download patch file Download diff file Expand all files Collapse all files

19
include/rotgen/impl/matrix_impl64.hpp
View file

@ -36,6 +36,22 @@ namespace rotgen
void resize(std::size_t new_rows, std::size_t new_cols);
void conservativeResize(std::size_t new_rows, std::size_t new_cols);
matrix_impl64 transpose() const;
matrix_impl64 conjugate() const;
matrix_impl64 adjoint() const;
void transposeInPlace();
void adjointInPlace();
double sum() const;
double prod() const;
double mean() const;
double maxCoeff() const;
double maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
double minCoeff() const;
double minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
double trace() const;
double& operator()(std::size_t i, std::size_t j);
double const& operator()(std::size_t i, std::size_t j) const;
@ -43,8 +59,11 @@ namespace rotgen
double const& operator()(std::size_t index) const;
matrix_impl64& operator+=(matrix_impl64 const& rhs);
matrix_impl64& operator-=(matrix_impl64 const& rhs);
matrix_impl64 operator-() const;
matrix_impl64& operator*=(matrix_impl64 const& rhs);
matrix_impl64& operator*=(double d);
matrix_impl64& operator/=(double d);
friend std::ostream& operator<<(std::ostream&,matrix_impl64 const&);
friend bool operator==(matrix_impl64 const& lhs, matrix_impl64 const& rhs);

52
include/rotgen/matrix.hpp
View file

@ -29,6 +29,8 @@ namespace rotgen
if constexpr(Cols != -1) assert(c == Cols && "Mismatched between dynamic and static column size");
}
matrix(parent const& base) : parent(base) {}
matrix(std::initializer_list<std::initializer_list<Scalar>> init) : parent(init)
{
if constexpr(Rows != -1) assert(init.size() == Rows && "Mismatched between dynamic and static row size");
@ -58,6 +60,25 @@ namespace rotgen
parent::conservativeResize(new_rows, new_cols);
}
matrix transpose() const
{
return matrix(static_cast<parent const&>(*this).transpose());
}
matrix conjugate() const
{
return matrix(static_cast<parent const&>(*this).conjugate());
}
matrix adjoint() const
{
return matrix(static_cast<parent const&>(*this).adjoint());
}
void transposeInPlace() { parent::transposeInPlace(); }
void adjointInPlace() { parent::adjointInPlace(); }
friend bool operator==(matrix const& lhs, matrix const& rhs)
{
return static_cast<parent const&>(lhs) == static_cast<parent const&>(rhs);
@ -69,6 +90,17 @@ namespace rotgen
return *this;
}
matrix& operator-=(matrix const& rhs)
{
static_cast<parent&>(*this) -= static_cast<parent const&>(rhs);
return *this;
}
matrix operator-() const
{
return matrix(static_cast<parent const&>(*this).operator-());
}
matrix& operator*=(matrix const& rhs)
{
static_cast<parent&>(*this) *= static_cast<parent const&>(rhs);
@ -80,6 +112,12 @@ namespace rotgen
static_cast<parent&>(*this) *= rhs;
return *this;
}
matrix& operator/=(double rhs)
{
static_cast<parent&>(*this) /= rhs;
return *this;
}
};
template<typename S, int R, int C, int O, int MR, int MC>
@ -89,6 +127,13 @@ namespace rotgen
return that += rhs;
}
template<typename S, int R, int C, int O, int MR, int MC>
matrix<S,R,C,O,MR,MC> operator-(matrix<S,R,C,O,MR,MC> const& lhs, matrix<S,R,C,O,MR,MC> const& rhs)
{
matrix<S,R,C,O,MR,MC> that(lhs);
return that -= rhs;
}
template<typename S, int R, int C, int O, int MR, int MC>
matrix<S,R,C,O,MR,MC> operator*(matrix<S,R,C,O,MR,MC> const& lhs, matrix<S,R,C,O,MR,MC> const& rhs)
{
@ -108,4 +153,11 @@ namespace rotgen
{
return rhs * lhs;
}
template<typename S, int R, int C, int O, int MR, int MC>
matrix<S,R,C,O,MR,MC> operator/(matrix<S,R,C,O,MR,MC> const& lhs, double rhs)
{
matrix<S,R,C,O,MR,MC> that(lhs);
return that /= rhs;
}
}

63
src/matrix_impl64.cpp
View file

@ -76,6 +76,50 @@ namespace rotgen
const double* matrix_impl64::data() const { return storage_->data.data(); }
matrix_impl64 matrix_impl64::transpose() const
{
matrix_impl64 result(*this);
result.storage_->data.transposeInPlace();
return result;
}
matrix_impl64 matrix_impl64::conjugate() const
{
matrix_impl64 result(*this);
result.storage_->data = storage_->data.conjugate();
return result;
}
matrix_impl64 matrix_impl64::adjoint() const
{
matrix_impl64 result(*this);
result.storage_->data.adjointInPlace();
return result;
}
void matrix_impl64::transposeInPlace()
{
storage_->data.transposeInPlace();
}
void matrix_impl64::adjointInPlace()
{
storage_->data.adjointInPlace();
}
double matrix_impl64::sum() const { return storage_->data.sum(); }
double matrix_impl64::prod() const { return storage_->data.prod(); }
double matrix_impl64::mean() const { return storage_->data.mean(); }
double matrix_impl64::minCoeff() const { return storage_->data.minCoeff(); }
double matrix_impl64::minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.minCoeff(row, col); }
double matrix_impl64::maxCoeff() const { return storage_->data.maxCoeff(); }
double matrix_impl64::maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.maxCoeff(row, col); }
double matrix_impl64::trace() const { return storage_->data.trace(); }
//==================================================================================================
// Operators
//==================================================================================================
@ -95,6 +139,19 @@ namespace rotgen
return *this;
}
matrix_impl64& matrix_impl64::operator-=(matrix_impl64 const& rhs)
{
storage_->data -= rhs.storage_->data;
return *this;
}
matrix_impl64 matrix_impl64::operator-() const
{
matrix_impl64 result(*this);
result.storage_->data = -result.storage_->data;
return result;
}
matrix_impl64& matrix_impl64::operator*=(matrix_impl64 const& rhs)
{
storage_->data *= rhs.storage_->data;
@ -106,4 +163,10 @@ namespace rotgen
storage_->data *= s;
return *this;
}
matrix_impl64& matrix_impl64::operator/=(double s)
{
storage_->data /= s;
return *this;
}
}

2
test/CMakeLists.txt
View file

@ -12,7 +12,7 @@ rotgen_setup_test_targets()
add_library(rotgen_test INTERFACE)
target_compile_options(rotgen_test INTERFACE -std=c++20 -Werror -Wall -Wextra -Wshadow -Wunused-variable)
target_include_directories( rotgen_test INTERFACE ${PROJECT_SOURCE_DIR}/test)
target_link_libraries(rotgen_test INTERFACE rotgen)
target_link_libraries(rotgen_test INTERFACE rotgen Eigen3::Eigen)
##======================================================================================================================
## Unit Tests registration

121
test/basic/arithmetic_functions.cpp Normal file
View file

@ -0,0 +1,121 @@
//==================================================================================================
/*
ROTGEN - Runtime Overlay for Eigen
Copyright : CODE RECKONS
SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#define TTS_MAIN
#include <rotgen/matrix.hpp>
#include "tts.hpp"
#include <Eigen/Dense>
template <typename MatrixType>
struct MatrixDescriptor
{
std::size_t rows, cols;
std::function<void(MatrixType &, std::size_t, std::size_t)> init_fn;
};
template <typename MatrixType>
void test_matrix_unary_ops(std::size_t rows, std::size_t cols,
const std::function<void(MatrixType &, std::size_t, std::size_t)> &init_fn)
{
MatrixType original(rows, cols);
MatrixType transposed_matrix(cols, rows);
for (std::size_t r = 0; r < rows; ++r)
for (std::size_t c = 0; c < cols; ++c)
init_fn(original, r, c);
for (std::size_t r = 0; r < rows; ++r)
for (std::size_t c = 0; c < cols; ++c)
transposed_matrix(c, r) = original(r, c);
TTS_EQUAL(original.transpose(), transposed_matrix);
TTS_EQUAL(original.conjugate(), original);
TTS_EQUAL(original.adjoint(), transposed_matrix);
original.transposeInPlace();
TTS_EQUAL(original, transposed_matrix);
original.transposeInPlace();
original.adjointInPlace();
TTS_EQUAL(original, transposed_matrix);
}
template <typename MatrixType>
void test_matrix_scalar_reductions(std::size_t rows, std::size_t cols,
const std::function<void(MatrixType&, std::size_t, std::size_t)>& init_fn)
{
using EigenMatrix = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
MatrixType original(rows, cols);
EigenMatrix ref(rows, cols);
for (std::size_t r = 0; r < rows; ++r)
for (std::size_t c = 0; c < cols; ++c)
{
init_fn(original, r, c);
ref(r, c) = original(r, c);
}
TTS_EQUAL(original.sum(), ref.sum());
TTS_EQUAL(original.prod(), ref.prod());
TTS_EQUAL(original.mean(), ref.mean());
TTS_EQUAL(original.maxCoeff(), ref.maxCoeff());
TTS_EQUAL(original.minCoeff(), ref.minCoeff());
TTS_EQUAL(original.trace(), ref.trace());
std::ptrdiff_t row, col, ref_row, ref_col;
double cmin = original.minCoeff(&row, &col);
double rmin = ref.minCoeff(&ref_row, &ref_col);
TTS_EQUAL(cmin, rmin);
TTS_EQUAL(row, ref_row);
TTS_EQUAL(col, ref_col);
double cmax = original.maxCoeff(&row, &col);
double rmax = ref.maxCoeff(&ref_row, &ref_col);
TTS_EQUAL(cmax, rmax);
TTS_EQUAL(row, ref_row);
TTS_EQUAL(col, ref_col);
}
TTS_CASE("Matrix unary operations: transpose, adjoint, conjugate")
{
std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
{3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + 3 * c - 2.5; }},
{4, 9, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r*r + 3.12 * c + 6.87; }},
{2, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 1.1 * (r - c); }},
{1, 5, [](auto &m, std::size_t r, std::size_t c) { m(r, c) = 9.99; }},
{4, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
{3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c) ? 1.0 : 0.0; }},
{2, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r + c) * 1e-10; }},
{2, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r + 1) * 1e+10; }},
};
for (auto const &desc : test_matrices)
test_matrix_unary_ops<rotgen::matrix<double>>(desc.rows, desc.cols, desc.init_fn);
};
TTS_CASE("Basic arithmetic reduction operations")
{
std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
{3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + c; }},
{3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
{2, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -r -c*c - 1234; }},
{4, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 7.0; }},
{1, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 42.0; }},
{4, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = std::sin(r + c); }},
{1, 5, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -1.5 * r + 2.56 * c; }},
{5, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
};
for (const auto& [rows, cols, init_fn] : test_matrices)
test_matrix_scalar_reductions<rotgen::matrix<double>>(rows, cols, init_fn);
};

254
test/basic/operators.cpp
View file

@ -9,23 +9,257 @@
#include <rotgen/matrix.hpp>
#include "tts.hpp"
TTS_CASE("Check operator*")
template <typename MatrixType>
void test_matrix_scalar_multiplication(std::size_t rows, std::size_t cols, double scalar, auto init_fn)
{
rotgen::matrix<double> a(2,2);
rotgen::matrix<double> ref(2,2);
MatrixType a(rows, cols);
MatrixType ref(rows, cols);
for(std::size_t r=0;r<a.rows();r++)
for (std::size_t r = 0; r < rows; ++r)
{
for(std::size_t c=0;c<a.cols();c++)
for (std::size_t c = 0; c < cols; ++c)
{
a(r,c) = (1+c) + 10*(1+r);
ref(r,c) = ((1+c) + 10*(1+r)) * 10.5;
init_fn(a, r, c);
ref(r, c) = a(r, c) * scalar;
}
}
TTS_EQUAL(a * 10.5, ref);
TTS_EQUAL(10.5 * a, ref);
TTS_EQUAL(a * scalar, ref);
TTS_EQUAL(scalar * a, ref);
a *= 10.5;
a *= scalar;
TTS_EQUAL(a, ref);
}
template <typename MatrixType>
void test_matrix_scalar_division(std::size_t rows, std::size_t cols, double scalar, auto init_fn)
{
MatrixType a(rows, cols);
MatrixType ref(rows, cols);
for (std::size_t r = 0; r < rows; ++r)
{
for (std::size_t c = 0; c < cols; ++c)
{
init_fn(a, r, c);
ref(r, c) = a(r, c) / scalar;
}
}
TTS_EQUAL(a / scalar, ref);
a /= scalar;
TTS_EQUAL(a, ref);
}
template <typename MatrixType>
void test_matrix_multiplication(std::size_t n, std::size_t m, std::size_t p, auto a_init_fn, auto b_init_fn)
{
MatrixType a(n, m);
MatrixType b(m, p);
MatrixType ref(n, p);
for (std::size_t r = 0; r < n; ++r)
for (std::size_t c = 0; c < m; ++c)
a_init_fn(a, r, c);
for (std::size_t r = 0; r < m; ++r)
for (std::size_t c = 0; c < p; ++c)
b_init_fn(b, r, c);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < p; ++j)
{
ref(i, j) = 0;
for (std::size_t k = 0; k < m; ++k)
ref(i, j) += a(i, k) * b(k, j);
}
TTS_EQUAL(a * b, ref);
a *= b;
TTS_EQUAL(a, ref);
}
template <typename MatrixType>
void test_matrix_addition(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn)
{
MatrixType a(rows, cols);
MatrixType b(rows, cols);
MatrixType ref(rows, cols);
for (std::size_t r = 0; r < rows; ++r)
{
for (std::size_t c = 0; c < cols; ++c)
{
a_init_fn(a, r, c);
b_init_fn(b, r, c);
ref(r, c) = a(r,c) + b(r,c);
}
}
TTS_EQUAL(a + b, ref);
TTS_EQUAL(b + a, ref);
a += b;
TTS_EQUAL(a, ref);
}
template <typename MatrixType>
void test_matrix_substraction(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn)
{
MatrixType a(rows, cols);
MatrixType b(rows, cols);
MatrixType ref(rows, cols);
MatrixType a_minus_ref(rows, cols);
MatrixType b_minus_ref(rows, cols);
for (std::size_t r = 0; r < rows; ++r)
{
for (std::size_t c = 0; c < cols; ++c)
{
a_init_fn(a, r, c);
b_init_fn(b, r, c);
ref(r, c) = a(r,c) - b(r,c);
a_minus_ref(r, c) = -a(r,c);
b_minus_ref(r, c) = -b(r,c);
}
}
MatrixType a_unary = -a;
MatrixType b_unary = -b;
TTS_EQUAL(a - b, ref);
TTS_EQUAL(a_unary, a_minus_ref);
TTS_EQUAL(-a, a_minus_ref);
TTS_EQUAL(-b, b_minus_ref);
TTS_EQUAL(-(-a), a);
TTS_EQUAL(-(-b), b);
a -= b;
TTS_EQUAL(a, ref);
};
TTS_CASE("Check matrix * scalar and scalar * matrix with default values")
{
test_matrix_scalar_multiplication<rotgen::matrix<double>>(2, 2, 10.5,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = (1 + c) + 10 * (1 + r); });
};
TTS_CASE("Check matrix * scalar with zero scalar multiplication")
{
test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, 0.0,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 5 * c - r; });
};
TTS_CASE("Check matrix * scalar with one scalar multiplication")
{
test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, 1,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 3.3*r - 6; });
};
TTS_CASE("Check matrix - scalar with negative scalar multiplication")
{
test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, -36.2,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = r * r - c + 3.9; });
};
TTS_CASE("Check static matrix - scalar with float scalar multiplication")
{
test_matrix_scalar_multiplication<rotgen::matrix<double, 3, 4>>(3, 4, 5.6,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 1.2*r+4.5*c-6; });
};
TTS_CASE("Check matrix / scalar with default values")
{
test_matrix_scalar_division<rotgen::matrix<double>>(6, 7, 10.5,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = (1 + c) + 10 * (1 + r); });
};
TTS_CASE("Check matrix * scalar with one scalar multiplication")
{
test_matrix_scalar_division<rotgen::matrix<double>>(3, 2, 1,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 3.3*r - 4.5*c + 1.1; });
};
TTS_CASE("Check matrix - scalar with negative scalar multiplication")
{
test_matrix_scalar_division<rotgen::matrix<double>>(2, 7, -36.2,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 3.4 * r * r - c + 3.9; });
};
TTS_CASE("Check static matrix - scalar with float scalar multiplication")
{
test_matrix_scalar_division<rotgen::matrix<double, 3, 4>>(3, 4, 5.6,
[](auto& a, std::size_t r, std::size_t c)
{ a(r, c) = 1.2*r+4.5*c-6; });
};
TTS_CASE("Matrix multiplication")
{
test_matrix_multiplication<rotgen::matrix<int>>(2, 3, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r + c; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r * c; });
};
TTS_CASE("Matrix multiplication with floats")
{
test_matrix_multiplication<rotgen::matrix<double>>(3, 6, 3,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.4 * r + 5.4 * c*c - 5.2; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r * c + 5 * r - 3.33; });
};
TTS_CASE("Matrix multiplication with zero matrix")
{
test_matrix_multiplication<rotgen::matrix<double>>(3, 4, 5,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = (-4) * r + 3*c; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
};
TTS_CASE("Matrix multiplication with itself")
{
test_matrix_multiplication<rotgen::matrix<double>>(3, 3, 3,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3*r - 7*r*c + 14.4; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3*r - 7*r*c + 14.4; });
};
TTS_CASE("Matrix multiplication with identity matrix")
{
test_matrix_multiplication<rotgen::matrix<double>>(4, 4, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r*r*r - c + 12; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = (r==c); });
};
TTS_CASE("Matrix addition")
{
test_matrix_addition<rotgen::matrix<double>>(3, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 9.9*r*r*r - 6*c -12; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.1*r + 4.2*c - 12.3; });
};
TTS_CASE("Matrix addition with zero matrix")
{
test_matrix_addition<rotgen::matrix<double>>(3, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 2*r*r + c + 3.4; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
};
TTS_CASE("Matrix subtraction")
{
test_matrix_substraction<rotgen::matrix<double>>(3, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 6.78*r - 5.2*c - 0.01; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.1*r + 33.456*c*c; });
};
TTS_CASE("Matrix subtraction with zero matrix")
{
test_matrix_substraction<rotgen::matrix<double>>(3, 4,
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r + c*c*c*c - 56.6; },
[](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
};