From b6fcd4b3415b3c1f824a1d065669cfc875f120cc Mon Sep 17 00:00:00 2001 From: Joel Falcou Date: Sun, 28 Sep 2025 16:15:15 +0200 Subject: [PATCH] Implements some missing functions See merge request oss/rotgen!28 --- CMakeLists.txt | 1 + README.md | 6 +- cmake/options.cmake | 2 +- include/rotgen/common/ref.hpp | 83 ++++++++++++++++++-- include/rotgen/concepts.hpp | 28 +++++-- include/rotgen/config.hpp | 5 ++ include/rotgen/dynamic/block.hpp | 11 ++- include/rotgen/dynamic/map.hpp | 57 ++++++++++---- include/rotgen/dynamic/matrix.hpp | 17 ++-- include/rotgen/dynamic/svd.hpp | 36 +++++++++ include/rotgen/extract.hpp | 6 ++ include/rotgen/fixed/block.hpp | 42 ++++++++-- include/rotgen/fixed/map.hpp | 26 +++++-- include/rotgen/fixed/matrix.hpp | 9 ++- include/rotgen/fixed/svd.hpp | 88 +++++++++++++++++++++ include/rotgen/functions.hpp | 77 +++++++++++-------- include/rotgen/impl/map_indirect.hpp | 24 ++++-- include/rotgen/impl/map_model.hpp | 11 ++- include/rotgen/impl/svd.hpp | 68 ++++++++++++++++ include/rotgen/impl/svd_model.hpp | 44 +++++++++++ include/rotgen/operators.hpp | 26 +++++++ include/rotgen/rotgen.hpp | 2 + src/block_model.cpp | 8 +- src/map.cpp | 1 + src/map_indirect.cpp | 24 ++++-- src/map_model.cpp | 32 +++++++- src/svd.cpp | 61 +++++++++++++++ src/svd_model.cpp | 111 +++++++++++++++++++++++++++ test/unit/block/extract.cpp | 5 +- test/unit/common/norms.hpp | 24 ++---- test/unit/functions/svd.cpp | 79 +++++++++++++++++++ test/unit/matrix/constructors.cpp | 25 ++++-- test/unit/matrix/generators.cpp | 10 +-- test/unit/matrix/inverse.cpp | 62 +++++++++++++++ 34 files changed, 972 insertions(+), 139 deletions(-) create mode 100644 include/rotgen/dynamic/svd.hpp create mode 100644 include/rotgen/fixed/svd.hpp create mode 100644 include/rotgen/impl/svd.hpp create mode 100644 include/rotgen/impl/svd_model.hpp create mode 100644 src/svd.cpp create mode 100644 src/svd_model.cpp create mode 100644 test/unit/functions/svd.cpp create mode 100644 test/unit/matrix/inverse.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 476ba59..c942806 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,6 +44,7 @@ if(ROTGEN_FORCE_DYNAMIC) src/matrix.cpp src/block.cpp src/info.cpp + src/svd.cpp src/operators.cpp ) else() diff --git a/README.md b/README.md index 824a9a5..bc0a964 100644 --- a/README.md +++ b/README.md @@ -32,7 +32,7 @@ While the command `cmake -DROTGEN_FORCE_DYNAMIC=ON` will output: ```bash -- ==================== Rotgen Configuration Summary ==================== -- Configuration mode: DYNAMIC --- Reason : No static size options were provided +-- Reason : No static size options were provided -- Preprocessor flags: -DROTGEN_FORCE_DYNAMIC -- ====================================================================== -- Configuring done (0.0s) @@ -64,7 +64,7 @@ Supported scalar types are `float` and `double`. + Data generation interface (i.e `matrix::Zero`, etc...) + Reductions (sum, prod, norms) + Resizing and conservative resizing - + ## Project status ### Current Releases + 08/2025 - Version 0.0.1beta @@ -74,7 +74,7 @@ Supported scalar types are `float` and `double`. + Support `array` and related operations. + Support non-trivial indexing. + Precompile and provide free functions API for common linear algebra solvers. - + ## License **ROTGEN** is licensed under the Boost Software License diff --git a/cmake/options.cmake b/cmake/options.cmake index f1b1153..c47c01f 100644 --- a/cmake/options.cmake +++ b/cmake/options.cmake @@ -58,7 +58,7 @@ function(print_configuration_summary FORCE_DYNAMIC_VAR FORCE_CONFIG_REASON MAX_S if(${FORCE_DYNAMIC_VAR}) message(STATUS " Configuration mode: DYNAMIC") if(${FORCE_CONFIG_REASON}) - message(STATUS " Reason : No static size options were provided") + message(STATUS " Reason : No static size options were provided") endif() else() message(STATUS " Configuration mode: STATIC") diff --git a/include/rotgen/common/ref.hpp b/include/rotgen/common/ref.hpp index d58dd0d..9e21c92 100644 --- a/include/rotgen/common/ref.hpp +++ b/include/rotgen/common/ref.hpp @@ -17,13 +17,14 @@ namespace rotgen { // Primary template: mutable ref - template + template class ref : private map { public: - using parent = map; - using value_type = typename T::value_type; - using rotgen_tag = void; + using parent = map; + using value_type = typename T::value_type; + using rotgen_tag = void; + using rotgen_ref_tag = void; static constexpr int storage_order = T::storage_order; static constexpr int RowsAtCompileTime = T::RowsAtCompileTime; @@ -42,6 +43,7 @@ namespace rotgen using parent::prod; using parent::mean; using parent::trace; + using parent::transpose; using parent::cwiseAbs; using parent::cwiseAbs2; using parent::cwiseInverse; @@ -49,6 +51,7 @@ namespace rotgen using parent::maxCoeff; using parent::minCoeff; using parent::norm; + using parent::normalize; using parent::squaredNorm; using parent::lpNorm; using parent::operator+=; @@ -59,6 +62,10 @@ namespace rotgen using parent::Constant; using parent::Random; using parent::Identity; + using parent::setZero; + using parent::setConstant; + using parent::setRandom; + using parent::setIdentity; using parent::operator=; @@ -74,6 +81,12 @@ namespace rotgen static_assert((O & 1) == storage_order, "ref: Incompatible storage layout"); } + template + ref(block&& b) : parent(b.data(), b.rows(), b.cols(), stride_type{b.outerStride(),b.innerStride()}) + { + static_assert((Ref::storage_order & 1) == storage_order, "ref: Incompatible storage layout"); + } + template ref(block& b) : parent(b.data(), b.rows(), b.cols(), stride_type{b.outerStride(),b.innerStride()}) { @@ -120,6 +133,7 @@ namespace rotgen using parent::prod; using parent::mean; using parent::trace; + using parent::transpose; using parent::cwiseAbs; using parent::cwiseAbs2; using parent::cwiseInverse; @@ -127,6 +141,7 @@ namespace rotgen using parent::maxCoeff; using parent::minCoeff; using parent::norm; + using parent::normalize; using parent::squaredNorm; using parent::lpNorm; using parent::operator+=; @@ -201,18 +216,36 @@ namespace rotgen return lhs.base() + rhs.base(); } + template + auto operator+=(ref lhs, ref rhs) -> decltype(lhs.base() += rhs.base()) + { + return lhs.base() += rhs.base(); + } + template auto operator-(ref lhs, ref rhs) { return lhs.base() - rhs.base(); } + template + auto operator-=(ref lhs, ref rhs) -> decltype(lhs.base() -= rhs.base()) + { + return lhs.base() -= rhs.base(); + } + template auto operator*(ref lhs, ref rhs) { return lhs.base() * rhs.base(); } + template + auto operator*=(ref lhs, ref rhs) -> decltype(lhs.base() *= rhs.base()) + { + return lhs.base() *= rhs.base(); + } + template auto operator*(ref lhs, std::convertible_to auto s) { @@ -297,6 +330,18 @@ namespace rotgen return lhs / s; } + template + auto inverse(ref lhs) -> decltype(lhs.base().inverse()) + { + return lhs.base().inverse(); + } + + template + auto cross(ref lhs, ref rhs) -> decltype(lhs.base().cross(rhs.base())) + { + return lhs.base().cross(rhs.base()); + } + //------------------------------------------------------------------------------------------- // Convert entity/eigen types to a proper ref so we can write less function overloads template struct generalize; @@ -314,10 +359,28 @@ namespace rotgen template struct generalize { static constexpr bool is_const = std::is_const_v; - using base = matrix; + using base = matrix; using type = std::conditional_t, ref>; }; + template + struct generalize> + { + using type = ref; + }; + + template + struct generalize const> + { + using type = ref; + }; + + template + typename T::parent& base_of(T& a) + { + return a.base(); + } + template typename T::parent const& base_of(T const& a) { @@ -335,12 +398,18 @@ namespace rotgen { static constexpr bool is_const = std::is_const_v; using concrete_type = decltype(std::declval().eval()); - using base = matrix; + using base = matrix; using type = std::conditional_t, ref>; }; template - T const& base_of(T const& a) + auto const& base_of(T const& a) + { + return a; + } + + template + auto& base_of(T& a) { return a; } diff --git a/include/rotgen/concepts.hpp b/include/rotgen/concepts.hpp index 3c26e53..5f6003c 100644 --- a/include/rotgen/concepts.hpp +++ b/include/rotgen/concepts.hpp @@ -9,6 +9,24 @@ namespace rotgen::concepts { + //================================================================================================ + //! @brief Check if a type is a Rotgen reference. + //================================================================================================ + template + concept reference = requires { typename std::remove_cvref_t::rotgen_ref_tag; }; + + //================================================================================================ + //! @brief Check if a type is a vector. + //================================================================================================ + template + concept vector = T::IsVectorAtCompileTime; + + //================================================================================================ + //! @brief Check if a type is a 3D vector. + //================================================================================================ + template + concept vectorND = vector && (T::RowsAtCompileTime == N || T::ColsAtCompileTime == N); + //================================================================================================ //! @brief Check if a type is a ROTGEN type. //================================================================================================ @@ -25,12 +43,12 @@ namespace rotgen::concepts template concept eigen_compatible = requires(T const& a) { - typename T::Scalar; + typename std::remove_cvref_t::Scalar; - { T::RowsAtCompileTime } -> std::convertible_to; - { T::ColsAtCompileTime } -> std::convertible_to; - { T::MaxRowsAtCompileTime } -> std::convertible_to; - { T::MaxColsAtCompileTime } -> std::convertible_to; + { std::remove_cvref_t::RowsAtCompileTime } -> std::convertible_to; + { std::remove_cvref_t::ColsAtCompileTime } -> std::convertible_to; + { std::remove_cvref_t::MaxRowsAtCompileTime } -> std::convertible_to; + { std::remove_cvref_t::MaxColsAtCompileTime } -> std::convertible_to; { a.eval() }; }; diff --git a/include/rotgen/config.hpp b/include/rotgen/config.hpp index eb91cae..927990a 100644 --- a/include/rotgen/config.hpp +++ b/include/rotgen/config.hpp @@ -22,6 +22,11 @@ namespace rotgen return w == -1 ? 0 : w / 8; }(); + inline constexpr int ComputeFullU = 0x04; + inline constexpr int ComputeThinU = 0x08; + inline constexpr int ComputeFullV = 0x10; + inline constexpr int ComputeThinV = 0x20; + inline constexpr int Dynamic = -1; inline constexpr int Infinity = -1; inline constexpr int AutoAlign = 0; diff --git a/include/rotgen/dynamic/block.hpp b/include/rotgen/dynamic/block.hpp index 24d7146..370d7b9 100644 --- a/include/rotgen/dynamic/block.hpp +++ b/include/rotgen/dynamic/block.hpp @@ -15,6 +15,9 @@ namespace rotgen { + template + class ref; + template class block : public find_block { @@ -55,7 +58,7 @@ namespace rotgen block(Ref const& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) requires( !requires{typename Ref::rotgen_block_tag; } && is_immutable) - : parent(r,i0,j0,ni,nj) + : parent(r.base(),i0,j0,ni,nj) {} block(Ref const& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) @@ -65,7 +68,7 @@ namespace rotgen block(Ref const& r, std::size_t i0, std::size_t j0) requires(!requires{typename Ref::rotgen_block_tag; } && Rows != -1 && Cols != -1 && is_immutable) - : parent(r,i0,j0,Rows,Cols) + : parent(r.base(),i0,j0,Rows,Cols) {} block(Ref const& r, std::size_t i0, std::size_t j0) @@ -75,7 +78,7 @@ namespace rotgen block(Ref& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) requires( !requires{typename Ref::rotgen_block_tag; } && !is_immutable) - : parent(r,i0,j0,ni,nj) + : parent(r.base(),i0,j0,ni,nj) {} block(Ref& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) @@ -85,7 +88,7 @@ namespace rotgen block(Ref& r, std::size_t i0, std::size_t j0) requires(!requires{typename Ref::rotgen_block_tag; } && Rows != -1 && Cols != -1 && !is_immutable) - : parent(r,i0,j0,Rows,Cols) + : parent(r.base(),i0,j0,Rows,Cols) {} block(Ref& r, std::size_t i0, std::size_t j0) diff --git a/include/rotgen/dynamic/map.hpp b/include/rotgen/dynamic/map.hpp index fe45aed..6f19b7f 100644 --- a/include/rotgen/dynamic/map.hpp +++ b/include/rotgen/dynamic/map.hpp @@ -136,13 +136,33 @@ namespace rotgen concrete_type cwiseAbs2() const { return concrete_type(base().cwiseAbs2()); } concrete_type cwiseInverse() const { return concrete_type(base().cwiseInverse()); } concrete_type cwiseSqrt() const { return concrete_type(base().cwiseSqrt()); } - concrete_type cwiseMin (map const& rhs) const { return concrete_type(base().cwiseMin(rhs.base())); } - concrete_type cwiseMax (map const& rhs) const { return concrete_type(base().cwiseMax(rhs.base())); } - concrete_type cwiseQuotient (map const& rhs) const { return concrete_type(base().cwiseQuotient(rhs.base())); } + concrete_type cwiseMin(map const& rhs) const { return concrete_type(base().cwiseMin(rhs.base())); } + concrete_type cwiseMax(map const& rhs) const { return concrete_type(base().cwiseMax(rhs.base())); } + concrete_type cwiseQuotient(map const& rhs) const { return concrete_type(base().cwiseQuotient(rhs.base())); } concrete_type cwiseProduct(map const& rhs) const { return concrete_type(base().cwiseProduct(rhs.base())); } concrete_type cwiseMin(value_type s) const { return concrete_type(base().cwiseMin(s)); } concrete_type cwiseMax(value_type s) const { return concrete_type(base().cwiseMax(s)); } + concrete_type inverse() const { return concrete_type(base().inverse()); } + + concrete_type cross(map const& other) const + { + concrete_type that; + if constexpr(RowsAtCompileTime==-1) + { + that(0,0) = (*this)(1,0) * other(2,0) - (*this)(2,0) * other(1,0); + that(1,0) = (*this)(2,0) * other(0,0) - (*this)(0,0) * other(2,0); + that(2,0) = (*this)(0,0) * other(1,0) - (*this)(1,0) * other(0,0); + } + else + { + that(0,0) = (*this)(0,1) * other(0,2) - (*this)(0,2) * other(0,1); + that(0,1) = (*this)(0,2) * other(0,0) - (*this)(0,0) * other(0,2); + that(0,2) = (*this)(0,0) * other(0,1) - (*this)(0,1) * other(0,0); + } + return that; + } + void normalize() requires(!is_immutable && IsVectorAtCompileTime) { parent::normalize(); @@ -150,15 +170,17 @@ namespace rotgen void transposeInPlace() requires(!is_immutable) { parent::transposeInPlace(); } void adjointInPlace() requires(!is_immutable) { parent::adjointInPlace(); } - map& operator+=(map const& rhs) requires(!is_immutable) + template + map& operator+=(map const& rhs) requires(!is_immutable) { - base() += static_cast(rhs); + base() += rhs.base(); return *this; } - map& operator-=(map const& rhs) requires(!is_immutable) + template + map& operator-=(map const& rhs) requires(!is_immutable) { - base() -= static_cast(rhs); + base() -= rhs.base(); return *this; } @@ -167,9 +189,10 @@ namespace rotgen return concrete_type(base().operator-()); } - map& operator*=(map const& rhs) requires(!is_immutable) + template + map& operator*=(map const& rhs) requires(!is_immutable) { - base() *= static_cast(rhs); + base() *= rhs.base(); return *this; } @@ -297,22 +320,26 @@ namespace rotgen template typename map::concrete_type operator+(map const& lhs, map const& rhs) { - using concrete_type = typename map::concrete_type; - return concrete_type(lhs.base().add(map(rhs))); + // REPLICATE TODO + using map_type = map; + using concrete_type = typename map_type::concrete_type; + return concrete_type(map_type(lhs).base().add(map_type(rhs))); } template typename map::concrete_type operator-(map const& lhs, map const& rhs) { - using concrete_type = typename map::concrete_type; - return concrete_type(lhs.base().sub(map(rhs))); + using map_type = map; + using concrete_type = typename map_type::concrete_type; + return concrete_type(map_type(lhs).base().sub(map_type(rhs))); } template typename map::concrete_type operator*(map const& lhs, map const& rhs) { - using concrete_type = typename map::concrete_type; - return concrete_type(lhs.base().mul(map(rhs))); + using map_type = map; + using concrete_type = typename map_type::concrete_type; + return concrete_type(map_type(lhs).base().mul(map_type(rhs))); } template diff --git a/include/rotgen/dynamic/matrix.hpp b/include/rotgen/dynamic/matrix.hpp index 3274632..2ed3e52 100644 --- a/include/rotgen/dynamic/matrix.hpp +++ b/include/rotgen/dynamic/matrix.hpp @@ -29,6 +29,7 @@ namespace rotgen static constexpr auto storage_order = Opts & 1; static constexpr Index RowsAtCompileTime = Rows; static constexpr Index ColsAtCompileTime = Cols; + static constexpr bool IsCompileTimeSized = Rows != -1 && Cols != -1; static constexpr bool IsVectorAtCompileTime = (RowsAtCompileTime == 1) || (ColsAtCompileTime == 1); static constexpr int Options = Opts; static constexpr bool IsRowMajor = (Opts & RowMajor) == RowMajor; @@ -37,7 +38,8 @@ namespace rotgen matrix() : parent(Rows==-1?0:Rows,Cols==-1?0:Cols) {} - matrix(Index r, Index c) : parent(r, c) + matrix(Index r, Index c) requires(!IsCompileTimeSized) + : parent(r, c) { if constexpr(Rows != -1) assert(r == Rows && "Mismatched between dynamic and static row size"); if constexpr(Cols != -1) assert(c == Cols && "Mismatched between dynamic and static column size"); @@ -49,6 +51,11 @@ namespace rotgen matrix(Scalar v) requires(Rows == 1 && Cols == 1) : parent(1,1,{v}) {} + matrix(Scalar v0, Scalar v1, auto... vs) + requires((Rows == (2+sizeof...(vs)) && Cols == 1) || (Rows == 1 && Cols == (2+sizeof...(vs)))) + : parent(Rows,Cols,{v0,v1,static_cast(vs)...}) + {} + matrix(parent const& base) : parent(base) {} matrix(std::initializer_list> init) : parent(init) @@ -129,9 +136,9 @@ namespace rotgen return matrix(base().normalized()); } - matrix transpose() const + matrix transpose() const { - return matrix(base().transpose()); + return matrix(base().transpose()); } matrix conjugate() const @@ -139,9 +146,9 @@ namespace rotgen return matrix(base().conjugate()); } - matrix adjoint() const + matrix adjoint() const { - return matrix(base().adjoint()); + return matrix(base().adjoint()); } void normalize() requires(IsVectorAtCompileTime) diff --git a/include/rotgen/dynamic/svd.hpp b/include/rotgen/dynamic/svd.hpp new file mode 100644 index 0000000..e40bd4b --- /dev/null +++ b/include/rotgen/dynamic/svd.hpp @@ -0,0 +1,36 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#pragma once + +#include + +namespace rotgen +{ + template + struct svd : find_svd + { + using parent = find_svd; + + using m_type = matrix; + using d_type = matrix; + + svd(M const& m, int options = ComputeThinU | ComputeThinV) : parent(m,options) {} + + int rank() const { return parent::rank(); } + + m_type U() const { return parent::U(); } + m_type D() const { return parent::D(); } + m_type V() const { return parent::V(); } + d_type singular_values() const { return parent::singular_values(); } + + m_type U(int r) const { return parent::U(r); } + m_type D(int r) const { return parent::D(r); } + m_type V(int r) const { return parent::V(r); } + m_type singular_values(int r) const { return parent::singular_values(r); } + }; +} \ No newline at end of file diff --git a/include/rotgen/extract.hpp b/include/rotgen/extract.hpp index 77c3535..02186a4 100644 --- a/include/rotgen/extract.hpp +++ b/include/rotgen/extract.hpp @@ -30,6 +30,12 @@ namespace rotgen detail::validate_extract(e,i0,j0,ni,nj); return block(e, i0, j0, ni, nj); } + template + auto extract(Entity const& e, Index i0, Index j0, Index ni, Index nj) + { + detail::validate_extract(e,i0,j0,ni,nj); + return block(e, i0, j0, ni, nj); + } template requires(NI!=-1 && NJ!=-1) diff --git a/include/rotgen/fixed/block.hpp b/include/rotgen/fixed/block.hpp index 8293af0..d16b6df 100644 --- a/include/rotgen/fixed/block.hpp +++ b/include/rotgen/fixed/block.hpp @@ -12,6 +12,9 @@ namespace rotgen { + template + class ref; + namespace detail { template @@ -25,6 +28,12 @@ namespace rotgen >; }; + template + struct compute_block_type,Rows,Cols,Inner,IsConst> + : compute_block_type::parent,Rows,Cols,Inner,IsConst> + { + }; + template using block_type = typename compute_block_type::type; } @@ -46,7 +55,7 @@ namespace rotgen using value_type = typename parent::value_type; using Index = typename parent::Index; - static constexpr int storage_order = (ForceStorageOrder == -1) ? Ref::storage_order : ForceStorageOrder; + static constexpr int storage_order = (ForceStorageOrder == -1) ? Ref::storage_order : ForceStorageOrder; static constexpr bool is_immutable = std::is_const_v; static constexpr bool IsRowMajor = parent::IsRowMajor; @@ -57,7 +66,6 @@ namespace rotgen template using as_concrete_type = as_concrete_t; - static constexpr auto Flags = Ref::Flags; static constexpr int Options = Ref::Options; static constexpr int RowsAtCompileTime = Rows; static constexpr int ColsAtCompileTime = Cols; @@ -72,26 +80,48 @@ namespace rotgen block& operator=(const block&) = default; block& operator=(block&&) = default; + // Constructs from regular Ref block(Ref const& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) - requires(is_immutable) + requires(is_immutable && !concepts::reference) : parent(r.base(),i0,j0,ni,nj) {} block(Ref const& r, std::size_t i0, std::size_t j0) - requires(Rows != -1 && Cols != -1 && is_immutable) + requires(Rows != -1 && Cols != -1 && is_immutable && !concepts::reference) : parent(r.base(),i0,j0,Rows,Cols) {} block(Ref& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) - requires(!is_immutable) + requires(!is_immutable && !concepts::reference) : parent(r.base(),i0,j0,ni,nj) {} block(Ref& r, std::size_t i0, std::size_t j0) - requires(Rows != -1 && Cols != -1 && !is_immutable) + requires(Rows != -1 && Cols != -1 && !is_immutable && !concepts::reference) : parent(r.base(),i0,j0,Rows,Cols) {} + // Constructs from rotgen::ref + block(Ref const& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) + requires(is_immutable && concepts::reference) + : parent(r.base().base(),i0,j0,ni,nj) + {} + + block(Ref const& r, std::size_t i0, std::size_t j0) + requires(Rows != -1 && Cols != -1 && is_immutable && concepts::reference) + : parent(r.base().base(),i0,j0,Rows,Cols) + {} + + block(Ref& r, std::size_t i0, std::size_t j0, std::size_t ni, std::size_t nj) + requires(!is_immutable && concepts::reference) + : parent(r.base().base(),i0,j0,ni,nj) + {} + + block(Ref& r, std::size_t i0, std::size_t j0) + requires(Rows != -1 && Cols != -1 && !is_immutable && concepts::reference) + : parent(r.base().base(),i0,j0,Rows,Cols) + {} + template block(block const& other) : parent(other.base()) {} diff --git a/include/rotgen/fixed/map.hpp b/include/rotgen/fixed/map.hpp index 619734a..e099cca 100644 --- a/include/rotgen/fixed/map.hpp +++ b/include/rotgen/fixed/map.hpp @@ -39,7 +39,6 @@ namespace rotgen using value_type = typename std::remove_const_t::value_type; using concrete_type = typename std::remove_const_t::concrete_type; - static constexpr auto Flags = Ref::Flags; static constexpr Index RowsAtCompileTime = Ref::RowsAtCompileTime; static constexpr Index ColsAtCompileTime = Ref::ColsAtCompileTime; static constexpr bool IsVectorAtCompileTime = Ref::IsVectorAtCompileTime; @@ -131,21 +130,24 @@ namespace rotgen value_type operator()(Index i) const requires(IsVectorAtCompileTime) { return base().data()[i]; } value_type operator[](Index i) const requires(IsVectorAtCompileTime) { return (*this)(i); } - map& operator+=(map const& rhs) + template + map& operator+=(map const& rhs) requires(!is_immutable) { base() += rhs.base(); return *this; } - map& operator-=(map const& rhs) + template + map& operator-=(map const& rhs) requires(!is_immutable) { base() -= rhs.base(); return *this; } - map& operator*=(map const& rhs) + template + map& operator*=(map const& rhs) requires(!is_immutable) { - base() *= rhs; + base() *= rhs.base(); return *this; } @@ -221,6 +223,18 @@ namespace rotgen else return base().cwiseSqrt(); } + auto cross(map const& rhs) const + { + if constexpr(!use_expression_templates) return concrete_type{parent::cross(rhs.base())}; + else return base().cross(rhs.base()); + } + + auto inverse() const + { + if constexpr(use_expression_templates) return base().inverse(); + else return as_concrete_type(base().inverse()); + } + auto normalized() const requires(IsVectorAtCompileTime) { if constexpr(use_expression_templates) return base().normalized(); @@ -253,7 +267,7 @@ namespace rotgen void transposeInPlace() { base().transposeInPlace(); } void adjointInPlace() { base().adjointInPlace(); } - auto qr_solve(map const& rhs) const + auto qr_solve(auto const& rhs) const { return concrete_type(base().colPivHouseholderQr().solve(rhs.base())); }; diff --git a/include/rotgen/fixed/matrix.hpp b/include/rotgen/fixed/matrix.hpp index 6ceccfe..35b8447 100644 --- a/include/rotgen/fixed/matrix.hpp +++ b/include/rotgen/fixed/matrix.hpp @@ -38,9 +38,9 @@ namespace rotgen using concrete_type = matrix; using concrete_dynamic_type = matrix; - static constexpr auto Flags = parent::Flags; static constexpr int RowsAtCompileTime = Rows; static constexpr int ColsAtCompileTime = Cols; + static constexpr bool IsCompileTimeSized = Rows != -1 && Cols != -1; static constexpr bool IsVectorAtCompileTime = (RowsAtCompileTime == 1) || (ColsAtCompileTime == 1); static constexpr int Options = parent::Options; static constexpr bool IsRowMajor = parent::IsRowMajor; @@ -55,7 +55,7 @@ namespace rotgen matrix() requires(has_static_storage) {} matrix() requires(!has_static_storage) : parent(Rows > 0 ? Rows : 0, Cols > 0 ? Cols : 0){} - matrix(Index r, Index c) : parent(r, c) {} + matrix(Index r, Index c) requires(!IsCompileTimeSized) : parent(r, c) {} matrix(const matrix& other) = default; matrix(matrix&& other) = default; @@ -94,6 +94,11 @@ namespace rotgen (*this)(i) = first[i]; } + matrix(Scalar v0, Scalar v1, auto... vs) + requires((Rows == (2+sizeof...(vs)) && Cols == 1) || (Rows == 1 && Cols == (2+sizeof...(vs)))) + : matrix({v0,v1,static_cast(vs)...}) + {} + matrix(concepts::entity auto const& other) : parent(other.base()) {} diff --git a/include/rotgen/fixed/svd.hpp b/include/rotgen/fixed/svd.hpp new file mode 100644 index 0000000..ff9cb0d --- /dev/null +++ b/include/rotgen/fixed/svd.hpp @@ -0,0 +1,88 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#pragma once + +#include +#include + +namespace rotgen +{ + template + struct svd + { + using base = Eigen::JacobiSVD; + using u_type = typename base::MatrixUType; + using v_type = typename base::MatrixVType; + using d_type = typename base::SingularValuesType; + + svd(M const& m, int options = ComputeThinU | ComputeThinV) : svd_(m.base(),options) {} + + int rank() const { return svd_.rank(); } + + auto singular_values() const + { + if constexpr(!use_expression_templates) return as_concrete_t{svd_.singularValues()}; + else return svd_.singularValues(); + } + + auto U() const + { + if constexpr(!use_expression_templates) return as_concrete_t{svd_.matrixU()}; + else return svd_.matrixU(); + } + + auto V() const + { + if constexpr(!use_expression_templates) return as_concrete_t{svd_.matrixV()}; + else return svd_.matrixV(); + } + + auto D() const + { + auto d = svd_.singularValues().asDiagonal(); + if constexpr(!use_expression_templates) + return as_concrete_t{d.toDenseMatrix ()}; + else + return d; + } + + + auto singular_values(int r) const + { + auto that = svd_.singularValues().head(r); + if constexpr(!use_expression_templates) return as_concrete_t{that}; + else return svd_.singularValues(); + } + + auto U(int r) const + { + auto that = svd_.matrixU().leftCols(r); + if constexpr(!use_expression_templates) return as_concrete_t{that}; + else return that; + } + + auto V(int r) const + { + auto that = svd_.matrixV().leftCols(r); + if constexpr(!use_expression_templates) return as_concrete_t{that}; + else return that; + } + + auto D(int r) const + { + auto d = svd_.singularValues().head(r).asDiagonal(); + if constexpr(!use_expression_templates) + return as_concrete_t{d.toDenseMatrix ()}; + else + return d; + } + + private: + base svd_; + }; +} \ No newline at end of file diff --git a/include/rotgen/functions.hpp b/include/rotgen/functions.hpp index ba6a81e..54c58d0 100644 --- a/include/rotgen/functions.hpp +++ b/include/rotgen/functions.hpp @@ -12,47 +12,43 @@ namespace rotgen //----------------------------------------------------------------------------------------------- // Infos & Shape //----------------------------------------------------------------------------------------------- - std::size_t rows(concepts::entity auto const& arg) { return arg.rows(); } - std::size_t cols(concepts::entity auto const& arg) { return arg.cols(); } - std::size_t size(concepts::entity auto const& arg) { return arg.size(); } + std::size_t rows(auto const& m) requires(requires{ m.rows(); }){ return m.rows(); } + std::size_t cols(auto const& m) requires(requires{ m.cols(); }){ return m.cols(); } + std::size_t size(auto const& m) requires(requires{ m.size(); }){ return m.size(); } - template - void resize(E& arg, int sz) requires requires{arg.resize(sz);} - { - arg.resize(sz); - } + void resize(auto& a, int sz) requires requires{a.resize(sz);} { a.resize(sz); } - template - void resize(E& arg, int new_rows, int new_cols) + template + void resize(E& a, int r, int c) requires(E::RowsAtCompileTime == -1 && E::ColsAtCompileTime == -1) { - arg.resize(new_rows, new_cols); + a.resize(r, c); } - template - void conservativeResize(E& arg, int sz) requires requires{arg.conservativeResize(sz);} + template + void conservativeResize(E& a, int sz) requires requires{a.conservativeResize(sz);} { - arg.conservativeResize(sz); + a.conservativeResize(sz); } - template - void conservativeResize(E& arg, int new_rows, int new_cols) + template + void conservativeResize(E& a, int r, int c) requires(E::RowsAtCompileTime == -1 && E::ColsAtCompileTime == -1) { - arg.conservativeResize(new_rows, new_cols); + a.conservativeResize(r, c); } //----------------------------------------------------------------------------------------------- // Global operations //----------------------------------------------------------------------------------------------- - decltype(auto) normalized(concepts::entity auto const& arg) { return arg.normalized(); } - decltype(auto) transpose (concepts::entity auto const& arg) { return arg.transpose(); } - decltype(auto) conjugate (concepts::entity auto const& arg) { return arg.conjugate(); } - decltype(auto) adjoint (concepts::entity auto const& arg) { return arg.adjoint(); } + decltype(auto) normalized (auto const& m) requires(requires{ m.normalized(); }) { return m.normalized(); } + decltype(auto) transpose (auto const& m) requires(requires{ m.transpose(); }) { return m.transpose(); } + decltype(auto) conjugate (auto const& m) requires(requires{ m.conjugate(); }) { return m.conjugate(); } + decltype(auto) adjoint (auto const& m) requires(requires{ m.adjoint(); }) { return m.adjoint(); } - void normalize(concepts::entity auto& arg) { arg.normalize(); } - void transposeInPlace(concepts::entity auto& arg) { arg.transposeInPlace(); } - void adjointInPlace(concepts::entity auto& arg) { arg.adjointInPlace(); } + void normalize(auto& a) requires(requires{ a.normalize(); }) { a.normalize(); } + void transposeInPlace(auto& a) requires(requires{ a.transposeInPlace(); }) { a.transposeInPlace(); } + void adjointInPlace(auto& a) requires(requires{ a.adjointInPlace(); }) { a.adjointInPlace(); } //----------------------------------------------------------------------------------------------- // Component-wise functions @@ -146,13 +142,8 @@ namespace rotgen auto prod(concepts::entity auto const& arg) { return arg.prod(); } auto mean(concepts::entity auto const& arg) { return arg.mean(); } - auto maxCoeff(auto const& arg) - requires( requires{ arg.maxCoeff(); } ) - { - return arg.maxCoeff(); - } - - auto minCoeff(concepts::entity auto const& arg) { return arg.minCoeff(); } + auto maxCoeff(auto const& arg) requires( requires{ arg.maxCoeff(); } ) { return arg.maxCoeff(); } + auto minCoeff(auto const& arg) requires( requires{ arg.minCoeff(); } ) { return arg.minCoeff(); } template auto maxCoeff(concepts::entity auto const& arg, IndexType* row, IndexType* col) @@ -166,8 +157,8 @@ namespace rotgen return arg.minCoeff(row, col); } - template - auto lpNorm(concepts::entity auto const& arg) + template + auto lpNorm(T const& arg) requires( requires{arg.template lpNorm

();} ) { static_assert(P == 1 || P == 2 || P == Infinity); return arg.template lpNorm

(); @@ -316,4 +307,24 @@ namespace rotgen { return T::Constant(r,c,v); } + + //----------------------------------------------------------------------------------------------- + // Others + // TODO: Adapt other functions ot behave as inverse and limit code in non-ref/non-map containers + //----------------------------------------------------------------------------------------------- + template auto inverse(A const& a) + { + if constexpr(!use_expression_templates) return inverse(generalize_t(a)); + else if constexpr(requires{typename A::rotgen_tag;}) return base_of(a).inverse(); + else return a.inverse(); + } + + template L, concepts::vectorND<3> R> auto cross(L const& l, R const& r) + { + if constexpr(!use_expression_templates) + return cross(generalize_t(l),generalize_t(r)); + else if constexpr(requires{typename L::rotgen_tag;} || requires{typename R::rotgen_tag;} ) + return base_of(l).cross(base_of(r)); + else return l.cross(r); + } } \ No newline at end of file diff --git a/include/rotgen/impl/map_indirect.hpp b/include/rotgen/impl/map_indirect.hpp index a11c1fe..b7dba56 100644 --- a/include/rotgen/impl/map_indirect.hpp +++ b/include/rotgen/impl/map_indirect.hpp @@ -1,17 +1,21 @@ #define SIZE 64 #define TYPE double - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) #include #undef CLASSNAME #undef SOURCENAME + #undef CLASSCONSTNAME - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) #include #undef CLASSNAME #undef SOURCENAME + #undef CLASSCONSTNAME #undef SIZE #undef TYPE @@ -19,17 +23,21 @@ #define SIZE 32 #define TYPE float - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) #include #undef CLASSNAME #undef SOURCENAME + #undef CLASSCONSTNAME - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) #include #undef CLASSNAME #undef SOURCENAME + #undef CLASSCONSTNAME #undef SIZE #undef TYPE diff --git a/include/rotgen/impl/map_model.hpp b/include/rotgen/impl/map_model.hpp index 7c2c74c..2f2015e 100644 --- a/include/rotgen/impl/map_model.hpp +++ b/include/rotgen/impl/map_model.hpp @@ -18,15 +18,13 @@ class ROTGEN_EXPORT CLASSNAME CLASSNAME(TYPE CONST* ptr, Index r, Index c, stride s); CLASSNAME(CLASSNAME const& other); + CLASSNAME(CLASSNAME&&) noexcept; #if !defined(USE_CONST) CLASSNAME& operator=(CLASSNAME const& other); + CLASSNAME& operator=(CLASSNAME&& other); #endif - CLASSNAME(CLASSNAME&&) noexcept; - - CLASSNAME& operator=(CLASSNAME&&) noexcept; - ~CLASSNAME(); Index rows() const; @@ -85,8 +83,11 @@ class ROTGEN_EXPORT CLASSNAME #if !defined(USE_CONST) CLASSNAME& operator+=(CLASSNAME const& rhs); + CLASSNAME& operator+=(CLASSCONSTNAME const& rhs); CLASSNAME& operator-=(CLASSNAME const& rhs); + CLASSNAME& operator-=(CLASSCONSTNAME const& rhs); CLASSNAME& operator*=(CLASSNAME const& rhs); + CLASSNAME& operator*=(CLASSCONSTNAME const& rhs); CLASSNAME& operator*=(TYPE d); CLASSNAME& operator/=(TYPE d); #endif @@ -98,6 +99,8 @@ class ROTGEN_EXPORT CLASSNAME SOURCENAME mul(TYPE s) const; SOURCENAME div(TYPE s) const; + SOURCENAME inverse() const; + friend ROTGEN_EXPORT std::ostream& operator<<(std::ostream&,CLASSNAME const&); const TYPE* data() const; diff --git a/include/rotgen/impl/svd.hpp b/include/rotgen/impl/svd.hpp new file mode 100644 index 0000000..ad95258 --- /dev/null +++ b/include/rotgen/impl/svd.hpp @@ -0,0 +1,68 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#pragma once + +#include +#include +#include +#include + +namespace rotgen +{ + class matrix_impl64_row; + class matrix_impl64_col; + class matrix_impl32_row; + class matrix_impl32_col; + + #define SIZE 64 + #define TYPE double + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #include + #undef CLASSNAME + #undef SOURCENAME + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #include + #undef CLASSNAME + #undef SOURCENAME + + #undef SIZE + #undef TYPE + + #define SIZE 32 + #define TYPE float + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #include + #undef CLASSNAME + #undef SOURCENAME + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #include + #undef CLASSNAME + #undef SOURCENAME + + #undef SIZE + #undef TYPE + + + template struct find_svd_impl; + + template<> struct find_svd_impl { using type = svd_impl32_col; }; + template<> struct find_svd_impl { using type = svd_impl32_row; }; + template<> struct find_svd_impl { using type = svd_impl64_col; }; + template<> struct find_svd_impl { using type = svd_impl64_row; }; + + template + using find_svd = typename find_svd_impl::type; +} \ No newline at end of file diff --git a/include/rotgen/impl/svd_model.hpp b/include/rotgen/impl/svd_model.hpp new file mode 100644 index 0000000..83f428f --- /dev/null +++ b/include/rotgen/impl/svd_model.hpp @@ -0,0 +1,44 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== + +//================================================================================================== +/* + This file is a X-File to generate various svd_impl_* declarations variant +*/ +//================================================================================================== +class ROTGEN_EXPORT CLASSNAME +{ + public: + CLASSNAME(SOURCENAME const& m, int options); + CLASSNAME(CLASSNAME const& other); + CLASSNAME(CLASSNAME&&) noexcept; + + CLASSNAME& operator=(CLASSNAME const& other); + CLASSNAME& operator=(CLASSNAME&&) noexcept; + + ~CLASSNAME(); + + int rank() const; + SOURCENAME U() const; + SOURCENAME D() const; + SOURCENAME singular_values() const; + SOURCENAME V() const; + + SOURCENAME U(int) const; + SOURCENAME D(int) const; + SOURCENAME singular_values(int) const; + SOURCENAME V(int) const; + + private: + struct payload; + std::unique_ptr storage_; + + public: + std::unique_ptr& storage() { return storage_; } + std::unique_ptr const& storage() const { return storage_; } +}; \ No newline at end of file diff --git a/include/rotgen/operators.hpp b/include/rotgen/operators.hpp index cc8f5d8..2ff082b 100644 --- a/include/rotgen/operators.hpp +++ b/include/rotgen/operators.hpp @@ -68,4 +68,30 @@ namespace rotgen if constexpr(!use_expression_templates) return generalize_t(a) / b; else return base_of(a) / b; } + + //------------------------------------------------------------------------------------------------ + // Compounds operators across types + template + auto operator+=(A& a, B const& b) + requires(concepts::entity || concepts::entity) + { + if constexpr(!use_expression_templates) return generalize_t(a) += generalize_t(b); + else return base_of(a) += base_of(b); + } + + template + auto operator-=(A& a, B const& b) + requires(concepts::entity || concepts::entity) + { + if constexpr(!use_expression_templates) return generalize_t(a) -= generalize_t(b); + else return base_of(a) -= base_of(b); + } + + template + auto operator*=(A& a, B const& b) + requires(concepts::entity || concepts::entity) + { + if constexpr(!use_expression_templates) return generalize_t(a) *= generalize_t(b); + else return base_of(a) *= base_of(b); + } } \ No newline at end of file diff --git a/include/rotgen/rotgen.hpp b/include/rotgen/rotgen.hpp index 9207428..1787d10 100644 --- a/include/rotgen/rotgen.hpp +++ b/include/rotgen/rotgen.hpp @@ -16,10 +16,12 @@ #include #include #include +#include #else #include #include #include +#include #endif #include diff --git a/src/block_model.cpp b/src/block_model.cpp index 7bbf3a2..5c02ce4 100644 --- a/src/block_model.cpp +++ b/src/block_model.cpp @@ -11,8 +11,6 @@ This file is a X-File to generate various block_impl_* definitions variant */ //================================================================================================== -#define STR(text) STR_I(text) -#define STR_I(...) #__VA_ARGS__ //================================================================================================== // Internal payload @@ -107,7 +105,11 @@ struct CLASSNAME::payload return *this; } - CLASSNAME& CLASSNAME::operator=(CLASSNAME&&) noexcept = default; + CLASSNAME& CLASSNAME::operator=(CLASSNAME&& o) noexcept + { + if (this != &o) { storage_->data = std::move(o.storage_->data); } + return *this; + } #endif diff --git a/src/map.cpp b/src/map.cpp index f129684..c7d5d41 100644 --- a/src/map.cpp +++ b/src/map.cpp @@ -9,6 +9,7 @@ #include #include #include +#include namespace rotgen { diff --git a/src/map_indirect.cpp b/src/map_indirect.cpp index d271163..05b153c 100644 --- a/src/map_indirect.cpp +++ b/src/map_indirect.cpp @@ -2,19 +2,23 @@ #define TYPE double #define STORAGE_ORDER Eigen::ColMajor - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) #include "map_model.cpp" #undef CLASSNAME + #undef CLASSCONSTNAME #undef SOURCENAME #undef STORAGE_ORDER #define STORAGE_ORDER Eigen::RowMajor - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) #include "map_model.cpp" #undef CLASSNAME #undef SOURCENAME + #undef CLASSCONSTNAME #undef STORAGE_ORDER #undef SIZE @@ -24,18 +28,22 @@ #define TYPE float #define STORAGE_ORDER Eigen::ColMajor - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_col) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) #include "map_model.cpp" #undef CLASSNAME + #undef CLASSCONSTNAME #undef SOURCENAME #undef STORAGE_ORDER #define STORAGE_ORDER Eigen::RowMajor - #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) - #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #define CLASSNAME ROTGEN_MATRIX_NAME(BASENAME,SIZE,_row) + #define CLASSCONSTNAME ROTGEN_MATRIX_NAME(map_const_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) #include "map_model.cpp" #undef CLASSNAME + #undef CLASSCONSTNAME #undef SOURCENAME #undef STORAGE_ORDER diff --git a/src/map_model.cpp b/src/map_model.cpp index c2e1df6..bec9a87 100644 --- a/src/map_model.cpp +++ b/src/map_model.cpp @@ -35,10 +35,15 @@ if (this != &o) storage_->data = o.storage_->data; return *this; } + + CLASSNAME& CLASSNAME::operator=(CLASSNAME&& o) + { + if (this != &o) storage_->data = std::move(o.storage_->data); + return *this; + } #endif CLASSNAME::CLASSNAME(CLASSNAME&&) noexcept = default; - CLASSNAME& CLASSNAME::operator=(CLASSNAME&&) noexcept = default; CLASSNAME::~CLASSNAME() = default; @@ -160,6 +165,13 @@ return result; } + SOURCENAME CLASSNAME::inverse() const + { + SOURCENAME result; + result.storage()->assign(storage_->data.inverse().eval()); + return result; + } + #if !defined(USE_CONST) void CLASSNAME::normalize() { @@ -267,18 +279,36 @@ return *this; } + CLASSNAME& CLASSNAME::operator+=(CLASSCONSTNAME const& rhs) + { + storage_->data += rhs.storage()->data; + return *this; + } + CLASSNAME& CLASSNAME::operator-=(CLASSNAME const& rhs) { storage_->data -= rhs.storage_->data; return *this; } + CLASSNAME& CLASSNAME::operator-=(CLASSCONSTNAME const& rhs) + { + storage_->data -= rhs.storage()->data; + return *this; + } + CLASSNAME& CLASSNAME::operator*=(CLASSNAME const& rhs) { storage_->data *= rhs.storage_->data; return *this; } + CLASSNAME& CLASSNAME::operator*=(CLASSCONSTNAME const& rhs) + { + storage_->data *= rhs.storage()->data; + return *this; + } + CLASSNAME& CLASSNAME::operator*=(TYPE s) { storage_->data *= s; diff --git a/src/svd.cpp b/src/svd.cpp new file mode 100644 index 0000000..d85ff3e --- /dev/null +++ b/src/svd.cpp @@ -0,0 +1,61 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#include +#include +#include +#include +#include +#include +#include + +namespace rotgen +{ + #define SIZE 64 + #define TYPE double + #define STORAGE_ORDER Eigen::ColMajor + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #include "svd_model.cpp" + #undef CLASSNAME + #undef SOURCENAME + #undef STORAGE_ORDER + + #define STORAGE_ORDER Eigen::RowMajor + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #include "svd_model.cpp" + #undef CLASSNAME + #undef SOURCENAME + #undef STORAGE_ORDER + + #undef SIZE + #undef TYPE + + #define SIZE 32 + #define TYPE float + #define STORAGE_ORDER Eigen::ColMajor + + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_col) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col) + #include "svd_model.cpp" + #undef CLASSNAME + #undef SOURCENAME + #undef STORAGE_ORDER + + #define STORAGE_ORDER Eigen::RowMajor + #define CLASSNAME ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row) + #define SOURCENAME ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row) + #include "svd_model.cpp" + #undef CLASSNAME + #undef SOURCENAME + #undef STORAGE_ORDER + + #undef SIZE + #undef TYPE +} diff --git a/src/svd_model.cpp b/src/svd_model.cpp new file mode 100644 index 0000000..deb3082 --- /dev/null +++ b/src/svd_model.cpp @@ -0,0 +1,111 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== + +//================================================================================================== +/* + This file is a X-File to generate various svd_impl_* definitions variant +*/ +//================================================================================================== + +//================================================================================================== +// Internal payload +//================================================================================================== +struct CLASSNAME::payload +{ + using matrix_type = Eigen::Matrix; + Eigen::JacobiSVD data; + + payload(SOURCENAME const& src, int options) : data(src.storage()->data,options) {} +}; + +//================================================================================================== +// Constructors & Special Members +//================================================================================================== +CLASSNAME::CLASSNAME(SOURCENAME const& m, int options) +: storage_(std::make_unique(m,options)) +{} + + +CLASSNAME::CLASSNAME(CLASSNAME const& o) +: storage_(std::make_unique(*o.storage_)) +{ +} + +CLASSNAME::CLASSNAME(CLASSNAME&&) noexcept = default; + +CLASSNAME& CLASSNAME::operator=(CLASSNAME const& o) +{ + if (this != &o) storage_->data = o.storage_->data; + return *this; +} + +CLASSNAME& CLASSNAME::operator=(CLASSNAME&& o) noexcept +{ + if (this != &o) storage_->data = std::move(o.storage_->data); + return *this; +} + +CLASSNAME::~CLASSNAME() = default; + +int CLASSNAME::rank() const { return storage_->data.rank(); } + +SOURCENAME CLASSNAME::U() const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.matrixU()); + return result; +} + +SOURCENAME CLASSNAME::U(int r) const +{ + SOURCENAME result; + result.storage()->assign( storage_->data.matrixU().leftCols(r) ); + return result; +} + +SOURCENAME CLASSNAME::singular_values() const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.singularValues()); + return result; +} + +SOURCENAME CLASSNAME::singular_values(int r) const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.singularValues().head(r)); + return result; +} + +SOURCENAME CLASSNAME::D() const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.singularValues().asDiagonal()); + return result; +} + +SOURCENAME CLASSNAME::D(int r) const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.singularValues().head(r).asDiagonal()); + return result; +} + +SOURCENAME CLASSNAME::V() const +{ + SOURCENAME result; + result.storage()->assign(storage_->data.matrixV()); + return result; +} + +SOURCENAME CLASSNAME::V(int r) const +{ + SOURCENAME result; + result.storage()->assign( storage_->data.matrixV().leftCols(r) ); + return result; +} diff --git a/test/unit/block/extract.cpp b/test/unit/block/extract.cpp index d820a75..dc1e7c5 100644 --- a/test/unit/block/extract.cpp +++ b/test/unit/block/extract.cpp @@ -30,7 +30,10 @@ MatrixType make_initialized_matrix(rotgen::tests::matrix_block_test_case() , ref.template lpNorm<1>()); - TTS_EQUAL(input.template lpNorm<2>() , ref.template lpNorm<2>()); - TTS_EQUAL(input.template lpNorm() , ref.template lpNorm()); - - TTS_EQUAL(norm(input) , ref.norm()); - TTS_EQUAL(squaredNorm(input) , ref.squaredNorm()); - TTS_EQUAL(lpNorm<1>(input) , ref.template lpNorm<1>()); - TTS_EQUAL(lpNorm<2>(input) , ref.template lpNorm<2>()); - TTS_EQUAL(lpNorm(input) , ref.template lpNorm()); + TTS_EQUAL(rotgen::norm(input) , ref.norm()); + TTS_EQUAL(rotgen::squaredNorm(input) , ref.squaredNorm()); + TTS_EQUAL(rotgen::lpNorm<1>(input) , ref.template lpNorm<1>()); + TTS_EQUAL(rotgen::lpNorm<2>(input) , ref.template lpNorm<2>()); + TTS_EQUAL(rotgen::lpNorm(input) , ref.template lpNorm()); if constexpr(T::IsVectorAtCompileTime) { @@ -40,14 +34,10 @@ namespace rotgen::tests mat_t norm_ref(input.rows(), input.cols()); prepare([&](auto r, auto c) { return e_norm(r,c); }, norm_ref); - TTS_EQUAL(input.normalized(), norm_ref); - TTS_EQUAL(normalized(input), norm_ref); - auto m_norm = input; - m_norm.normalize(); - TTS_EQUAL(m_norm, norm_ref); + TTS_EQUAL(rotgen::normalized(input), norm_ref); auto f_norm = input; - normalize(f_norm); + rotgen::normalize(f_norm); TTS_EQUAL(f_norm, norm_ref); } } diff --git a/test/unit/functions/svd.cpp b/test/unit/functions/svd.cpp new file mode 100644 index 0000000..6874817 --- /dev/null +++ b/test/unit/functions/svd.cpp @@ -0,0 +1,79 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#include "unit/tests.hpp" +#include + +TTS_CASE_TPL("SVD decomposition - Dynamic case", rotgen::tests::types) +( tts::type< tts::types> ) +{ + int rank, i = 5; + auto eps = std::numeric_limits::epsilon(); + + auto m = rotgen::matrix::Random(5,5); + auto decomp = rotgen::svd(m); + + do + { + rank = decomp.rank(); + + auto u = decomp.U(rank); + auto d = decomp.singular_values(rank); + auto dd = decomp.D(rank); + auto v = decomp.V(rank); + + TTS_EQUAL(rank, i); + + auto rec = (u * dd * rotgen::transpose(v)); + auto error = m - rec; + + TTS_LESS_EQUAL(rotgen::maxCoeff(rotgen::abs(error)) / eps, 16.) + << "Result:\n" << rec << "\n" + << "Residuals:\n" << error << "\n"; + + // Reduce rank by duplicating one column + i--; + col(m,i) = col(m,0); + decomp = rotgen::svd(m); + + }while(rank != 1); +}; + +TTS_CASE_TPL("SVD decomposition - Static case", rotgen::tests::types) +( tts::type< tts::types> ) +{ + int rank, i = 5; + auto eps = std::numeric_limits::epsilon(); + + auto m = rotgen::matrix::Random(); + auto decomp = rotgen::svd(m); + + do + { + rank = decomp.rank(); + + auto u = decomp.U(rank); + auto d = decomp.singular_values(rank); + auto dd = decomp.D(rank); + auto v = decomp.V(rank); + + TTS_EQUAL(rank, i); + + auto rec = (u * dd * rotgen::transpose(v)); + auto error = m - rec; + + TTS_LESS_EQUAL(rotgen::maxCoeff(rotgen::abs(error)) / eps, 16.) + << "Result:\n" << rec << "\n" + << "Residuals:\n" << error << "\n"; + + // Reduce rank by duplicating one column + i--; + col(m,i) = col(m,0); + decomp = rotgen::svd(m); + + }while(rank != 1); +}; \ No newline at end of file diff --git a/test/unit/matrix/constructors.cpp b/test/unit/matrix/constructors.cpp index 418b8bd..2e9db4f 100644 --- a/test/unit/matrix/constructors.cpp +++ b/test/unit/matrix/constructors.cpp @@ -35,13 +35,28 @@ TTS_CASE_TPL("Dynamic matrix constructor with row and columns", rotgen::tests::t TTS_EQUAL(matrix.cols(), rotgen::Index{5}); }; -TTS_CASE_TPL("Static matrix constructor with row and columns", rotgen::tests::types) -( tts::type< tts::types> ) +TTS_CASE_TPL("Static matrix constructor with row and columns", float, double) +( tts::type ) { - rotgen::matrix matrix(6, 11); + rotgen::matrix v2(6, 11); + rotgen::matrix w2(6, 11); - TTS_EQUAL(matrix.rows(), rotgen::Index{6}); - TTS_EQUAL(matrix.cols(), rotgen::Index{11}); + TTS_EQUAL(v2(0), T{6}); + TTS_EQUAL(v2(1), T{11}); + + TTS_EQUAL(w2(0), T{6}); + TTS_EQUAL(w2(1), T{11}); + + rotgen::matrix v3(6, 11, 125); + rotgen::matrix w3(6, 11, 125); + + TTS_EQUAL(v3(0), T{6}); + TTS_EQUAL(v3(1), T{11}); + TTS_EQUAL(v3(2), T{125}); + + TTS_EQUAL(w3(0), T{6}); + TTS_EQUAL(w3(1), T{11}); + TTS_EQUAL(w3(2), T{125}); }; TTS_CASE_TPL("Copy constructor produces identical but independent matrix", rotgen::tests::types) diff --git a/test/unit/matrix/generators.cpp b/test/unit/matrix/generators.cpp index 3b3a758..410b92d 100644 --- a/test/unit/matrix/generators.cpp +++ b/test/unit/matrix/generators.cpp @@ -51,7 +51,7 @@ TTS_CASE_TPL("Test zero", rotgen::tests::types) test_value(setZero(matrix{} ), 1, 1, 0); test_value(setZero(matrix{} ), 10, 10, 0); test_value(setZero(matrix{3, 4}), 3, 4, 0); - test_value(setZero(matrix{7, 5} ), 7, 5, 0); + test_value(setZero(matrix{} ), 7, 5, 0); test_value(setZero(matrix{9, 3} ), 9, 3, 0); test_value(setZero(matrix{2, 3} ), 2, 3, 0); }; @@ -72,7 +72,7 @@ TTS_CASE_TPL("Test ones", rotgen::tests::types) test_value(setOnes(matrix{} ), 1, 1, 1); test_value(setOnes(matrix{} ), 10, 10, 1); test_value(setOnes(matrix{3, 4}), 3, 4, 1); - test_value(setOnes(matrix{7, 5} ), 7, 5, 1); + test_value(setOnes(matrix{} ), 7, 5, 1); test_value(setOnes(matrix{9, 3} ), 9, 3, 1); test_value(setOnes(matrix{2, 3} ), 2, 3, 1); }; @@ -93,7 +93,7 @@ TTS_CASE_TPL("Test constant", rotgen::tests::types) test_value(setConstant(matrix{} , T(5.12)), 1, 1, T(5.12)); test_value(setConstant(matrix{} , T(5.12)), 10, 10, T(5.12)); test_value(setConstant(matrix{3, 4}, T(5.12)), 3, 4, T(5.12)); - test_value(setConstant(matrix{7, 5} , T(5.12)), 7, 5, T(5.12)); + test_value(setConstant(matrix{} , T(5.12)), 7, 5, T(5.12)); test_value(setConstant(matrix{9, 3} , T(5.12)), 9, 3, T(5.12)); test_value(setConstant(matrix{2, 3} , T(5.12)), 2, 3, T(5.12)); }; @@ -114,7 +114,7 @@ TTS_CASE_TPL("Test random", rotgen::tests::types) test_random(setRandom(matrix{} ), 1, 1); test_random(setRandom(matrix{} ), 10, 10); test_random(setRandom(matrix{3, 4}), 3, 4); - test_random(setRandom(matrix{7, 5} ), 7, 5); + test_random(setRandom(matrix{} ), 7, 5); test_random(setRandom(matrix{9, 3} ), 9, 3); test_random(setRandom(matrix{2, 3} ), 2, 3); }; @@ -135,7 +135,7 @@ TTS_CASE_TPL("Test identity", rotgen::tests::types) test_identity(setIdentity(matrix{} ), 1, 1); test_identity(setIdentity(matrix{} ), 10, 10); test_identity(setIdentity(matrix{3, 4}), 3, 4); - test_identity(setIdentity(matrix{7, 5} ), 7, 5); + test_identity(setIdentity(matrix{} ), 7, 5); test_identity(setIdentity(matrix{9, 3} ), 9, 3); test_identity(setIdentity(matrix{2, 3} ), 2, 3); }; \ No newline at end of file diff --git a/test/unit/matrix/inverse.cpp b/test/unit/matrix/inverse.cpp new file mode 100644 index 0000000..3850964 --- /dev/null +++ b/test/unit/matrix/inverse.cpp @@ -0,0 +1,62 @@ +//================================================================================================== +/* + ROTGEN - Runtime Overlay for Eigen + Copyright : CODE RECKONS + SPDX-License-Identifier: BSL-1.0 +*/ +//================================================================================================== +#include "unit/tests.hpp" +#include "unit/common/arithmetic.hpp" +#include + +TTS_CASE_TPL("Test dynamic matrix inverse", rotgen::tests::types) +( tts::type< tts::types> ) +{ + using mat_t = rotgen::matrix; + auto eps = std::numeric_limits::epsilon(); + + auto const cases = rotgen::tests::generate_matrix_references(); + for (const auto& [r, c, fn] : cases) + { + if(r == c) + { + auto input = mat_t::Random(r, c); + auto inv = rotgen::inverse(input); + + auto rec = input * inv; + auto id = mat_t::Identity(rotgen::rows(rec),rotgen::cols(rec)); + auto error = rec - id; + + TTS_LESS_EQUAL(rotgen::maxCoeff(rotgen::abs(error)) / eps, 64.) + << "Result:\n" << rec << "\n" + << "Residuals:\n" << error << "\n"; + } + } +}; + +TTS_CASE_TPL("Test static matrix inverse", rotgen::tests::types) +( tts::type< tts::types> ) +{ + using mat_t = rotgen::matrix; + auto eps = std::numeric_limits::epsilon(); + auto const cases = rotgen::tests::generate_static_matrix_references(); + + auto process = [&](D const&) + { + if constexpr(D::rows == D::cols) + { + auto input = rotgen::matrix::Random(); + auto inv = rotgen::inverse(input); + + auto rec = input * inv; + auto id = mat_t::Identity(rotgen::rows(rec),rotgen::cols(rec)); + auto error = rec - id; + + TTS_LESS_EQUAL(rotgen::maxCoeff(rotgen::abs(error)) / eps, 64.) + << "Result:\n" << rec << "\n" + << "Residuals:\n" << error << "\n"; + } + }; + + std::apply([&](auto const&... d) { (process(d),...);}, cases); +};