diff options
| -rw-r--r-- | CMakeLists.txt | 1 | ||||
| -rw-r--r-- | include/mm/mmmatrix.hpp | 374 | ||||
| -rw-r--r-- | mmmatrix.cpp | 1 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 32 | ||||
| -rw-r--r-- | test/matrix_example.cpp | 54 |
5 files changed, 462 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index ed126da..7a14759 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,6 +12,7 @@ list(INSERT CMAKE_MODULE_PATH 0 ${CMAKE_SOURCE_DIR}/cmake) # build libmm add_library(mm STATIC ${CMAKE_CURRENT_SOURCE_DIR}/mmvec.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/mmmatrix.cpp ) add_library(MM::mm ALIAS mm) diff --git a/include/mm/mmmatrix.hpp b/include/mm/mmmatrix.hpp new file mode 100644 index 0000000..5285429 --- /dev/null +++ b/include/mm/mmmatrix.hpp @@ -0,0 +1,374 @@ +/* mmmatrix.hpp + * Part of Mathematical library built (ab)using Modern C++ 17 abstractions. + * + * This library is not intended to be _performant_, it does not contain + * hand written SMID / SSE / AVX optimizations. It is instead an example + * of highly inefficient (but abstract!) code, where matrices can contain any + * data type. + * + * Naoki Pross <naopross@thearcway.org> + * 2018 ~ 2019 + */ +#pragma once + +#include <iostream> +#include <iomanip> +#include <cstring> +#include <cassert> +#include <initializer_list> +#include <array> + +namespace mm { + template<typename T, std::size_t Rows, std::size_t Cols> + class basic_matrix; + + /* specialization of basic_matrx for Cols = 1 */ + template<typename T, std::size_t Rows> + class row_vec; + + /* specialization of basic_matrx for Rows = 1 */ + template<typename T, std::size_t Cols> + class col_vec; + + /* shorter name for basic_matrix */ + template<typename T, std::size_t Rows, std::size_t Cols> + class matrix; + + /* specialization of basic_matrix for Rows == Cols */ + template<typename T, std::size_t N> + class square_matrix; + + // template<typename T, std::size_t N> + // class diag_matrix; +} + +template<typename T, std::size_t Rows, std::size_t Cols> +class mm::basic_matrix { +public: + using type = T; + + template<typename U, std::size_t ORows, std::size_t OCols> + friend class mm::basic_matrix; + + static constexpr std::size_t rows = Rows; + static constexpr std::size_t cols = Cols; + + basic_matrix(); + + // from initializer_list + basic_matrix(std::initializer_list<std::initializer_list<T>> l); + + // copyable and movable + basic_matrix(const basic_matrix<T, Rows, Cols>& other); + basic_matrix(basic_matrix<T, Rows, Cols>&& other); + + // copy from another matrix + template<std::size_t ORows, std::size_t OCols> + basic_matrix(const basic_matrix<T, ORows, OCols>& other); + + // access data + T& at(std::size_t row, std::size_t col); + const T& at(std::size_t row, std::size_t col) const; + // allows to access a matrix M at row j col k with M[j][k] + auto operator[](std::size_t index); + + void swap_rows(std::size_t x, std::size_t y); + void swap_cols(std::size_t x, std::size_t y); + + // mathematical operations + virtual basic_matrix<T, Cols, Rows> transposed() const; + inline basic_matrix<T, Cols, Rows> td() const { return transposed(); } + + // bool is_invertible() const; + // basic_matrix<T, Rows, Cols> inverse() const; + + + /// downcast to square matrix + static inline constexpr bool is_square() { return (Rows == Cols); } + inline constexpr square_matrix<T, Rows> to_square() const { + static_assert(is_square()); + return static_cast<square_matrix<T, Rows>>(*this); + } + + + /// downcast to row_vector + static inline constexpr bool is_row_vec() { return (Cols == 1); } + inline constexpr row_vec<T, Rows> to_row_vec() const { + static_assert(is_row_vec()); + return static_cast<row_vec<T, Rows>>(*this); + } + + /// downcast to col_vector + static inline constexpr bool is_col_vec() { return (Rows == 1); } + inline constexpr col_vec<T, Cols> to_col_vec() const { + static_assert(is_col_vec()); + return static_cast<col_vec<T, Cols>>(*this); + } + +protected: + template<typename ConstIterator> + basic_matrix(ConstIterator begin, ConstIterator end); + +private: + std::array<T, Rows * Cols> data; +}; + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols>::basic_matrix() { + std::fill(data.begin(), data.end(), 0); +} + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols>::basic_matrix( + std::initializer_list<std::initializer_list<T>> l +) { + assert(l.size() == Rows); + auto data_it = data.begin(); + + for (auto&& row : l) { + data_it = std::copy(row.begin(), row.end(), data_it); + } +} + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols>::basic_matrix( + const mm::basic_matrix<T, Rows, Cols>& other +) : data(other.data) {} + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols>::basic_matrix( + mm::basic_matrix<T, Rows, Cols>&& other +) : data(std::forward<decltype(other.data)>(other.data)) {} + +template<typename T, std::size_t Rows, std::size_t Cols> +template<std::size_t ORows, std::size_t OCols> +mm::basic_matrix<T, Rows, Cols>::basic_matrix( + const mm::basic_matrix<T, ORows, OCols>& other +) { + static_assert((ORows <= Rows), + "cannot copy a taller matrix into a smaller one" + ); + + static_assert((OCols <= Cols), + "cannot copy a larger matrix into a smaller one" + ); + + std::fill(data.begin(), data.end(), 0); + for (unsigned row = 0; row < Rows; row++) + for (unsigned col = 0; col < Cols; col++) + this->at(row, col) = other.at(row, col); +} + +/* protected construtor */ +template<typename T, std::size_t Rows, std::size_t Cols> +template<typename ConstIterator> +mm::basic_matrix<T, Rows, Cols>::basic_matrix( + ConstIterator begin, ConstIterator end +) { + assert(static_cast<unsigned>(std::distance(begin, end)) >= ((Rows * Cols))); + std::copy(begin, end, data.begin()); +} + + +/* member functions */ + +template<typename T, std::size_t Rows, std::size_t Cols> +T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) { + assert(row < Rows); // "out of row bound" + assert(col < Cols); // "out of column bound" + + return data[row * Cols + col]; +} + +template<typename T, std::size_t Rows, std::size_t Cols> +const T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) const { + assert(row < Rows); // "out of row bound" + assert(col < Cols); // "out of column bound" + + return data[row * Cols + col]; +} + +template<typename T, std::size_t Rows, std::size_t Cols> +auto mm::basic_matrix<T, Rows, Cols>::operator[](std::size_t index) { + if constexpr (is_row_vec() || is_col_vec()) { + return data.at(index); + } else { + return row_vec<T, Rows>( + data.cbegin() + (index * Cols), + data.cbegin() + ((index + 1) * Cols) + 1 + ); + } +} + +template<typename T, std::size_t Rows, std::size_t Cols> +void mm::basic_matrix<T, Rows, Cols>::swap_rows(std::size_t x, std::size_t y) { + if (x == y) + return; + + for (unsigned col = 0; col < Cols; col++) + std::swap(this->at(x, col), this->at(y, col)); +} + +template<typename T, std::size_t Rows, std::size_t Cols> +void mm::basic_matrix<T, Rows, Cols>::swap_cols(std::size_t x, std::size_t y) { + if (x == y) + return; + + for (unsigned row = 0; row < rows; row++) + std::swap(this->at(row, x), this->at(row, y)); +} + +template<typename T, std::size_t M, std::size_t N> +mm::basic_matrix<T, N, M> mm::basic_matrix<T, M, N>::transposed() const { + mm::basic_matrix<T, N, M> result; + + for (unsigned row = 0; row < M; row++) + for (unsigned col = 0; col < N; col++) + result.at(col, row) = this->at(row, col); + + return result; +} + + +/* operator overloading */ +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols> operator+( + const mm::basic_matrix<T, Rows, Cols>& a, + const mm::basic_matrix<T, Rows, Cols>& b +) { + mm::basic_matrix<T, Rows, Cols> result; + + for (unsigned row = 0; row < Rows; row++) + for (unsigned col = 0; col < Cols; col++) + result.at(row, col) = a.at(row, col) + b.at(row, col); + + return result; +} + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols> operator*( + const mm::basic_matrix<T, Rows, Cols>& m, + const T& scalar +) { + mm::basic_matrix<T, Rows, Cols> result; + for (unsigned row = 0; row < Rows; row++) + for (unsigned col = 0; col < Cols; col++) + result.at(row, col) = m.at(row, col) * scalar; + + return result; +} + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols> operator*( + const T& scalar, + const mm::basic_matrix<T, Rows, Cols>& m +) { + return m * scalar; +} + +template<typename T, std::size_t M, std::size_t P1, std::size_t P2, std::size_t N> +mm::basic_matrix<T, M, N> operator*( + const mm::basic_matrix<T, M, P1>& a, + const mm::basic_matrix<T, P2, N>& b +) { + static_assert(P1 == P2, "invalid matrix multiplication"); + mm::basic_matrix<T, M, N> result; + + // TODO: use a more efficient algorithm + for (unsigned row = 0; row < M; row++) + for (unsigned col = 0; col < N; col++) + for (unsigned k = 0; k < P1; k++) + result.at(row, col) = a.at(row, k) * b.at(k, col); + + return result; +} + + +template<typename T, std::size_t Rows, std::size_t Cols> +mm::basic_matrix<T, Rows, Cols> operator-( + const mm::basic_matrix<T, Rows, Cols>& a, + const mm::basic_matrix<T, Rows, Cols>& b +) { + return a + (static_cast<T>(-1) * b); +} + +template<typename T, std::size_t Rows, std::size_t Cols, unsigned NumW = 3> +std::ostream& operator<<(std::ostream& os, const mm::basic_matrix<T, Rows, Cols>& m) { + for (unsigned row = 0; row < Rows; row++) { + os << "[ "; + for (unsigned col = 0; col < (Cols -1); col++) { + os << std::setw(NumW) << m.at(row, col) << ", "; + } + os << std::setw(NumW) << m.at(row, (Cols -1)) << " ]\n"; + } + + return os; +} + + + + +/* row vector specialization */ +template<typename T, std::size_t Rows> +class mm::row_vec : public mm::basic_matrix<T, Rows, 1> { +public: + using mm::basic_matrix<T, Rows, 1>::basic_matrix; +}; + +/* column vector specialization */ +template<typename T, std::size_t Cols> +class mm::col_vec : public mm::basic_matrix<T, 1, Cols> { +public: + using mm::basic_matrix<T, 1, Cols>::basic_matrix; +}; + +/* general specialization (alias) */ +template<typename T, std::size_t Rows, std::size_t Cols> +class mm::matrix : public mm::basic_matrix<T, Rows, Cols> { +public: + using mm::basic_matrix<T, Rows, Cols>::basic_matrix; +}; + +/* square matrix specializaiton */ +template<typename T, std::size_t N> +class mm::square_matrix : public mm::basic_matrix<T, N, N> { +public: + using mm::basic_matrix<T, N, N>::basic_matrix; + + /// in place transpose + void transpose(); + inline void t() { transpose(); } + + T trace(); + inline T tr() { return trace(); } + + /// in place inverse + void invert(); + + // get the identity of size N + static inline constexpr square_matrix<T, N> identity() { + square_matrix<T, N> i; + for (unsigned row = 0; row < N; row++) + for (unsigned col = 0; col < N; col++) + i.at(row, col) = (row == col) ? 1 : 0; + + return i; + } +}; + + +template<typename T, std::size_t N> +void mm::square_matrix<T, N>::transpose() { + for (unsigned row = 0; row < N; row++) + for (unsigned col = 0; col < row; col++) + std::swap(this->at(row, col), this->at(col, row)); +} + +template<typename T, std::size_t N> +T mm::square_matrix<T, N>::trace() { + T sum = 0; + for (unsigned i = 0; i < N; i++) + sum += this->at(i, i); + + return sum; +} diff --git a/mmmatrix.cpp b/mmmatrix.cpp new file mode 100644 index 0000000..c0ce9eb --- /dev/null +++ b/mmmatrix.cpp @@ -0,0 +1 @@ +#include "mm/mmmatrix.hpp" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e80e737..89e7c47 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,3 +1,4 @@ +# vectors example add_executable(vec_example vec_example.cpp ) @@ -26,3 +27,34 @@ target_link_libraries(vec_example PRIVATE MM::mm ) + + +# matrix example +add_executable(matrix_example + matrix_example.cpp +) + +target_compile_options(matrix_example + PRIVATE + $<$<CXX_COMPILER_ID:GNU>: + -pedantic -Wall -Wextra -Wcast-qual -Wcast-align -Wpointer-arith + -Winit-self -Wshadow -Wswitch-enum -Wredundant-decls -Wfloat-equal + -Wundef -Wvla -Wconversion -Wstrict-aliasing + > + $<$<CXX_COMPILER_ID:MSVC>:/W4> +) + +target_compile_features(matrix_example + PRIVATE + cxx_std_17 +) + +target_include_directories(matrix_example + PRIVATE + ${MM_INCLUDE_DIRS} +) + +target_link_libraries(matrix_example + PRIVATE + MM::mm +) diff --git a/test/matrix_example.cpp b/test/matrix_example.cpp new file mode 100644 index 0000000..e7f6580 --- /dev/null +++ b/test/matrix_example.cpp @@ -0,0 +1,54 @@ +#include "mm/mmmatrix.hpp" + +#include <iostream> +#include <complex> + +int main(int argc, char *argv[]) { + std::cout << "MxN dimensional (int) matrices" << std::endl; + mm::matrix<int, 3, 2> a {{1, 2}, {3, 4}, {5, 6}}; + mm::matrix<int, 3, 2> b {{4, 3}, {9, 1}, {2, 5}}; + mm::matrix<int, 2, 4> c {{1, 2, 3, 4}, {5, 6, 7, 8}}; + + std::cout << "a = \n" << a; + std::cout << "b = \n" << b; + std::cout << "c = \n" << c; + std::cout << std::endl; + + // access elements + std::cout << "Access elements" << std::endl; + std::cout << "a.at(2,0) = " << a.at(2, 0) << std::endl; + std::cout << "a[2][0] = " << a[2][0] << std::endl;; + std::cout << std::endl; + + // basic operations + std::cout << "Basic operations" << std::endl; + std::cout << "a + b = \n" << a + b; + std::cout << "a - b = \n" << a - b; + std::cout << "a * c = \n" << a * c; + std::cout << "a * 2 = \n" << a * 2; + std::cout << "2 * a = \n" << 2 * a; + std::cout << "a.td() = \n" << a.td(); // or a.trasposed(); + std::cout << std::endl; + + // special matrices + mm::square_matrix<std::complex<int>, 2> f {{{2, 3}, {1, 4}}, {{6, 1}, {-3, 4}}}; + + std::cout << "Square matrix" << std::endl; + std::cout << "f = \n" << f; + + std::cout << "tr(f) = " << f.tr() /* or f.trace() */ << std::endl; + + f.t(); + std::cout << "after in place transpose f.t(), f = \n" << f; + std::cout << std::endl; + + + auto identity = mm::square_matrix<int, 3>::identity(); + + std::cout << "Identity matrix" << std::endl; + std::cout << "I = \n" << identity; + std::cout << std::endl; + + + return 0; +} |