Implements reductions visitors broadcasting

Co-authored-by: kallore <kkaspar@codereckons.com>

See merge request oss/rotgen!6

include/rotgen/impl/matrix_impl64.hpp

@@ -52,6 +52,10 @@ namespace rotgen
     double minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
     double trace() const;
 
+    double squaredNorm() const;
+    double norm() const;
+    double lpNorm(int p) const;
+
     double&       operator()(std::size_t i, std::size_t j);
     double const& operator()(std::size_t i, std::size_t j) const;
 

include/rotgen/matrix.hpp

@@ -13,6 +13,8 @@
 
 namespace rotgen
 {
+  inline constexpr int Infinity = -1;
+
   template< typename Scalar, int Rows    = -1  , int Cols    = -1
           , int Options = 0, int MaxRows = Rows, int MaxCols = Cols
           >
@@ -170,6 +172,13 @@ namespace rotgen
       if constexpr(Cols != -1) assert(cols == Cols &&  "Mismatched between dynamic and static column size");
       return parent::Identity(rows, cols);
     }
+
+    template<int P>
+    double lpNorm() const {
+      assert(P == 1 || P == 2 || P == Infinity);
+      return parent::lpNorm(P);
+    }
+
   };
 
   template<typename S, int R, int C, int O, int MR, int MC>

src/matrix_impl64.cpp

@@ -122,6 +122,18 @@ namespace rotgen
 
   double matrix_impl64::trace() const { return storage_->data.trace(); }
 
+  double matrix_impl64::squaredNorm() const { return storage_->data.squaredNorm(); }
+  double matrix_impl64::norm() const { return storage_->data.norm(); }
+
+  double matrix_impl64::lpNorm(int p) const {
+    if (p == 1)
+      return storage_->data.lpNorm<1>();
+    else if (p == 2)
+      return storage_->data.lpNorm<2>();
+    else
+      return storage_->data.lpNorm<Eigen::Infinity>();
+  }
+
   //==================================================================================================
   // Operators
   //==================================================================================================

test/basic/norms.cpp

@@ -0,0 +1,61 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#define TTS_MAIN
+#include <rotgen/matrix.hpp>
+#include "tts.hpp"
+#include <Eigen/Dense>
+
+template <typename MatrixType>
+struct MatrixDescriptor
+{
+  std::size_t rows, cols;
+  std::function<void(MatrixType &, std::size_t, std::size_t)> init_fn;
+};
+
+template <typename MatrixType>
+void test_matrix_norms(std::size_t rows, std::size_t cols,
+                       const std::function<void(MatrixType &, std::size_t, std::size_t)> &init_fn)
+{
+  MatrixType matrix(rows, cols);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+    for (std::size_t c = 0; c < cols; ++c)
+    {
+      init_fn(matrix, r, c);
+      ref(r, c) = matrix(r, c);
+    }
+
+  constexpr double epsilon = 1e-10;
+
+  TTS_RELATIVE_EQUAL(matrix.norm(), ref.norm(), epsilon);
+  TTS_RELATIVE_EQUAL(matrix.squaredNorm(), ref.squaredNorm(), epsilon);
+
+  TTS_RELATIVE_EQUAL(matrix.template lpNorm<1>(), ref.lpNorm<1>(), epsilon);
+  TTS_RELATIVE_EQUAL(matrix.template lpNorm<2>(), ref.lpNorm<2>(), epsilon);
+  TTS_RELATIVE_EQUAL(matrix.template lpNorm<-1>(), ref.lpNorm<Eigen::Infinity>(), epsilon);
+}
+
+TTS_CASE("Matrix norm-related operations")
+{
+  std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
+      {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + c; }},
+      {2, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
+      {2, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -r*r*r - c*c - 1.23; }},
+      {17, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r*c + 0.98; }},
+      {1, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 42.0; }},
+      {10, 11, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = std::tan(r + c); }},
+      {1, 5, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -1.5*r + 2.56*c + 3.33; }},
+      {7, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
+      {0, 0, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
+  };
+
+  for (const auto& [rows, cols, init_fn] : test_matrices)
+    test_matrix_norms<rotgen::matrix<double>>(rows, cols, init_fn);
+
+};