Create compatibility between vectors and matrices

include/cotila/cotila.h

@@ -7,7 +7,6 @@
 #include <cotila/matrix/utility.h>
 #include <cotila/scalar/math.h>
 #include <cotila/vector/math.h>
-#include <cotila/vector/operators.h>
 #include <cotila/vector/utility.h>
 #include <cotila/vector/vector.h>
 

include/cotila/detail/tmp.h

@@ -19,6 +19,15 @@ template <typename T, T v> struct all_same_value<T, v> {
   static constexpr T value = v;
 };
 
+template<typename T>
+static auto test_callable(int)
+        -> decltype(void(std::declval<T>()(std::declval<std::size_t>())), std::true_type{});
+template<typename T>
+static auto test_callable(long) -> std::false_type;
+
+template<typename T>
+struct is_callable_with_a_size_t : decltype(detail::test_callable<T>(0)){};
+
 } // namespace detail
 } // namespace cotila
 

include/cotila/matrix/math.h

@@ -7,10 +7,9 @@
 #include<algorithm>
 
 #include <cotila/scalar/math.h>
-#include <cotila/vector/vector.h>
-#include <cotila/vector/math.h>
 #include <cotila/matrix/matrix.h>
 #include <cotila/matrix/utility.h>
+#include <cotila/vector/math.h>
 #include <cotila/detail/assert.h>
 
 namespace cotila {
@@ -31,7 +30,7 @@ constexpr matrix<T, M, N> conj(const matrix<T, M, N> &m) {
   return elementwise(cotila::conj<T>, m);
 }
 
-/** @brief computes the elementwise real 
+/** @brief computes the elementwise real
  *  @param m an \f$ M \times N \f$ matrix of type T
  *  @return an \f$ M \times N \f$ matrix \f$\textbf{m}\f$ of type T such that
  *  \f$ \left(\textbf{m}}_{ij} = \mathbb{R}\{\textbf{m}_{ij}\},\ \forall i,j \f$
@@ -44,7 +43,7 @@ real(const matrix<T, M, N> &m) {
   return elementwise([](auto i) { return std::real(i); }, m);
 }
 
-/** @brief computes the elementwise imag 
+/** @brief computes the elementwise imag
  *  @param m an \f$ M \times N \f$ matrix of type T
  *  @return an \f$ M \times N \f$ matrix \f$\textbf{m}\f$ of type T such that
  *  \f$ \left(\textbf{m}}_{ij} = \mathbb{I}\{\textbf{m}_{ij}\},\ \forall i,j \f$
@@ -57,6 +56,47 @@ imag(const matrix<T, M, N> &m) {
   return elementwise([](auto i) { return std::imag(i); }, m);
 }
 
+/** @brief computes the elementwise absolute value
+ *  @param v an N-vector of type T
+ *  @return an N-vector \f$ \begin{bmatrix} \lvert v_1 \rvert & \ldots & \lvert v_N \rvert \end{bmatrix} \f$ of type T
+ *
+ *  Computes the elementwise absolute value of a matrix.
+ */
+template <typename T, std::size_t N, std::size_t M>
+constexpr matrix<detail::remove_complex_t<T>, N, M> abs(const matrix<T, N, M> &v) {
+  return elementwise(abs<T>, v);
+}
+
+/** @brief computes the sum of elements of a matrix
+ *  @param m an \f$ N \times M \f$ matrix of type T
+ *  @return a scalar \f$ \sum\limits_{i} m_i \f$ of type T
+ *
+ *  Computes the sum of the elements of a matrix.
+ */
+template <typename T, std::size_t N, std::size_t M> constexpr T sum(const matrix<T, N, M> &v) {
+  return accumulate(v, T{}, std::plus<T>());
+}
+
+/** @brief computes the minimum valued element
+ *  @param m an \f$ N \times M \f$ matrix of type T
+ *  @return a scalar \f$ v_i \f$ of type T where \f$ v_i \leq v_j,\ \forall j \f$
+ *
+ *  Computes the minimum valued element of a matrix.
+ */
+template <typename T, std::size_t N, std::size_t M> constexpr T min(const matrix<T, N, M> &m) {
+  return accumulate(m, m[0], [](T a, T b) { return std::min(a, b); });
+}
+
+/** @brief computes the maximum valued element
+ *  @param m an \f$ N \times M \f$ matrix of type T
+ *  @return a scalar \f$ v_i \f$ of type T where \f$ v_i \geq v_j,\ \forall j \f$
+ *
+ *  Computes the maximum valued element of a matrix.
+ */
+template <typename T, std::size_t N, std::size_t M> constexpr T max(const matrix<T, N, M> &m) {
+  return accumulate(m, m[0], [](T a, T b) { return std::max(a, b); });
+}
+
 /** @brief computes the transpose
  *  @param m an \f$ M \times N \f$ matrix of type T
  *  @return an \f$ N \times M \f$ matrix \f$ \textbf{m}^{\mathrm{T}} \f$ of type T such that
@@ -269,7 +309,7 @@ constexpr T det(const matrix<T, M, M> &m) {
 template <typename T, std::size_t M>
 constexpr matrix<T, M, M> inverse(const matrix<T, M, M> &m) {
   if (rank(m) < M)
-    throw "matrix is not invertible";
+    throw std::logic_error("matrix is not invertible");
   return submat<M, M>(rref(horzcat(m, identity<T, M>)), 0, M);
 }
 
@@ -285,6 +325,16 @@ constexpr T trace(const matrix<T, M, M> &m) {
   return sum(generate<M>([&m](std::size_t i){ return m[i][i]; }));
 }
 
+/** @brief computes the elementwise square root
+ *  @param m an \f$ N \times M \f$ matrix of type T
+ *  @return an \f$ M \times N \f$ matrix of type T, the element wise
+ *  square root of \f$ \textbf{m} \f$
+ *  Computes the elementwise square root of a matrix.
+ */
+template <typename T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sqrt(const matrix<T, N, M> &v) {
+  return elementwise(static_cast<T (*)(T)>(sqrt), v);
+}
 /** }@*/
 
 } // namespace cotila

include/cotila/matrix/matrix.h

@@ -6,7 +6,7 @@
 #define COTILA_MATRIX_MATRIX_H_
 
 #include <array>
-#include <cotila/vector/utility.h>
+#include <cotila/matrix/utility.h>
 #include <cotila/vector/vector.h>
 #include <cotila/detail/assert.h>
 #include <tuple>
@@ -31,6 +31,7 @@ template <typename T, std::size_t N, std::size_t M> struct matrix {
   using size_type = std::size_t;
   static constexpr size_type column_size = N; ///< Number of rows
   static constexpr size_type row_size = M;    ///< Number of columns
+  static constexpr size_type size = N * M;    ///< Number of elements
 
   /** @name Element access */
   ///@{
@@ -42,7 +43,7 @@ template <typename T, std::size_t N, std::size_t M> struct matrix {
    */
   constexpr vector<T, M> row(std::size_t i) const {
     if (i >= N)
-      throw "index out of range";
+      throw std::out_of_range("Row index is out-of-range");
     return generate<M>([i, this](std::size_t j) { return arrays[i][j]; });
   }
 
@@ -54,10 +55,37 @@ template <typename T, std::size_t N, std::size_t M> struct matrix {
    */
   constexpr vector<T, N> column(std::size_t i) const {
     if (i >= M)
-      throw "index out of range";
+      throw std::out_of_range("Column index is out-of-range");
     return generate<N>([i, this](std::size_t j) { return arrays[j][i]; });
   }
 
+  /**
+   * @brief acess specified element
+   * @param row selected line
+   * @param column selected row
+   * @return the selected scalar element
+   *
+   * Allow an access independent to the operator[],
+   * ease the writing of generic algorithm
+   */
+  constexpr T &at(std::size_t row, std::size_t column = 0) {
+    return at_impl(*this, row, column);
+  }
+
+  /**
+   * @copydoc at
+   */
+  constexpr const T &at(std::size_t row, std::size_t column = 0) const {
+    return at_impl(*this, row, column);
+  }
+
+  template <typename Self>
+  static constexpr auto &at_impl(Self &self, std::size_t row, std::size_t column = 0) {
+    if (row >= N || column >= M)
+      throw std::out_of_range("Indexes out of range");
+    return self.arrays[row][column];
+  } ///< @private
+
   /** @brief access specified element
    *  @param i index of the row
    *  @return pointer to the specified row
@@ -67,19 +95,82 @@ template <typename T, std::size_t N, std::size_t M> struct matrix {
    *  row pointer.  For a matrix `m`, accessing the element in the 5th row
    *  and 3rd column can be done with `m[5][3]`.
    */
-  constexpr T *operator[](std::size_t i) { return arrays[i]; }
+  template <std::size_t L = M>
+  typename std::enable_if<L == 1, T&>::type
+  constexpr operator[](std::size_t i) { return arrays[i][0]; }
+
+  /// @copydoc operator[]
+  template <std::size_t L = M>
+  typename std::enable_if<L == 1, const T&>::type
+  constexpr operator[](std::size_t i) const { return arrays[i][0]; }
 
   /// @copydoc operator[]
-  constexpr T const *operator[](std::size_t i) const { return arrays[i]; }
+  template <std::size_t L = M>
+  typename std::enable_if<L != 1, T*>::type
+  constexpr operator[](std::size_t i) { return arrays[i]; }
+
+  /// @copydoc operator[]
+  template <std::size_t L = M>
+  typename std::enable_if<L != 1, const T*>::type
+  constexpr operator[](std::size_t i) const { return arrays[i]; }
+  ///@}
+
+  /** @name Iterators */
+  ///@{
+  /** @brief returns an iterator to the beginning
+   *  @return an iterator to the beginning
+   *
+   *  Returns an iterator to the beginning of the matrix.
+   */
+  constexpr T *begin() noexcept { return *arrays; }
+
+  /** @brief returns an iterator to the end
+   *  @return an iterator past the end
+   *
+   *  Returns an iterator to the end of the matrix.
+   */
+  constexpr T *end() noexcept { return *(arrays + size); }
+
+  /// @copydoc begin
+  constexpr const T *cbegin() const noexcept { return *arrays; }
+  /// @copydoc begin
+  constexpr const T *begin() const noexcept { return *arrays; }
+
+  /// @copydoc end
+  constexpr const T *cend() const noexcept { return *(arrays + size); }
+  /// @copydoc end
+  constexpr const T *end() const noexcept { return *(arrays + size); }
   ///@}
 
+  constexpr operator vector<T, N>() const {
+    vector<T, N> ret{};
+    for (size_type i = 0; i < column_size; ++i)
+            ret[i] = arrays[i][0];
+    return ret;
+  }
+
+  template <typename U, std::size_t O, typename Container>
+  friend constexpr matrix<U, O, 1> from_initializer(Container &);
+
   T arrays[N][M]; ///< @private
 };
 
 /** \addtogroup matrix
  *  @{
  */
 
+template <typename T, std::size_t N, typename Container>
+constexpr matrix<T, N, 1> from_initializer(Container & init) {
+  if(std::size(init) != N)
+    throw std::invalid_argument("The initializer has not the same size than the vector");
+
+  matrix<T, N, 1> ret{};
+  for (std::size_t i = 0; i < N; ++i)
+    ret.arrays[i][0] = std::data(init)[i];
+
+  return ret;
+} ///< @private
+
 /** @name cotila::matrix deduction guides */
 ///@{
 

include/cotila/matrix/operators.h

@@ -21,7 +21,7 @@ constexpr bool operator==(const matrix<T, N, M> &a,
                           const matrix<T, N, M> &b) {
   for (std::size_t i = 0; i < N; ++i) {
     for (std::size_t j = 0; j < M; ++j) {
-      if (a[i][j] != b[i][j])
+      if (a.at(i, j) != b.at(i, j))
         return false;
     }
   }
@@ -70,7 +70,7 @@ constexpr matrix<T, N, M> operator+(T a, const matrix<T, N, M> &m) {
  *  @param b an \f$ N \times M \f$ matrix of type T
  *  @return \f$ \textbf{a} + \textbf{b} \f$ such that \f$ \left(\textbf{a} + \textbf{b}\right)_{ij} = \textbf{a}_{ij} + \textbf{b}_{ij} \f$
  *
- *  Computes the vector sum.
+ *  Computes the matrix sum.
  */
 template <typename T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> operator+(const matrix<T, N, M> &a,
@@ -83,7 +83,7 @@ constexpr matrix<T, N, M> operator+(const matrix<T, N, M> &a,
  *  @param a a scalar of type T
  *  @return \f$ \textbf{m}a \f$ such that \f$ \left(\textbf{m} a\right)_{ij} = \textbf{m}_{ij} a \f$
  *
- *  Computes the sum of a matrix and a scalar.
+ *  Computes the product of a matrix and a scalar.
  */
 template <typename T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> operator*(const matrix<T, N, M> &m, T a) {
@@ -95,7 +95,7 @@ constexpr matrix<T, N, M> operator*(const matrix<T, N, M> &m, T a) {
  *  @param m an \f$ N \times M \f$ matrix of type T
  *  @return \f$ a\textbf{m} \f$ such that \f$ \left(a\textbf{m}\right)_{ij} = a\textbf{m}_{ij} \f$
  *
- *  Computes the sum of a matrix and a scalar.
+ *  Computes the product of a matrix and a scalar.
  */
 template <typename T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> operator*(T a, const matrix<T, N, M> &m) {
@@ -107,7 +107,7 @@ constexpr matrix<T, N, M> operator*(T a, const matrix<T, N, M> &m) {
  *  @param b an \f$ N \times M \f$ matrix of type T
  *  @return \f$ \textbf{a} \circ \textbf{b} \f$ such that \f$ \left(\textbf{a} \circ \textbf{b}\right)_{ij} = \textbf{a}_{ij} \textbf{b}_{ij} \f$
  *
- *  Computes the Hadamard, or elementwise, product of two vectors.
+ *  Computes the Hadamard, or elementwise, product of two matrices.
  */
 template <typename T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> operator*(const matrix<T, N, M> &a,