Feat/non member functions

Co-authored-by: Karen <kkaspar@codereckons.com>
Co-authored-by: Joel FALCOU <jfalcou@codereckons.com>

See merge request oss/rotgen!7

include/rotgen/functions.hpp

@@ -0,0 +1,142 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+#include <rotgen/matrix.hpp>
+
+namespace rotgen
+{
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  std::size_t rows(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.rows();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  std::size_t cols(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.cols();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  std::size_t size(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.size();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  void resize(rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix, int new_rows, int new_cols)
+  requires(Rows == -1 && Cols == -1)
+  {
+    matrix.resize(new_rows, new_cols);
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  void conservativeResize(rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix, int new_rows, int new_cols)
+  requires(Rows == -1 && Cols == -1)
+  {
+    matrix.conservativeResize(new_rows, new_cols);
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  rotgen::matrix<Scalar, Cols, Rows, Options, MaxCols, MaxRows> transpose(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.transpose();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> conjugate(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.conjugate();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  rotgen::matrix<Scalar, Cols, Rows, Options, MaxCols, MaxRows> adjoint(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.adjoint();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  void transposeInPlace(rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    matrix.transposeInPlace();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  void adjointInPlace(rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    matrix.adjointInPlace();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double squaredNorm(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.squaredNorm();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double norm(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.norm();
+  }
+
+  template< int P, typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double lpNorm(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    static_assert(P == 1 || P == 2 || P == Infinity);
+    return matrix.template lp_norm<P>();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double sum(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.sum();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double prod(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.prod();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double mean(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.mean();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double maxCoeff(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.maxCoeff();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double maxCoeff(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix, std::ptrdiff_t* row, std::ptrdiff_t* col)
+  {
+    return matrix.maxCoeff(row, col);
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double minCoeff(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.minCoeff();
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double minCoeff(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix, std::ptrdiff_t* row, std::ptrdiff_t* col)
+  {
+    return matrix.minCoeff(row, col);
+  }
+
+  template< typename Scalar, int Rows = -1, int Cols = 1, int Options = 0, int MaxRows = Rows, int MaxCols = Cols >
+  double trace(const rotgen::matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> &matrix)
+  {
+    return matrix.trace();
+  }
+
+}

include/rotgen/impl/matrix_model.hpp

@@ -52,7 +52,7 @@ class MATRIX
 
   TYPE squaredNorm()  const;
   TYPE norm()         const;
-  TYPE lpNorm(int p)  const;
+  TYPE lp_norm(int p) const;
 
   TYPE&       operator()(std::size_t i, std::size_t j);
   TYPE const& operator()(std::size_t i, std::size_t j) const;

include/rotgen/matrix.hpp

@@ -169,10 +169,10 @@ namespace rotgen
     }
 
     template<int P>
-    Scalar lpNorm() const
+    Scalar lp_norm() const
     {
-      assert(P == 1 || P == 2 || P == Infinity);
-      return parent::lpNorm(P);
+      static_assert(P == 1 || P == 2 || P == Infinity);
+      return parent::lp_norm(P);
     }
 
   };

src/matrix_model.cpp

@@ -132,7 +132,7 @@ TYPE MATRIX::maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return s
 TYPE MATRIX::squaredNorm()  const { return storage_->data.squaredNorm(); }
 TYPE MATRIX::norm()         const { return storage_->data.norm(); }
 
-TYPE MATRIX::lpNorm(int p) const
+TYPE MATRIX::lp_norm(int p) const
 {
   if (p == 1)       return storage_->data.lpNorm<1>();
   else if (p == 2)  return storage_->data.lpNorm<2>();

test/CMakeLists.txt

@@ -18,3 +18,4 @@ target_link_libraries(rotgen_test INTERFACE rotgen Eigen3::Eigen)
 ## Unit Tests registration
 ##======================================================================================================================
 rotgen_glob_unit(PATTERN "unit/basic/*.cpp"  INTERFACE rotgen_test)
+rotgen_glob_unit(PATTERN "unit/free_functions/*.cpp"  INTERFACE rotgen_test)

test/unit/basic/norms.cpp

@@ -40,8 +40,8 @@ TTS_CASE_TPL("Matrix norm-related operations", rotgen::tests::types)
 
     TTS_EQUAL(matrix.norm(), ref.norm());
     TTS_EQUAL(matrix.squaredNorm(), ref.squaredNorm());
-    TTS_EQUAL(matrix.template lpNorm<1>()               , ref.template lpNorm<1>());
-    TTS_EQUAL(matrix.template lpNorm<2>()               , ref.template lpNorm<2>());
-    TTS_EQUAL(matrix.template lpNorm<rotgen::Infinity>(), ref.template lpNorm<Eigen::Infinity>());
+    TTS_EQUAL(matrix.template lp_norm<1>()               , ref.template lpNorm<1>());
+    TTS_EQUAL(matrix.template lp_norm<2>()               , ref.template lpNorm<2>());
+    TTS_EQUAL(matrix.template lp_norm<rotgen::Infinity>(), ref.template lpNorm<Eigen::Infinity>());
   }
 };

test/unit/free_functions/arithmetic_functions.cpp

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

test/unit/free_functions/functions.cpp

@@ -0,0 +1,92 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#define TTS_MAIN
+#include <rotgen/matrix.hpp>
+#include <rotgen/functions.hpp>
+#include "tts.hpp"
+#include <functional>
+#include <array>
+
+
+
+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, std::size_t N>
+void test_matrix_sizes(std::size_t rows, std::size_t cols,
+                       const std::function<void(MatrixType &, std::size_t, std::size_t)> &init_fn,
+                       const std::array<std::pair<int, int>, N>& resize_dimensions)
+{
+  MatrixType matrix(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);
+
+  TTS_EQUAL(rotgen::size(matrix), rows*cols);
+
+  for (std::size_t i = 0; i < N; ++i) {
+    int r = resize_dimensions[i].first;
+    int c = resize_dimensions[i].second;
+
+    rotgen::resize(matrix, r, c);
+
+    TTS_EQUAL(rotgen::rows(matrix), static_cast<size_t>(r));
+    TTS_EQUAL(rotgen::cols(matrix), static_cast<size_t>(c));
+    TTS_EQUAL(rotgen::size(matrix), static_cast<size_t>(r*c));
+  }
+
+  rotgen::conservativeResize(matrix, rows, cols);
+  TTS_EQUAL(rotgen::size(matrix), rows*cols);
+
+  int i = 1;
+  for(std::size_t r=0;r<rows;++r)
+    for(std::size_t c=0;c<cols;++c)
+      matrix(r, c) = i++;
+
+  rotgen::conservativeResize(matrix, rows + 3, cols + 2);
+  TTS_EQUAL(rotgen::size(matrix), (rows + 3)*(cols + 2));
+
+  i = 1;
+  for(std::size_t r=0;r<rows;++r)
+    for(std::size_t c=0;c<cols;++c)
+      TTS_EQUAL(matrix(r,c),i++);
+}
+
+TTS_CASE("Matrix size-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_sizes<rotgen::matrix<double>>(rows, cols, init_fn, std::array{
+      std::pair{1, 2}, std::pair{11, 17}, std::pair{4, 5},
+      std::pair{9, 1}, std::pair{3, 8}, std::pair{22, 0},
+  });
+
+  rotgen::matrix<double, 3, 3> a;
+  TTS_EXPECT_NOT_COMPILES(a, { rotgen::resize(a, 4, 5); });
+  TTS_EXPECT_NOT_COMPILES(a, { rotgen::conservativeResize(a, 4, 5); });
+};
+
+
+
+

test/unit/free_functions/norms.cpp

@@ -0,0 +1,62 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#define TTS_MAIN
+#include <rotgen/matrix.hpp>
+#include <rotgen/functions.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(rotgen::norm(matrix), ref.norm(), epsilon);
+  TTS_RELATIVE_EQUAL(rotgen::squaredNorm(matrix), ref.squaredNorm(), epsilon);
+
+  TTS_RELATIVE_EQUAL(rotgen::lpNorm<1>(matrix), ref.lpNorm<1>(), epsilon);
+  TTS_RELATIVE_EQUAL(rotgen::lpNorm<2>(matrix), ref.lpNorm<2>(), epsilon);
+  TTS_RELATIVE_EQUAL(rotgen::lpNorm<-1>(matrix), 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);
+
+};