From a8b16d8c9729daf0047183ca660f9bc5fb613f74 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Sun, 10 Oct 2021 21:46:04 +0100 Subject: [PATCH 1/7] Create compatibility between vectors and matrices Vectors are now children of matrices, so they share the same set of properties. This hierarchy change allows the deletions of some duplicate code that was shared among the two classes. Some algorithms have been rewritten to accept both vectors and matrices. --- include/cotila/cotila.h | 1 - include/cotila/detail/tmp.h | 9 ++ include/cotila/matrix/math.h | 58 +++++++++++- include/cotila/matrix/matrix.h | 91 ++++++++++++++++++- include/cotila/matrix/operators.h | 2 +- include/cotila/matrix/utility.h | 87 ++++++++++++++++-- include/cotila/scalar/math.h | 10 ++- include/cotila/vector/math.h | 89 ------------------ include/cotila/vector/operators.h | 144 ------------------------------ include/cotila/vector/utility.h | 119 +----------------------- include/cotila/vector/vector.h | 80 +++++------------ test/vector_test.h | 2 +- 12 files changed, 264 insertions(+), 428 deletions(-) delete mode 100644 include/cotila/vector/operators.h diff --git a/include/cotila/cotila.h b/include/cotila/cotila.h index 6623a4c..5c7ec10 100644 --- a/include/cotila/cotila.h +++ b/include/cotila/cotila.h @@ -7,7 +7,6 @@ #include #include #include -#include #include #include diff --git a/include/cotila/detail/tmp.h b/include/cotila/detail/tmp.h index de7108b..d47bbbe 100644 --- a/include/cotila/detail/tmp.h +++ b/include/cotila/detail/tmp.h @@ -19,6 +19,15 @@ template struct all_same_value { static constexpr T value = v; }; +template +static auto test_callable(int) + -> decltype(void(std::declval()(std::declval())), std::true_type{}); +template +static auto test_callable(long) -> std::false_type; + +template +struct is_callable_with_a_size_t : decltype(detail::test_callable(0)){}; + } // namespace detail } // namespace cotila diff --git a/include/cotila/matrix/math.h b/include/cotila/matrix/math.h index c69f042..b4fb568 100644 --- a/include/cotila/matrix/math.h +++ b/include/cotila/matrix/math.h @@ -7,10 +7,9 @@ #include #include -#include -#include #include #include +#include #include namespace cotila { @@ -31,7 +30,7 @@ constexpr matrix conj(const matrix &m) { return elementwise(cotila::conj, 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 &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 &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 +constexpr matrix, N, M> abs(const matrix &v) { + return elementwise(abs, v); +} + +/** @brief computes the sum of elements of a matrix + * @param m an \f$ M \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 constexpr T sum(const matrix &v) { + return accumulate(v, T{}, std::plus()); +} + +/** @brief computes the minimum valued element + * @param m an \f$ M \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 constexpr T min(const matrix &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$ M \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 constexpr T max(const matrix &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 @@ -285,6 +325,16 @@ constexpr T trace(const matrix &m) { return sum(generate([&m](std::size_t i){ return m[i][i]; })); } +/** @brief computes the elementwise square root + * @param m an \f$ M \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 +constexpr matrix sqrt(const matrix &v) { + return elementwise(static_cast(sqrt), v); +} /** }@*/ } // namespace cotila diff --git a/include/cotila/matrix/matrix.h b/include/cotila/matrix/matrix.h index 49cc2bb..3650401 100644 --- a/include/cotila/matrix/matrix.h +++ b/include/cotila/matrix/matrix.h @@ -6,7 +6,7 @@ #define COTILA_MATRIX_MATRIX_H_ #include -#include +#include #include #include #include @@ -31,6 +31,7 @@ template 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 */ ///@{ @@ -58,6 +59,27 @@ template struct matrix { return generate([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) { + if (row >= N || column >= M) + throw "index out of range"; + return arrays[row][column]; + } + + constexpr const T &at(std::size_t row, std::size_t column = 0) const { + if (row >= N || column >= M) + throw "index out of range"; + return arrays[row][column]; + } + /** @brief access specified element * @param i index of the row * @return pointer to the specified row @@ -67,12 +89,63 @@ template 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 + typename std::enable_if::type + constexpr operator[](std::size_t i) { return arrays[i][0]; } + + /// @copydoc operator[] + template + typename std::enable_if::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 + typename std::enable_if::type + constexpr operator[](std::size_t i) { return arrays[i]; } + + /// @copydoc operator[] + template + typename std::enable_if::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() const { + vector ret{}; + for (size_type i = 0; i < column_size; ++i) + ret[i] = arrays[i][0]; + return ret; + } + + template + friend constexpr matrix from_initializer(Container &); + T arrays[N][M]; ///< @private }; @@ -80,6 +153,18 @@ template struct matrix { * @{ */ +template +constexpr matrix from_initializer(Container & init) { + if(std::size(init) != N) + throw "The initializer has not the same size than the vector"; + + matrix 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 */ ///@{ diff --git a/include/cotila/matrix/operators.h b/include/cotila/matrix/operators.h index 26e84a6..325121b 100644 --- a/include/cotila/matrix/operators.h +++ b/include/cotila/matrix/operators.h @@ -21,7 +21,7 @@ constexpr bool operator==(const matrix &a, const matrix &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; } } diff --git a/include/cotila/matrix/utility.h b/include/cotila/matrix/utility.h index dea8593..7d9d912 100644 --- a/include/cotila/matrix/utility.h +++ b/include/cotila/matrix/utility.h @@ -1,7 +1,8 @@ #ifndef COTILA_MATRIX_UTILITY_H_ #define COTILA_MATRIX_UTILITY_H_ -#include +#include +#include #include namespace cotila { @@ -10,6 +11,23 @@ namespace cotila { * @{ */ + +/** @brief accumulates an operation across a matrix + * @param m an \f$ M \times M \f$ matrix of type T + * @param init the initial value + * @param f a function of type F that operates between U and matrix elements of type T + * @return \f$ f\left(f\left(\ldots f\left(\textrm{init}, \textbf{v}_1\right), \ldots\right), \textbf{v}_N \right) \f$ + * + * Accumulates an operation over the elements. This is equivalent to a functional fold. + */ +template +constexpr U accumulate(const matrix &m, U init, F &&f) { + U r = init; + for (const auto val : m) + r = std::apply(std::forward(f), std::forward_as_tuple(r, val)); + return r; +} + /** @brief applies a function elementwise between many matrices * @param f a function of type F that operates on many scalars of type T and returns a scalar of type U * @param m an \f$ N \times M \f$ matrix of type T @@ -30,14 +48,14 @@ constexpr matrix elementwise(F f, const matrix &m, matrix op_applied = {}; for (std::size_t i = 0; i < N; ++i) { for (std::size_t j = 0; j < M; ++j) { - op_applied[i][j] = - std::apply(f, std::forward_as_tuple(m[i][j], matrices[i][j]...)); + op_applied.at(i, j) = + std::apply(f, std::forward_as_tuple(m.at(i, j), matrices.at(i, j)...)); } } return op_applied; } -/** @brief casts a vector to another type +/** @brief casts a matrix to another type * @param m an \f$ N \times M \f$ matrix of type U * @return an \f$ N \times M \f$ matrix of type T containing the casted elements of \f$ \textbf{M} \f$ * @@ -54,29 +72,71 @@ constexpr matrix cast(const matrix &m) { * * Generates a matrix as a function of its indices. */ -template +template constexpr decltype(auto) generate(F &&f) { - matrix, N, M> generated = + + const auto global_f = [&f](std::size_t i, std::size_t j = 0) { + if constexpr (detail::is_callable_with_a_size_t::value) + return f(i); + else + return f(i, j); + }; + + matrix, N, M> generated = {}; for (std::size_t i = 0; i < N; ++i) { for (std::size_t j = 0; j < M; ++j) { - generated[i][j] = std::apply(f, std::forward_as_tuple(i, j)); + generated.at(i, j) = std::apply(global_f, std::forward_as_tuple(i, j)); } } return generated; } +/** @brief generates a matrix containing consecutive elements + * @param value the value of the first element of the matrix + * @return an \f$ M \times N \f$ matrix of type T, with consecutive elements spaced by 1. + * Generates a matrix containing consecutive elements spaced by 1. + */ +template +constexpr matrix iota(T value = T()) { + matrix seq = {}; + for (auto &x : seq) { + x = value; + value += 1; // equivalent to value++, see GCC Bug 91705 + } + return seq; +} + /** @brief generates a matrix containing a single value * @param value the scalar value of all elements * @return an \f$ N \times M \f$ matrix of type T such that \f$ \textbf{m}_{ij} = \textrm{value}\ \forall i,j \f$ * * Generates a matrix with all elements equal to a single value. */ -template +template constexpr matrix fill(T value) { return generate([value](std::size_t, std::size_t) { return value; }); } +/** @brief shifts matrix elements + * @param m an \f$ M \times M \f$ matrix of type T + * @param n the amount to shift each element + * @return an N-vector of type T \f$ \textbf{v} \gg n \f$ such that + * \f$ \left(\textbf{v} \gg n\right)_i = \textbf{v}_{(i + n)\ \textrm{mod}\ N} \f$ + * + * Rotates a matrix by shifting its elements vertically. + */ +template +constexpr matrix rotate(matrix m, int n) { + matrix rotated = {}; + // add N (the modulus) to n until it is positive + while (n < 0) + n += N; + for (std::size_t i = 0; i < N; ++i) + rotated[i] = m[(i + n) % N]; + return rotated; +} + /** @brief the matrix identity * * The matrix identity \f$ I_N \f$. @@ -227,6 +287,17 @@ constexpr matrix as_row(const vector &v) { return generate<1, N>([&v](auto, auto j) { return v[j]; }); } +/** @brief generates a matrix of equally spaced elements + * @param min the value of the first element of the matrix + * @param max the value of the last element of the matrix + * @return an \f$ M \times N \f$ matrix of type T, with equally spaced elements. + * Generates a matrix of equally spaced elements. + */ +template +constexpr matrix linspace(T min, T max) { + return ((max - min) / ((N * M) - 1)) * iota() + min; +} + /** }@*/ } // namespace cotila diff --git a/include/cotila/scalar/math.h b/include/cotila/scalar/math.h index c7ff5e8..9de4e26 100644 --- a/include/cotila/scalar/math.h +++ b/include/cotila/scalar/math.h @@ -47,7 +47,10 @@ constexpr float sqrt(float x) { return float(sqrt(double(x))); } * * Computes the absolute value. */ -template constexpr detail::remove_complex_t abs(T x) { +template +typename std::enable_if_t< + std::is_arithmetic_v>, detail::remove_complex_t> +constexpr abs(T x) { COTILA_DETAIL_ASSERT_ARITHMETIC(T); if constexpr (detail::is_complex_v) return sqrt(x.real() * x.real() + x.imag() * x.imag()); @@ -108,7 +111,10 @@ constexpr double nthroot(double x, int n) { * * Computes the complex conjugate. */ -template constexpr T conj(T x) { +template +typename std::enable_if_t< + std::is_arithmetic_v>, T> +constexpr conj(T x) { COTILA_DETAIL_ASSERT_ARITHMETIC(T); if constexpr (detail::is_complex_v) return {x.real(), -x.imag()}; diff --git a/include/cotila/vector/math.h b/include/cotila/vector/math.h index c69fdd5..c3d6d19 100644 --- a/include/cotila/vector/math.h +++ b/include/cotila/vector/math.h @@ -16,64 +16,6 @@ namespace cotila { * @{ */ -/** @brief computes the elementwise complex conjugate - * @param v an N-vector of type T - * @return an N-vector \f$ \overline{\textbf{v}} \f$ of type T such that - * \f$ \left(\overline{\textbf{v}}\right)_i = \overline{v_i} \f$ - * - * Computes the elementwise complex conjugate of a vector. - */ -template -constexpr vector conj(const vector &v) { - return elementwise(cotila::conj, v); -} - -/** @brief computes the elementwise square root - * @param v an N-vector of type T - * @return an N-vector \f$ \begin{bmatrix} \sqrt{v_1} & \ldots &\sqrt{v_N} \end{bmatrix} \f$ of type T - * - * Computes the elementwise square root of a vector. - */ -template -constexpr vector sqrt(const vector &v) { - return elementwise(static_cast(sqrt), v); -} - -/** @brief computes the elementwise real - * @param v an N-vector of type T - * @return an N-vector \f$ \textbf{u} \f$ of type T such that - * \f$ \textbf{u}_i = \mathbb{R}\{\textbf{v}_i\} \f$ - * - * Computes the elementwise real of a vector. - */ -template -constexpr vector, N> real(const vector &v) { - return elementwise([](auto i) { return std::real(i); }, v); -} - -/** @brief computes the elementwise imag - * @param v an N-vector of type T - * @return an N-vector \f$ \textbf{u} \f$ of type T such that - * \f$ \textbf{u}_i = \mathbb{I}\{\textbf{v}_i\} \f$ - * - * Computes the elementwise imag of a vector. - */ -template -constexpr vector, N> imag(const vector &v) { - return elementwise([](auto i) { return std::imag(i); }, v); -} - -/** @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 vector. - */ -template -constexpr vector, N> abs(const vector &v) { - return elementwise(abs, v); -} - /** @brief computes the dot product * @param a an N-vector of type T * @param b an N-vector of type T @@ -90,37 +32,6 @@ constexpr T dot(const vector &a, const vector &b) { return r; } -/** @brief computes the sum of elements - * @param v an N-vector of type T - * @return a scalar \f$ \sum\limits_{i} v_i \f$ of type T - * - * Computes the sum of the elements of a vector. - */ -template constexpr T sum(const vector &v) { - return accumulate(v, static_cast(0), std::plus()); -} - -/** @brief computes the minimum valued element - * @param v an N-vector 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 vector. - */ -template constexpr T min(const vector &v) { - return accumulate(v, v[0], [](T a, T b) { return std::min(a, b); }); -} - -/** @brief computes the maximum valued element - * @param v an N-vector 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 vector. - */ -template constexpr T max(const vector &v) { - return accumulate(v, v[0], [](T a, T b) { return std::max(a, b); }); -} - - /** @brief computes the index of the minimum valued element * @param v an N-vector of type T * @return an index \f$ i \f$ where \f$ v_i \leq v_j,\ \forall j \f$ diff --git a/include/cotila/vector/operators.h b/include/cotila/vector/operators.h deleted file mode 100644 index b75666a..0000000 --- a/include/cotila/vector/operators.h +++ /dev/null @@ -1,144 +0,0 @@ -#ifndef COTILA_VECTOR_OPERATORS_H_ -#define COTILA_VECTOR_OPERATORS_H_ - -#include - -namespace cotila { - -/** \addtogroup vector - * @{ - */ - -/** @brief checks equality of two vectors - * @param a an N-vector of type T - * @param b an N-vector of type T - * @return true if and only if \f$ \textbf{a}_i = \textbf{b}_i\ \forall i \in 1\ .. N \f$ - * - * Checks the equality of two vectors. - */ -template -constexpr bool operator==(const vector &a, const vector &b) { - for (std::size_t i = 0; i < vector::size; ++i) { - if (a[i] != b[i]) - return false; - } - return true; -} - -/** @brief checks inequality of two vectors - * @param a an N-vector of type T - * @param b an N-vector of type T - * @return false if and only if \f$ \textbf{a}_i = \textbf{b}_i\ \forall i \in 1\ .. N \f$ - * - * Checks the inequality of two vectors. - */ -template -constexpr bool operator!=(const vector &a, const vector &b) { - return !(a == b); -} - -/** @brief computes the sum of a vector and a scalar - * @param v an N-vector of type T - * @param a a scalar of type T - * @return \f$ \textbf{v} + a \f$ such that \f$ \left(\textbf{v} + a\right)_i = \textbf{v}_i + a \f$ - * - * Computes the sum of a vector and a scalar. - */ -template -constexpr vector operator+(const vector &v, T a) { - return elementwise([a](T x) { return x + a; }, v); -} - -/** @brief computes the sum of a vector and a scalar - * @param a a scalar of type T - * @param v an N-vector of type T - * @return \f$ a + \textbf{v} \f$ such that \f$ \left(a + \textbf{v}\right)_i = a + \textbf{v}_i \f$ - * - * Computes the sum of a vector and a scalar. - */ -template -constexpr vector operator+(T a, const vector &v) { - return elementwise([a](T x) { return a + x; }, v); -} - -/** @brief computes the vector sum - * @param a an N-vector of type T - * @param b an N-vector of type T - * @return \f$ \textbf{a} + \textbf{b} \f$ such that \f$ \left(\textbf{a} + \textbf{b}\right)_i = \textbf{a}_i + \textbf{b}_i \f$ - * - * Computes the vector sum. - */ -template -constexpr vector operator+(const vector &a, - const vector &b) { - return elementwise(std::plus(), a, b); -} - -/** @brief computes the product of a vector and a scalar - * @param v an N-vector of type T - * @param a a scalar of type T - * @return \f$ \textbf{v}a \f$ such that \f$ \left(\textbf{v}a\right)_i = \textbf{v}_i a \f$ - * - * Computes the product of a vector and a scalar. - */ -template -constexpr vector operator*(const vector &v, T a) { - return elementwise([a](T x) { return x * a; }, v); -} - -/** @brief computes the product of a vector and a scalar - * @param a a scalar of type T - * @param v an N-vector of type T - * @return \f$ a\textbf{v} \f$ such that \f$ \left(a\textbf{v}\right)_i = a\textbf{v}_i \f$ - * - * Computes the product of a vector and a scalar. - */ -template -constexpr vector operator*(T a, const vector &v) { - return elementwise([a](T x) { return a * x; }, v); -} - -/** @brief computes the Hadamard product - * @param a an N-vector of type T - * @param b an N-vector of type T - * @return \f$ \textbf{a} \circ \textbf{b} \f$ such that \f$ \left(\textbf{a} \circ \textbf{b}\right)_i = \textbf{a}_i \textbf{b}_i \f$ - * - * Computes the Hadamard, or elementwise, product of two vectors. - */ -template -constexpr vector operator*(const vector &a, - const vector &b) { - return elementwise(std::multiplies(), a, b); -} - -/** @brief computes the quotient between a vector and a scalar - * @param v an N-vector of type T - * @param a a scalar of type T - * @return \f$ \textbf{v}/a \f$ such that \f$ \left(\textbf{v}/a\right)_i = \frac{\textbf{v}_i}{a} \f$ - * - * Computes division between a vector and a scalar. - */ -template -constexpr vector operator/(T a, const vector &v) { - return elementwise([a](T x) { return a / x; }, v); -} - -/** @brief computes the elementwise vector quotient - * @param a an N-vector of type T - * @param b an N-vector of type T - * @return \f$ \textbf{a} \circ \textbf{b}' \f$ such that \f$ {\textbf{b}_i}' = \left(\textbf{b}_i\right)^{-1}\f$ - * and \f$ \left(\textbf{a} \circ \textbf{b}'\right)_i = \textbf{a}_i {\textbf{b}'}_i \f$ - * - * Computes elementwise division between two vectors - */ -template -constexpr vector operator/(const vector &a, - const vector &b) { - return elementwise(std::divides(), a, b); -} - -/** }@*/ - -} // namespace cotila - -#endif // COTILA_VECTOR_OPERATORS_H_ diff --git a/include/cotila/vector/utility.h b/include/cotila/vector/utility.h index e55543b..adfcdc5 100644 --- a/include/cotila/vector/utility.h +++ b/include/cotila/vector/utility.h @@ -10,123 +10,6 @@ namespace cotila { * @{ */ -/** @brief applies a function elementwise between many vectors - * @param f a function of type F that operates on many scalars of type T and returns a scalar of type U - * @param v an N-vector of type T - * @param vectors additional N-vectors of type T - * @return an N-vector of type T with elements described by \f$ f\left(\textbf{v}_i, \ldots\right) \f$ - * - * Applies a function elementwise between many vectors. - */ -template < - typename F, typename T, typename... Vectors, - typename U = std::invoke_result_t, - std::size_t N = - detail::all_same_value::value> -constexpr vector elementwise(F f, const vector &v, - const Vectors &... vectors) { - vector op_applied = {}; - for (std::size_t i = 0; i < N; ++i) - op_applied[i] = std::apply(f, std::forward_as_tuple(v[i], vectors[i]...)); - return op_applied; -} - -/** @brief accumulates an operation across a vector - * @param v an N-vector of type T - * @param init the initial value - * @param f a function of type F that operates between U and vector elements of type T - * @return \f$ f\left(f\left(\ldots f\left(\textrm{init}, \textbf{v}_1\right), \ldots\right), \textbf{v}_N \right) \f$ - * - * Accumulates an operation over the elements. This is equivalent to a functional fold. - */ -template -constexpr U accumulate(const vector &v, U init, F &&f) { - U r = init; - for (std::size_t i = 0; i < vector::size; ++i) - r = std::apply(std::forward(f), std::forward_as_tuple(r, v[i])); - return r; -} - -/** @brief casts a vector to another type - * @param v an N-vector of type U - * @return an N-vector of type T containing the casted elements of \f$ \textbf{v} \f$ - * - * Casts a vector to another type by `static_cast`ing each element. - */ -template -constexpr vector cast(const vector &v) { - return elementwise([](const U u) { return static_cast(u); }, v); -} - -/** @brief generates a vector containing consecutive elements - * @param value the value of the first element of the vector - * @return an N-vector of type T such that \f$ \textbf{v}_i = \textrm{start} + i - 1 \f$ - * - * Generates a vector containing consecutive elements spaced by 1. - */ -template -constexpr vector iota(T value = T()) { - vector seq = {}; - for (auto &x : seq) { - x = value; - value += 1; // equivalent to value++, see GCC Bug 91705 - } - return seq; -} - -/** @brief generates a vector of equally spaced elements - * @param min the value of the first element of the vector - * @param max the value of the last element of the vector - * @return an N-vector of type T with elements spaced by \f$ \frac{\textbf{v}_N - \textbf{v}_1}{N-1} \f$ - * - * Generates a vector of equally spaced elements. - */ -template -constexpr vector linspace(T min, T max) { - return ((max - min) / (N - 1)) * iota() + min; -} - -/** @brief generates a vector containing a single value - * @param value the scalar value of all elements - * @return an N-vector of type T such that \f$ \textbf{v}_i = \textrm{value}\ \forall i \in 1\ .. N\f$ - * - * Generates a vector with all elements equal to a single value. - */ -template constexpr vector fill(T value) { - vector filled = {}; - for (auto &x : filled) - x = value; - return filled; -} - -/** @brief generates a vector as a function of its index - * @param f a function that operates on an integer index - * @return an N-vector with type matching the return type of f such that \f$ \textbf{v}_i = f(i) \f$ - * - * Generates a vector as a function of its index. - */ -template constexpr decltype(auto) generate(F &&f) { - return elementwise(f, iota()); -} - -/** @brief shifts vector elements - * @param v an N-vector of type T - * @param n the amount to shift each element - * @return an N-vector of type T \f$ \textbf{v} \gg n \f$ such that - * \f$ \left(\textbf{v} \gg n\right)_i = \textbf{v}_{(i + n)\ \textrm{mod}\ N} \f$ - * - * Rotates a vector by shifting its elements. - */ -template -constexpr vector rotate(vector v, int n) { - vector rotated = {}; - // add N (the modulus) to n until it is positive - while (n < 0) - n += N; - for (std::size_t i = 0; i < N; ++i) - rotated[i] = v[(i + n) % N]; - return rotated; -} /** @brief slices a vector into a subvector * @param v an N-vector of type T @@ -152,7 +35,7 @@ constexpr vector slice(vector v, std::size_t start = 0) { * @return an N+M-vector \f$ \left[\textbf{a} \textbf{b}\right] \f$ such that * \f$ \left(\left[\textbf{a} \textbf{b}\right]\right)_i = \begin{cases} \textbf{a}_i & i < N\\ \textbf{b}_{i - N} & i \ge N \end{cases} \f$ * - * Slices a vector into a subvector. + * Concatenates two vectors. */ template constexpr vector concat(vector a, vector b) { diff --git a/include/cotila/vector/vector.h b/include/cotila/vector/vector.h index e6293fb..3c62277 100644 --- a/include/cotila/vector/vector.h +++ b/include/cotila/vector/vector.h @@ -6,68 +6,34 @@ #define COTILA_VECTOR_VECTOR_H_ #include +#include #include #include -#include -#include namespace cotila { -/** @brief A container representing a vector - * @tparam T scalar type to contain - * @tparam N size of the vector - * - * `cotila::vector` is a container representing a vector. - * - * It is an aggregate type similar to `std::array`, and can be initialized with - * aggregate initialization or with the `cotila::make_vector` function. - */ -template struct vector { - static_assert(N != 0, "vector must contain at least one element"); - COTILA_DETAIL_ASSERT_ARITHMETIC(T) - - using value_type = T; - using size_type = std::size_t; - static constexpr size_type size = N; ///< @brief size of the vector - - /** @name Element access */ - ///@{ - /** @brief access specified element - * @param i position of the scalar element - * @return the requested scalar element - * - * Returns a reference to the scalar element in position `i`, without bounds checking. - */ - constexpr T &operator[](size_type i) noexcept { return array[i]; } - - /// @copydoc operator[] - constexpr const T &operator[](size_type i) const noexcept { return array[i]; } - ///@} - - /** @name Iterators */ - ///@{ - /** @brief returns an iterator to the beginning - * @return an iterator to the beginning - * - * Returns an iterator to the beginning of the vector. - */ - constexpr T *begin() noexcept { return array; } - - /** @brief returns an iterator to the end - * @return an iterator past the end - * - * Returns an iterator to the end of the vector. - */ - constexpr T *end() noexcept { return array + N; } - - /// @copydoc begin - constexpr const T *cbegin() const noexcept { return array; } - - /// @copydoc end - constexpr const T *cend() const noexcept { return array + N; } - ///@} - - T array[N]; ///< @private + +template +struct matrix; + +template +constexpr matrix from_initializer(Container &); + +template +struct vector : public matrix { + constexpr vector() = default; + constexpr vector(const vector &other) = default; + constexpr vector(vector &&other) noexcept = default; + constexpr vector &operator=(const vector &other) = default; + constexpr vector &operator=(vector &&other) noexcept = default; + ~vector() = default; + + constexpr vector(std::initializer_list init) : + matrix{from_initializer(init)} {} + + constexpr operator matrix() { + return matrix(*this); + } }; /** \addtogroup vector diff --git a/test/vector_test.h b/test/vector_test.h index 4515fd8..d76d12d 100644 --- a/test/vector_test.h +++ b/test/vector_test.h @@ -56,7 +56,7 @@ static_assert(max_index(vector{1, 2, 3}) == 2, "max_index"); static_assert(iota<5>(0) == vector{0, 1, 2, 3, 4}, "iota"); -static_assert(iota<5, double>() == vector{0., 1., 2., 3., 4.}, "iota"); +static_assert(iota<5, 1, double>() == vector{0., 1., 2., 3., 4.}, "iota"); static_assert(cast(vector{0, 1, 2}) == vector{0., 1., 2.}, "cast"); From f494e542a221b0a7b358553ecc095adf70c6a407 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Sun, 10 Oct 2021 22:05:32 +0100 Subject: [PATCH 2/7] Factorize matrix::at const and non-const implementation --- include/cotila/matrix/matrix.h | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/include/cotila/matrix/matrix.h b/include/cotila/matrix/matrix.h index 3650401..76da6d7 100644 --- a/include/cotila/matrix/matrix.h +++ b/include/cotila/matrix/matrix.h @@ -43,7 +43,7 @@ template struct matrix { */ constexpr vector 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([i, this](std::size_t j) { return arrays[i][j]; }); } @@ -55,7 +55,7 @@ template struct matrix { */ constexpr vector 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([i, this](std::size_t j) { return arrays[j][i]; }); } @@ -69,16 +69,22 @@ template struct matrix { * ease the writing of generic algorithm */ constexpr T &at(std::size_t row, std::size_t column = 0) { - if (row >= N || column >= M) - throw "index out of range"; - return arrays[row][column]; + return get_real_by_name_impl(*this, row, column); } + /** + * @copydoc at + */ constexpr const T &at(std::size_t row, std::size_t column = 0) const { + return get_real_by_name_impl(*this, row, column); + } + + template + static constexpr auto &get_real_by_name_impl(Self &self, std::size_t row, std::size_t column = 0) { if (row >= N || column >= M) throw "index out of range"; - return arrays[row][column]; - } + return self.arrays[row][column]; + } ///< @private /** @brief access specified element * @param i index of the row From f4e039102104ef492b58a2c49e982604b7153c81 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Sun, 10 Oct 2021 23:07:18 +0100 Subject: [PATCH 3/7] Use type derived from std::exception to throw --- include/cotila/matrix/math.h | 2 +- include/cotila/matrix/matrix.h | 10 +++++----- include/cotila/matrix/utility.h | 3 ++- include/cotila/scalar/math.h | 4 ++-- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/include/cotila/matrix/math.h b/include/cotila/matrix/math.h index b4fb568..2322960 100644 --- a/include/cotila/matrix/math.h +++ b/include/cotila/matrix/math.h @@ -309,7 +309,7 @@ constexpr T det(const matrix &m) { template constexpr matrix inverse(const matrix &m) { if (rank(m) < M) - throw "matrix is not invertible"; + throw std::logic_error("matrix is not invertible"); return submat(rref(horzcat(m, identity)), 0, M); } diff --git a/include/cotila/matrix/matrix.h b/include/cotila/matrix/matrix.h index 76da6d7..35e7ee4 100644 --- a/include/cotila/matrix/matrix.h +++ b/include/cotila/matrix/matrix.h @@ -69,20 +69,20 @@ template struct matrix { * ease the writing of generic algorithm */ constexpr T &at(std::size_t row, std::size_t column = 0) { - return get_real_by_name_impl(*this, row, column); + return at_impl(*this, row, column); } /** * @copydoc at */ constexpr const T &at(std::size_t row, std::size_t column = 0) const { - return get_real_by_name_impl(*this, row, column); + return at_impl(*this, row, column); } template - static constexpr auto &get_real_by_name_impl(Self &self, std::size_t row, std::size_t column = 0) { + static constexpr auto &at_impl(Self &self, std::size_t row, std::size_t column = 0) { if (row >= N || column >= M) - throw "index out of range"; + throw std::out_of_range("Indexes out of range"); return self.arrays[row][column]; } ///< @private @@ -162,7 +162,7 @@ template struct matrix { template constexpr matrix from_initializer(Container & init) { if(std::size(init) != N) - throw "The initializer has not the same size than the vector"; + throw std::invalid_argument("The initializer has not the same size than the vector"); matrix ret{}; for (std::size_t i = 0; i < N; ++i) diff --git a/include/cotila/matrix/utility.h b/include/cotila/matrix/utility.h index 7d9d912..f4abe04 100644 --- a/include/cotila/matrix/utility.h +++ b/include/cotila/matrix/utility.h @@ -243,7 +243,8 @@ constexpr matrix vertcat(const matrix &a, const matrix constexpr matrix submat(const matrix &m, std::size_t a, std::size_t b){ - if ((a + P > M) || (b + Q > N)) throw "index out of range"; + if ((a + P > M) || (b + Q > N)) + throw std::out_of_range("Index out of range"); return generate([&m, &a, &b](std::size_t i, std::size_t j){ return m[a + i][b + j]; }); diff --git a/include/cotila/scalar/math.h b/include/cotila/scalar/math.h index 9de4e26..0abf528 100644 --- a/include/cotila/scalar/math.h +++ b/include/cotila/scalar/math.h @@ -23,7 +23,7 @@ namespace cotila { */ constexpr double sqrt(double x) { if (x < 0) - throw "sqrt argument must be positive"; + throw std::domain_error("sqrt argument must be positive"); double prev = 0; double est = (1 + x) / 2; while (prev != est) { @@ -94,7 +94,7 @@ constexpr double exponentiate(double x, int n) { */ constexpr double nthroot(double x, int n) { if (x < 0) - throw "nth root argument must be positive"; + throw std::domain_error("nth root argument must be positive"); double prev = -1; double est = 1; while (prev != est) { From 413be98dc37611cd8e326990a4ce5f73c4c7b47b Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Sun, 10 Oct 2021 23:11:35 +0100 Subject: [PATCH 4/7] Remove typos and update documentation --- include/cotila/matrix/math.h | 8 ++++---- include/cotila/matrix/operators.h | 8 ++++---- include/cotila/matrix/utility.h | 8 ++++---- include/cotila/vector/vector.h | 7 +++++++ 4 files changed, 19 insertions(+), 12 deletions(-) diff --git a/include/cotila/matrix/math.h b/include/cotila/matrix/math.h index 2322960..003f788 100644 --- a/include/cotila/matrix/math.h +++ b/include/cotila/matrix/math.h @@ -68,7 +68,7 @@ constexpr matrix, N, M> abs(const matrix &v } /** @brief computes the sum of elements of a matrix - * @param m an \f$ M \times M \f$ matrix of type T + * @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. @@ -78,7 +78,7 @@ template constexpr T sum(const matrix } /** @brief computes the minimum valued element - * @param m an \f$ M \times M \f$ matrix of type T + * @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. @@ -88,7 +88,7 @@ template constexpr T min(const matrix } /** @brief computes the maximum valued element - * @param m an \f$ M \times M \f$ matrix of type T + * @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. @@ -326,7 +326,7 @@ constexpr T trace(const matrix &m) { } /** @brief computes the elementwise square root - * @param m an \f$ M \times M \f$ matrix of type T + * @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. diff --git a/include/cotila/matrix/operators.h b/include/cotila/matrix/operators.h index 325121b..bca22ac 100644 --- a/include/cotila/matrix/operators.h +++ b/include/cotila/matrix/operators.h @@ -70,7 +70,7 @@ constexpr matrix operator+(T a, const matrix &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 constexpr matrix operator+(const matrix &a, @@ -83,7 +83,7 @@ constexpr matrix operator+(const matrix &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 constexpr matrix operator*(const matrix &m, T a) { @@ -95,7 +95,7 @@ constexpr matrix operator*(const matrix &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 constexpr matrix operator*(T a, const matrix &m) { @@ -107,7 +107,7 @@ constexpr matrix operator*(T a, const matrix &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 constexpr matrix operator*(const matrix &a, diff --git a/include/cotila/matrix/utility.h b/include/cotila/matrix/utility.h index f4abe04..96e651c 100644 --- a/include/cotila/matrix/utility.h +++ b/include/cotila/matrix/utility.h @@ -13,7 +13,7 @@ namespace cotila { /** @brief accumulates an operation across a matrix - * @param m an \f$ M \times M \f$ matrix of type T + * @param m an \f$ N \times M \f$ matrix of type T * @param init the initial value * @param f a function of type F that operates between U and matrix elements of type T * @return \f$ f\left(f\left(\ldots f\left(\textrm{init}, \textbf{v}_1\right), \ldots\right), \textbf{v}_N \right) \f$ @@ -119,10 +119,10 @@ constexpr matrix fill(T value) { } /** @brief shifts matrix elements - * @param m an \f$ M \times M \f$ matrix of type T + * @param m an \f$ N \times M \f$ matrix of type T * @param n the amount to shift each element - * @return an N-vector of type T \f$ \textbf{v} \gg n \f$ such that - * \f$ \left(\textbf{v} \gg n\right)_i = \textbf{v}_{(i + n)\ \textrm{mod}\ N} \f$ + * @return an \f$ N \times M \f$ matrix of type T such that + * \f$ \left(\textbf{m} \gg n\right)_{ij} = \textbf{v}_{kj} with k = (i + n)\ \textrm{mod}\ N \f$ * * Rotates a matrix by shifting its elements vertically. */ diff --git a/include/cotila/vector/vector.h b/include/cotila/vector/vector.h index 3c62277..4dca372 100644 --- a/include/cotila/vector/vector.h +++ b/include/cotila/vector/vector.h @@ -19,6 +19,13 @@ struct matrix; template constexpr matrix from_initializer(Container &); +/** @brief A container representing a vector + * @tparam T scalar type to contain + * @tparam N number of rows + * + * `cotila::vector` is a container representing a vector. + * It inherits all its properties from `̀€cotila::matrix`. + */ template struct vector : public matrix { constexpr vector() = default; From 5fd74b0d851972485fbb6258d18efb1eb73204f2 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Mon, 11 Oct 2021 01:45:01 +0100 Subject: [PATCH 5/7] Add a test for dot --- test/vector_test.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/vector_test.h b/test/vector_test.h index d76d12d..33e2673 100644 --- a/test/vector_test.h +++ b/test/vector_test.h @@ -54,6 +54,8 @@ static_assert(min_index(vector{1, 2, 3}) == 0, "min_index"); static_assert(max_index(vector{1, 2, 3}) == 2, "max_index"); +static_assert(dot(vector{1, 2, 3}, vector{4, 5, 6}) == 32, "dot"); + static_assert(iota<5>(0) == vector{0, 1, 2, 3, 4}, "iota"); static_assert(iota<5, 1, double>() == vector{0., 1., 2., 3., 4.}, "iota"); From 834533f1508c94a96a2173c8554265ffa0620c77 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Mon, 11 Oct 2021 01:56:56 +0100 Subject: [PATCH 6/7] Fix rotate for matrices --- include/cotila/matrix/utility.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/cotila/matrix/utility.h b/include/cotila/matrix/utility.h index 96e651c..832018f 100644 --- a/include/cotila/matrix/utility.h +++ b/include/cotila/matrix/utility.h @@ -133,7 +133,8 @@ constexpr matrix rotate(matrix m, int n) { while (n < 0) n += N; for (std::size_t i = 0; i < N; ++i) - rotated[i] = m[(i + n) % N]; + for (std::size_t j = 0; j < M; ++j) + rotated.at(i, j) = m.at((i + n) % N, j); return rotated; } From 700014adecbf692275d02cdefa7050101d11edd8 Mon Sep 17 00:00:00 2001 From: Lucas CHOLLET Date: Mon, 11 Oct 2021 01:57:13 +0100 Subject: [PATCH 7/7] Add test for rotate --- test/matrix_test.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/matrix_test.h b/test/matrix_test.h index 0ebd6ad..2c894d8 100644 --- a/test/matrix_test.h +++ b/test/matrix_test.h @@ -4,8 +4,7 @@ #include #include -namespace cotila { -namespace test { +namespace cotila::test { // Uniform initialization constexpr matrix m1 = { @@ -42,6 +41,12 @@ static_assert(cast(m1) == static_assert(fill<2, 2>(3.) == matrix{{{3., 3.}, {3., 3.}}}, "matrix fill"); +static_assert(rotate(m1, 1) == matrix {4, 5, 6, 7, 8, 9, 1, 2, 3}, "rotate left"); + +static_assert(rotate(m1, -1) == matrix {7, 8, 9, 1, 2, 3, 4, 5, 6}, "rotate right"); + +static_assert(rotate(m1, -1) == rotate(m1, 2), "rotate modulo"); + static_assert(transpose(m1) == matrix{{{1., 4., 7.}, {2., 5., 8.}, @@ -143,7 +148,6 @@ static_assert(imag(m1c) == matrix{{{1., 0., 0.}, {1., 0., 0.}, {1., 0., 0.}}}, "imag"); -} // namespace test } // namespace cotila #endif // COTILA_MATRIX_TEST_H_