Implements some missing functions

See merge request oss/rotgen!28

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()

include/rotgen/common/ref.hpp

@@ -17,13 +17,14 @@
 namespace rotgen
 {
   // Primary template: mutable ref
-  template<typename T, int Options = T::storage_order, typename Stride = stride>
+  template<typename T, int Options, typename Stride>
   class ref : private map<T, Options, Stride>
   {
     public:
     using parent          = map<T, Options, Stride>;
     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<typename Ref, int R, int C, bool I>
+    ref(block<Ref,R,C,I>&& 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<typename Ref, int R, int C, bool I>
     ref(block<Ref,R,C,I>& 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<typename A, int O, typename S, typename B, int P, typename T>
+  auto operator+=(ref<A,O,S> lhs, ref<B,P,T> rhs) -> decltype(lhs.base() += rhs.base())
+  {
+    return lhs.base() += rhs.base();
+  }
+
   template<typename A, int O, typename S, typename B, int P, typename T>
   auto operator-(ref<A,O,S> lhs, ref<B,P,T> rhs)
   {
     return lhs.base() - rhs.base();
   }
 
+  template<typename A, int O, typename S, typename B, int P, typename T>
+  auto operator-=(ref<A,O,S> lhs, ref<B,P,T> rhs) -> decltype(lhs.base() -= rhs.base())
+  {
+    return lhs.base() -= rhs.base();
+  }
+
   template<typename A, int O, typename S, typename B, int P, typename T>
   auto operator*(ref<A,O,S> lhs, ref<B,P,T> rhs)
   {
     return lhs.base() * rhs.base();
   }
 
+  template<typename A, int O, typename S, typename B, int P, typename T>
+  auto operator*=(ref<A,O,S> lhs, ref<B,P,T> rhs) -> decltype(lhs.base() *= rhs.base())
+  {
+    return lhs.base() *= rhs.base();
+  }
+
   template<typename A, int O, typename S>
   auto operator*(ref<A,O,S> lhs, std::convertible_to<typename A::value_type> auto s)
   {
@@ -297,6 +330,18 @@ namespace rotgen
     return lhs / s;
   }
 
+  template<typename A, int O, typename S>
+  auto inverse(ref<A,O,S> lhs) -> decltype(lhs.base().inverse())
+  {
+    return lhs.base().inverse();
+  }
+
+  template<typename A, int O, typename S, typename B, int P, typename T>
+  auto cross(ref<A,O,S> lhs, ref<B,P,T> 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<typename T> struct generalize;
@@ -314,10 +359,28 @@ namespace rotgen
   template<concepts::entity T> struct generalize<T>
   {
     static constexpr bool is_const = std::is_const_v<T>;
-    using base = matrix<typename T::value_type,Dynamic,Dynamic,T::storage_order>;
+    using base = matrix<typename T::value_type,T::RowsAtCompileTime,T::ColsAtCompileTime,T::storage_order>;
     using type = std::conditional_t<is_const,ref<base const>, ref<base>>;
   };
 
+  template<typename T, int O, typename S>
+  struct generalize<ref<T,O,S>>
+  {
+    using type = ref<T,O,S>;
+  };
+
+  template<typename T, int O, typename S>
+  struct generalize<ref<T,O,S> const>
+  {
+    using type = ref<T,O,S>;
+  };
+
+  template<concepts::entity T>
+  typename T::parent& base_of(T& a)
+  {
+    return a.base();
+  }
+
   template<concepts::entity T>
   typename T::parent const& base_of(T const& a)
   {
@@ -335,12 +398,18 @@ namespace rotgen
   {
     static constexpr bool is_const = std::is_const_v<T>;
     using concrete_type = decltype(std::declval<T>().eval());
-    using base = matrix<typename T::Scalar,Dynamic,Dynamic,concrete_type::Options&1>;
+    using base = matrix<typename T::Scalar,T::RowsAtCompileTime,T::ColsAtCompileTime,concrete_type::Options&1>;
     using type = std::conditional_t<is_const,ref<base const>, ref<base>>;
   };
 
   template<concepts::eigen_compatible T>
-  T const& base_of(T const& a)
+  auto const& base_of(T const& a)
+  {
+    return a;
+  }
+
+  template<concepts::eigen_compatible T>
+  auto& base_of(T& a)
   {
     return a;
   }

include/rotgen/concepts.hpp

@@ -9,6 +9,24 @@
 
 namespace rotgen::concepts
 {
+  //================================================================================================
+  //! @brief  Check if a type is a Rotgen reference.
+  //================================================================================================
+  template<typename T>
+  concept reference  = requires { typename std::remove_cvref_t<T>::rotgen_ref_tag; };
+
+  //================================================================================================
+  //! @brief  Check if a type is a vector.
+  //================================================================================================
+  template<typename T>
+  concept vector  =   T::IsVectorAtCompileTime;
+
+  //================================================================================================
+  //! @brief  Check if a type is a 3D vector.
+  //================================================================================================
+  template<typename T, int N>
+  concept vectorND  = vector<T> && (T::RowsAtCompileTime == N || T::ColsAtCompileTime == N);
+
   //================================================================================================
   //! @brief  Check if a type is a ROTGEN type.
   //================================================================================================
@@ -25,12 +43,12 @@ namespace rotgen::concepts
   template<typename T>
   concept eigen_compatible = requires(T const& a)
   {
-    typename T::Scalar;
+    typename std::remove_cvref_t<T>::Scalar;
 
-    { T::RowsAtCompileTime }    -> std::convertible_to<int>;
+    { std::remove_cvref_t<T>::RowsAtCompileTime }    -> std::convertible_to<int>;
-    { T::ColsAtCompileTime }    -> std::convertible_to<int>;
+    { std::remove_cvref_t<T>::ColsAtCompileTime }    -> std::convertible_to<int>;
-    { T::MaxRowsAtCompileTime } -> std::convertible_to<int>;
+    { std::remove_cvref_t<T>::MaxRowsAtCompileTime } -> std::convertible_to<int>;
-    { T::MaxColsAtCompileTime } -> std::convertible_to<int>;
+    { std::remove_cvref_t<T>::MaxColsAtCompileTime } -> std::convertible_to<int>;
 
     { a.eval() };
   };

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;

include/rotgen/dynamic/block.hpp

@@ -15,6 +15,9 @@
 
 namespace rotgen
 {
+  template<typename T, int Options = T::storage_order, typename Stride = stride>
+  class ref;
+
   template<typename Ref, int Rows = Dynamic, int Cols = Dynamic, bool Inner = false>
   class block : public find_block<Ref>
   {
@@ -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)

include/rotgen/dynamic/map.hpp

@@ -143,6 +143,26 @@ namespace rotgen
     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<typename R2, int O2, typename S2>
+    map& operator+=(map<R2,O2,S2> const& rhs) requires(!is_immutable)
     {
-      base() += static_cast<parent const&>(rhs);
+      base() += rhs.base();
       return *this;
     }
 
-    map& operator-=(map const& rhs) requires(!is_immutable)
+    template<typename R2, int O2, typename S2>
+    map& operator-=(map<R2,O2,S2> const& rhs) requires(!is_immutable)
     {
-      base() -= static_cast<parent const&>(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<typename R2, int O2, typename S2>
+    map& operator*=(map<R2,O2,S2> const& rhs) requires(!is_immutable)
     {
-      base() *= static_cast<parent const&>(rhs);
+      base() *= rhs.base();
       return *this;
     }
 
@@ -297,22 +320,26 @@ namespace rotgen
   template<typename R1, typename R2, int O1, typename S1, int O2, typename S2>
   typename map<R1,O1,S1>::concrete_type operator+(map<R1,O1,S1> const& lhs, map<R2,O2,S2> const& rhs)
   {
-    using concrete_type = typename map<R1,O1,S1>::concrete_type;
+    // REPLICATE TODO
-    return concrete_type(lhs.base().add(map<R1,O1,S1>(rhs)));
+    using map_type = map<R1 const,O1,S1>;
+    using concrete_type = typename map_type::concrete_type;
+    return concrete_type(map_type(lhs).base().add(map_type(rhs)));
   }
 
   template<typename R1, typename R2, int O1, typename S1, int O2, typename S2>
   typename map<R1,O1,S1>::concrete_type operator-(map<R1,O1,S1> const& lhs, map<R2,O2,S2> const& rhs)
   {
-    using concrete_type = typename map<R1,O1,S1>::concrete_type;
+    using map_type = map<R1 const,O1,S1>;
-    return concrete_type(lhs.base().sub(map<R1,O1,S1>(rhs)));
+    using concrete_type = typename map_type::concrete_type;
+    return concrete_type(map_type(lhs).base().sub(map_type(rhs)));
   }
 
   template<typename R1, typename R2, int O1, typename S1, int O2, typename S2>
   typename map<R1,O1,S1>::concrete_type operator*(map<R1,O1,S1> const& lhs, map<R2,O2,S2> const& rhs)
   {
-    using concrete_type = typename map<R1,O1,S1>::concrete_type;
+    using map_type = map<R1 const,O1,S1>;
-    return concrete_type(lhs.base().mul(map<R1,O1,S1>(rhs)));
+    using concrete_type = typename map_type::concrete_type;
+    return concrete_type(map_type(lhs).base().mul(map_type(rhs)));
   }
 
   template<typename R, int O, typename S>

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<Scalar>(vs)...})
+    {}
+
     matrix(parent const& base) : parent(base) {}
 
     matrix(std::initializer_list<std::initializer_list<Scalar>> init) : parent(init)
@@ -129,9 +136,9 @@ namespace rotgen
       return matrix(base().normalized());
     }
 
-    matrix transpose() const
+    matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,storage_order> transpose() const
     {
-      return matrix(base().transpose());
+      return matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,storage_order>(base().transpose());
     }
 
     matrix conjugate() const
@@ -139,9 +146,9 @@ namespace rotgen
       return matrix(base().conjugate());
     }
 
-    matrix adjoint() const
+    matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,storage_order> adjoint() const
     {
-      return matrix(base().adjoint());
+      return matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,storage_order>(base().adjoint());
     }
 
     void normalize()        requires(IsVectorAtCompileTime)

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 <rotgen/impl/svd.hpp>
+
+namespace rotgen
+{
+  template<typename M>
+  struct svd : find_svd<typename M::value_type,M::Options>
+  {
+    using parent = find_svd<typename M::value_type,M::Options>;
+
+    using m_type = matrix<typename M::value_type,Dynamic,Dynamic,M::Options>;
+    using d_type = matrix<typename M::value_type,M::ColsAtCompileTime,1>;
+
+    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);  }
+  };
+}

include/rotgen/extract.hpp

@@ -30,6 +30,12 @@ namespace rotgen
     detail::validate_extract(e,i0,j0,ni,nj);
     return block<Entity>(e, i0, j0, ni, nj);
   }
+  template<concepts::entity Entity>
+  auto extract(Entity const& e, Index i0, Index j0, Index ni, Index nj)
+  {
+    detail::validate_extract(e,i0,j0,ni,nj);
+    return block<Entity const>(e, i0, j0, ni, nj);
+  }
 
   template<Index NI, Index NJ, concepts::entity Entity>
   requires(NI!=-1 && NJ!=-1)

include/rotgen/fixed/block.hpp

@@ -12,6 +12,9 @@
 
 namespace rotgen
 {
+  template<typename T, int Options = T::storage_order, typename Stride = stride>
+  class ref;
+
   namespace detail
   {
     template<typename Ref, int Rows, int Cols, bool Inner, bool IsConst>
@@ -25,6 +28,12 @@ namespace rotgen
                                       >;
     };
 
+    template<typename T, int O, typename S, int Rows, int Cols, bool Inner, bool IsConst>
+    struct  compute_block_type<ref<T,O,S>,Rows,Cols,Inner,IsConst>
+          : compute_block_type<typename ref<T,O,S>::parent,Rows,Cols,Inner,IsConst>
+    {
+    };
+
     template<typename Ref, int Rows, int Cols, bool Inner, bool IsConst>
     using block_type = typename compute_block_type<Ref,Rows,Cols,Inner,IsConst>::type;
   }
@@ -57,7 +66,6 @@ namespace rotgen
     template<typename ET>
     using as_concrete_type = as_concrete_t<ET, matrix>;
 
-    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<Ref>)
     : 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<Ref>)
     : 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<Ref>)
     : 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<Ref>)
     : parent(r.base(),i0,j0,Rows,Cols)
     {}
 
+    // Constructs from rotgen::ref<T>
+    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<Ref>)
+    : 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<Ref>)
+    : 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<Ref>)
+    : 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<Ref>)
+    :  parent(r.base().base(),i0,j0,Rows,Cols)
+    {}
+
     template<typename B, Index R, Index C, bool I, int FS>
     block(block<B,R,C,I,FS> const& other) : parent(other.base())
     {}

include/rotgen/fixed/map.hpp

@@ -39,7 +39,6 @@ namespace rotgen
     using value_type    = typename std::remove_const_t<Ref>::value_type;
     using concrete_type = typename std::remove_const_t<Ref>::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<typename R2, int O2, typename S2>
+    map& operator+=(map<R2,O2,S2> const& rhs) requires(!is_immutable)
     {
       base() += rhs.base();
       return *this;
     }
 
-    map& operator-=(map const& rhs)
+    template<typename R2, int O2, typename S2>
+    map& operator-=(map<R2,O2,S2> const& rhs) requires(!is_immutable)
     {
       base() -= rhs.base();
       return *this;
     }
 
-    map& operator*=(map const& rhs)
+    template<typename R2, int O2, typename S2>
+    map& operator*=(map<R2,O2,S2> 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<decltype(base().inverse())>(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()));
     };

include/rotgen/fixed/matrix.hpp

@@ -38,9 +38,9 @@ namespace rotgen
     using concrete_type = matrix;
     using concrete_dynamic_type = matrix<value_type>;
 
-    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<Scalar>(vs)...})
+    {}
+
     matrix(concepts::entity auto const& other) : parent(other.base())
     {}
 

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 <Eigen/Dense>
+#include <Eigen/SVD>
+
+namespace rotgen
+{
+  template<typename M>
+  struct svd
+  {
+    using base    = Eigen::JacobiSVD<typename M::parent>;
+    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<d_type, matrix>{svd_.singularValues()};
+      else return svd_.singularValues();
+    }
+
+    auto U() const
+    {
+      if constexpr(!use_expression_templates) return as_concrete_t<u_type, matrix>{svd_.matrixU()};
+      else return svd_.matrixU();
+    }
+
+    auto V() const
+    {
+      if constexpr(!use_expression_templates) return as_concrete_t<v_type, matrix>{svd_.matrixV()};
+      else return svd_.matrixV();
+    }
+
+    auto D() const
+    {
+      auto d = svd_.singularValues().asDiagonal();
+      if constexpr(!use_expression_templates)
+        return as_concrete_t<decltype(d.toDenseMatrix()), matrix>{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<decltype(that), matrix>{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<decltype(that), matrix>{that};
+      else return that;
+    }
+
+    auto V(int r) const
+    {
+      auto that = svd_.matrixV().leftCols(r);
+      if constexpr(!use_expression_templates) return as_concrete_t<decltype(that), matrix>{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<decltype(d.toDenseMatrix()), matrix>{d.toDenseMatrix ()};
+      else
+        return d;
+    }
+
+    private:
+    base svd_;
+  };
+}

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 rows(auto const& m) requires(requires{ m.rows(); }){ return m.rows(); }
-  std::size_t cols(concepts::entity auto const& arg) { return arg.cols(); }
+  std::size_t cols(auto const& m) requires(requires{ m.cols(); }){ return m.cols(); }
-  std::size_t size(concepts::entity auto const& arg) { return arg.size(); }
+  std::size_t size(auto const& m) requires(requires{ m.size(); }){ return m.size(); }
 
-  template<concepts::entity E>
+  void resize(auto& a, int sz) requires requires{a.resize(sz);} { a.resize(sz); }
-  void resize(E& arg, int sz) requires requires{arg.resize(sz);}
-  {
-    arg.resize(sz);
-  }
 
-  template<concepts::entity E>
+  template<typename E>
-  void resize(E& arg, int new_rows, int new_cols)
+  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<concepts::entity E>
+  template<typename E>
-  void conservativeResize(E& arg, int sz) requires requires{arg.conservativeResize(sz);}
+  void conservativeResize(E& a, int sz) requires requires{a.conservativeResize(sz);}
   {
-    arg.conservativeResize(sz);
+    a.conservativeResize(sz);
   }
 
-  template<concepts::entity E>
+  template<typename E>
-  void conservativeResize(E& arg, int new_rows, int new_cols)
+  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) normalized (auto const& m) requires(requires{ m.normalized(); }) { return m.normalized();  }
-  decltype(auto) transpose (concepts::entity auto const& arg) { return arg.transpose();   }
+  decltype(auto) transpose  (auto const& m) requires(requires{ m.transpose();  }) { return m.transpose();   }
-  decltype(auto) conjugate (concepts::entity auto const& arg) { return arg.conjugate();   }
+  decltype(auto) conjugate  (auto const& m) requires(requires{ m.conjugate();  }) { return m.conjugate();   }
-  decltype(auto) adjoint   (concepts::entity auto const& arg) { return arg.adjoint();     }
+  decltype(auto) adjoint    (auto const& m) requires(requires{ m.adjoint();    }) { return m.adjoint();     }
 
-  void normalize(concepts::entity auto& arg)        { arg.normalize();        }
+  void normalize(auto& a)         requires(requires{ a.normalize(); })        { a.normalize(); }
-  void transposeInPlace(concepts::entity auto& arg) { arg.transposeInPlace(); }
+  void transposeInPlace(auto& a)  requires(requires{ a.transposeInPlace(); }) { a.transposeInPlace(); }
-  void adjointInPlace(concepts::entity auto& arg)   { arg.adjointInPlace();   }
+  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)
+  auto maxCoeff(auto const& arg) requires( requires{ arg.maxCoeff(); } ) { return arg.maxCoeff(); }
-  requires( requires{ arg.maxCoeff(); } )
+  auto minCoeff(auto const& arg) requires( requires{ arg.minCoeff(); } ) { return arg.minCoeff(); }
-  {
-    return arg.maxCoeff();
-  }
-
-  auto minCoeff(concepts::entity auto const& arg)     { return arg.minCoeff(); }
 
   template<std::integral IndexType>
   auto maxCoeff(concepts::entity auto const& arg, IndexType* row, IndexType* col)
@@ -166,8 +157,8 @@ namespace rotgen
     return arg.minCoeff(row, col);
   }
 
-  template<int P>
+  template<int P, typename T>
-  auto lpNorm(concepts::entity auto const& arg)
+  auto lpNorm(T const& arg) requires( requires{arg.template lpNorm<P>();} )
   {
     static_assert(P == 1 || P == 2 || P == Infinity);
     return arg.template lpNorm<P>();
@@ -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<typename A> auto inverse(A const& a)
+  {
+    if constexpr(!use_expression_templates)               return inverse(generalize_t<A const>(a));
+    else if constexpr(requires{typename A::rotgen_tag;})  return base_of(a).inverse();
+    else                                                  return a.inverse();
+  }
+
+  template<concepts::vectorND<3> L, concepts::vectorND<3> R> auto cross(L const& l, R const& r)
+  {
+    if constexpr(!use_expression_templates)
+      return cross(generalize_t<R const>(l),generalize_t<R const>(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);
+  }
 }

include/rotgen/impl/map_indirect.hpp

@@ -2,16 +2,20 @@
   #define TYPE double
 
   #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 <rotgen/impl/map_model.hpp>
   #undef CLASSNAME
   #undef SOURCENAME
+  #undef CLASSCONSTNAME
 
   #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 <rotgen/impl/map_model.hpp>
   #undef CLASSNAME
   #undef SOURCENAME
+  #undef CLASSCONSTNAME
 
   #undef SIZE
   #undef TYPE
@@ -20,16 +24,20 @@
   #define TYPE float
 
   #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 <rotgen/impl/map_model.hpp>
   #undef CLASSNAME
   #undef SOURCENAME
+  #undef CLASSCONSTNAME
 
   #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 <rotgen/impl/map_model.hpp>
   #undef CLASSNAME
   #undef SOURCENAME
+  #undef CLASSCONSTNAME
 
   #undef SIZE
   #undef TYPE

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;
 

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 <rotgen/detail/generators.hpp>
+#include <rotgen/config.hpp>
+#include <cstddef>
+#include <memory>
+
+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 <rotgen/impl/svd_model.hpp>
+  #undef CLASSNAME
+  #undef SOURCENAME
+
+  #define CLASSNAME     ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row)
+  #define SOURCENAME    ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include <rotgen/impl/svd_model.hpp>
+  #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 <rotgen/impl/svd_model.hpp>
+  #undef CLASSNAME
+  #undef SOURCENAME
+
+  #define CLASSNAME     ROTGEN_MATRIX_NAME(svd_impl,SIZE,_row)
+  #define SOURCENAME    ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include <rotgen/impl/svd_model.hpp>
+  #undef CLASSNAME
+  #undef SOURCENAME
+
+  #undef SIZE
+  #undef TYPE
+
+
+  template<typename Scalar,int Options> struct find_svd_impl;
+
+  template<> struct find_svd_impl<float , ColMajor> { using type = svd_impl32_col; };
+  template<> struct find_svd_impl<float , RowMajor> { using type = svd_impl32_row; };
+  template<> struct find_svd_impl<double, ColMajor> { using type = svd_impl64_col; };
+  template<> struct find_svd_impl<double, RowMajor> { using type = svd_impl64_row; };
+
+  template<typename Scalar,int Options>
+  using find_svd = typename find_svd_impl<Scalar,(Options & 1)>::type;
+}

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<payload> storage_;
+
+  public:
+  std::unique_ptr<payload>&       storage()       { return storage_; }
+  std::unique_ptr<payload> const& storage() const { return storage_; }
+};

include/rotgen/operators.hpp

@@ -68,4 +68,30 @@ namespace rotgen
     if constexpr(!use_expression_templates) return generalize_t<A const>(a) / b;
     else                                    return base_of(a) / b;
   }
+
+  //------------------------------------------------------------------------------------------------
+  // Compounds operators across types
+  template<typename A, typename B>
+  auto operator+=(A& a, B const& b)
+  requires(concepts::entity<A> || concepts::entity<B>)
+  {
+    if constexpr(!use_expression_templates) return generalize_t<A>(a) += generalize_t<B const>(b);
+    else                                    return base_of(a) += base_of(b);
+  }
+
+  template<typename A, typename B>
+  auto operator-=(A& a, B const& b)
+  requires(concepts::entity<A> || concepts::entity<B>)
+  {
+    if constexpr(!use_expression_templates) return generalize_t<A>(a) -= generalize_t<B const>(b);
+    else                                    return base_of(a) -= base_of(b);
+  }
+
+  template<typename A, typename B>
+  auto operator*=(A& a, B const& b)
+  requires(concepts::entity<A> || concepts::entity<B>)
+  {
+    if constexpr(!use_expression_templates) return generalize_t<A>(a) *= generalize_t<B const>(b);
+    else                                    return base_of(a) *= base_of(b);
+  }
 }

include/rotgen/rotgen.hpp

@@ -16,10 +16,12 @@
 #include <rotgen/dynamic/matrix.hpp>
 #include <rotgen/dynamic/block.hpp>
 #include <rotgen/dynamic/map.hpp>
+#include <rotgen/dynamic/svd.hpp>
 #else
 #include <rotgen/fixed/matrix.hpp>
 #include <rotgen/fixed/block.hpp>
 #include <rotgen/fixed/map.hpp>
+#include <rotgen/fixed/svd.hpp>
 #endif
 
 #include <rotgen/common/ref.hpp>

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
 
 

src/map.cpp

@@ -9,6 +9,7 @@
 #include <rotgen/impl/map.hpp>
 #include <rotgen/impl/payload.hpp>
 #include <Eigen/Dense>
+#include <Eigen/LU>
 
 namespace rotgen
 {

src/map_indirect.cpp

@@ -3,18 +3,22 @@
   #define STORAGE_ORDER Eigen::ColMajor
 
   #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 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
@@ -25,17 +29,21 @@
   #define STORAGE_ORDER Eigen::ColMajor
 
   #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 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
 

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;

src/svd.cpp

@@ -0,0 +1,61 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include <rotgen/detail/generators.hpp>
+#include <rotgen/impl/matrix.hpp>
+#include <rotgen/impl/payload.hpp>
+#include <rotgen/impl/svd.hpp>
+#include <rotgen/config.hpp>
+#include <Eigen/Dense>
+#include <Eigen/SVD>
+
+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
+}

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<TYPE,Eigen::Dynamic,Eigen::Dynamic,STORAGE_ORDER>;
+  Eigen::JacobiSVD<matrix_type> 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<payload>(m,options))
+{}
+
+
+CLASSNAME::CLASSNAME(CLASSNAME const& o)
+: storage_(std::make_unique<payload>(*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;
+}

test/unit/block/extract.cpp

@@ -30,7 +30,10 @@ MatrixType make_initialized_matrix(rotgen::tests::matrix_block_test_case<MatrixT
   auto[d,i0,j0,ni,nj] = matrix_construct;
   auto[r,c,fn] = d;
 
-  MatrixType matrix(r, c);
+  MatrixType matrix;
+
+  if constexpr(MatrixType::RowsAtCompileTime == -1 &&MatrixType::ColsAtCompileTime == -1)
+    rotgen::resize(matrix, r, c);
 
   for(rotgen::Index i = 0; i < r; ++i)
     for(rotgen::Index j = 0; j < c; ++j)

test/unit/common/norms.hpp

@@ -21,17 +21,11 @@ namespace rotgen::tests
 
     prepare([&](auto r, auto c) { return input(r,c); }, ref);
 
-    TTS_EQUAL(input.norm()                              , ref.norm());
+    TTS_EQUAL(rotgen::norm(input)                     , ref.norm());
-    TTS_EQUAL(input.squaredNorm()                       , ref.squaredNorm());
+    TTS_EQUAL(rotgen::squaredNorm(input)              , ref.squaredNorm());
-    TTS_EQUAL(input.template lpNorm<1>()                , ref.template lpNorm<1>());
+    TTS_EQUAL(rotgen::lpNorm<1>(input)                , ref.template lpNorm<1>());
-    TTS_EQUAL(input.template lpNorm<2>()                , ref.template lpNorm<2>());
+    TTS_EQUAL(rotgen::lpNorm<2>(input)                , ref.template lpNorm<2>());
-    TTS_EQUAL(input.template lpNorm<rotgen::Infinity>() , ref.template lpNorm<Eigen::Infinity>());
+    TTS_EQUAL(rotgen::lpNorm<rotgen::Infinity>(input) , ref.template lpNorm<Eigen::Infinity>());
-
-    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<rotgen::Infinity>(input) , ref.template lpNorm<Eigen::Infinity>());
 
     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(rotgen::normalized(input), norm_ref);
-      TTS_EQUAL(normalized(input), norm_ref);
-      auto m_norm = input;
-      m_norm.normalize();
-      TTS_EQUAL(m_norm, norm_ref);
 
       auto f_norm = input;
-      normalize(f_norm);
+      rotgen::normalize(f_norm);
       TTS_EQUAL(f_norm, norm_ref);
     }
   }

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 <rotgen/rotgen.hpp>
+
+TTS_CASE_TPL("SVD decomposition - Dynamic case", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  int rank, i = 5;
+  auto eps = std::numeric_limits<T>::epsilon();
+
+  auto m = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>::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)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  int rank, i = 5;
+  auto eps = std::numeric_limits<T>::epsilon();
+
+  auto m = rotgen::matrix<T,5,5,O::value>::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);
+};

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_CASE_TPL("Static matrix constructor with row and columns", float, double)
-<typename T, typename O>( tts::type< tts::types<T,O>> )
+<typename T>( tts::type<T> )
 {
-  rotgen::matrix<T,6,11,O::value> matrix(6, 11);
+  rotgen::matrix<T,1,2> v2(6, 11);
+  rotgen::matrix<T,2,1> w2(6, 11);
 
-  TTS_EQUAL(matrix.rows(), rotgen::Index{6});
+  TTS_EQUAL(v2(0), T{6});
-  TTS_EQUAL(matrix.cols(), rotgen::Index{11});
+  TTS_EQUAL(v2(1), T{11});
+
+  TTS_EQUAL(w2(0), T{6});
+  TTS_EQUAL(w2(1), T{11});
+
+  rotgen::matrix<T,1,3> v3(6, 11, 125);
+  rotgen::matrix<T,3,1> 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)

test/unit/matrix/generators.cpp

@@ -51,7 +51,7 @@ TTS_CASE_TPL("Test zero", rotgen::tests::types)
   test_value(setZero(matrix<T, 1, 1, O::value>{}              ),  1,  1, 0);
   test_value(setZero(matrix<T, 10, 10, O::value>{}            ), 10, 10, 0);
   test_value(setZero(matrix<T,Dynamic,Dynamic, O::value>{3, 4}),  3,  4, 0);
-  test_value(setZero(matrix<T, 7, 5, O::value>{7, 5}          ),  7,  5, 0);
+  test_value(setZero(matrix<T, 7, 5, O::value>{}              ),  7,  5, 0);
   test_value(setZero(matrix<T, 9,Dynamic>{9, 3}               ),  9,  3, 0);
   test_value(setZero(matrix<T,Dynamic, 3>{2, 3}               ),  2,  3, 0);
 };
@@ -72,7 +72,7 @@ TTS_CASE_TPL("Test ones", rotgen::tests::types)
   test_value(setOnes(matrix<T, 1, 1, O::value>{}              ),  1,  1, 1);
   test_value(setOnes(matrix<T, 10, 10, O::value>{}            ), 10, 10, 1);
   test_value(setOnes(matrix<T,Dynamic,Dynamic, O::value>{3, 4}),  3,  4, 1);
-  test_value(setOnes(matrix<T, 7, 5, O::value>{7, 5}          ),  7,  5, 1);
+  test_value(setOnes(matrix<T, 7, 5, O::value>{}              ),  7,  5, 1);
   test_value(setOnes(matrix<T, 9,Dynamic>{9, 3}               ),  9,  3, 1);
   test_value(setOnes(matrix<T,Dynamic, 3>{2, 3}               ),  2,  3, 1);
 };
@@ -93,7 +93,7 @@ TTS_CASE_TPL("Test constant", rotgen::tests::types)
   test_value(setConstant(matrix<T, 1, 1, O::value>{}              , T(5.12)),  1,  1, T(5.12));
   test_value(setConstant(matrix<T, 10, 10, O::value>{}            , T(5.12)), 10, 10, T(5.12));
   test_value(setConstant(matrix<T,Dynamic,Dynamic, O::value>{3, 4}, T(5.12)),  3,  4, T(5.12));
-  test_value(setConstant(matrix<T, 7, 5, O::value>{7, 5}          , T(5.12)),  7,  5, T(5.12));
+  test_value(setConstant(matrix<T, 7, 5, O::value>{}              , T(5.12)),  7,  5, T(5.12));
   test_value(setConstant(matrix<T, 9,Dynamic>{9, 3}               , T(5.12)),  9,  3, T(5.12));
   test_value(setConstant(matrix<T,Dynamic, 3>{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<T, 1, 1, O::value>{}              ),  1,  1);
   test_random(setRandom(matrix<T, 10, 10, O::value>{}            ), 10, 10);
   test_random(setRandom(matrix<T,Dynamic,Dynamic, O::value>{3, 4}),  3,  4);
-  test_random(setRandom(matrix<T, 7, 5, O::value>{7, 5}          ),  7,  5);
+  test_random(setRandom(matrix<T, 7, 5, O::value>{}              ),  7,  5);
   test_random(setRandom(matrix<T, 9,Dynamic>{9, 3}               ),  9,  3);
   test_random(setRandom(matrix<T,Dynamic, 3>{2, 3}               ),  2,  3);
 };
@@ -135,7 +135,7 @@ TTS_CASE_TPL("Test identity", rotgen::tests::types)
   test_identity(setIdentity(matrix<T, 1, 1, O::value>{}              ),  1,  1);
   test_identity(setIdentity(matrix<T, 10, 10, O::value>{}            ), 10, 10);
   test_identity(setIdentity(matrix<T,Dynamic,Dynamic, O::value>{3, 4}),  3,  4);
-  test_identity(setIdentity(matrix<T, 7, 5, O::value>{7, 5}          ),  7,  5);
+  test_identity(setIdentity(matrix<T, 7, 5, O::value>{}              ),  7,  5);
   test_identity(setIdentity(matrix<T, 9,Dynamic>{9, 3}               ),  9,  3);
   test_identity(setIdentity(matrix<T,Dynamic, 3>{2, 3}               ),  2,  3);
 };

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 <rotgen/rotgen.hpp>
+
+TTS_CASE_TPL("Test dynamic matrix inverse", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto eps = std::numeric_limits<T>::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)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto eps = std::numeric_limits<T>::epsilon();
+  auto const cases = rotgen::tests::generate_static_matrix_references();
+
+  auto process = [&]<typename D>(D const&)
+  {
+    if constexpr(D::rows == D::cols)
+    {
+      auto input = rotgen::matrix<T,D::rows,D::cols,O::value>::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);
+};