matrix.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_MATRIX_HPP
#define NUMERICAL_MATRIX_HPP
 
#include <algorithm>
#include <cmath>
#include <complex> // so std::abs resolves for a complex scalar, regardless of include order
#include <cstddef>
#include <functional>
#include <initializer_list>
#include <utility>
#include <vector>
 
#include "gemm.hpp"
#include "parallel.hpp"
#include "vector.hpp"
 
namespace numerical {
 
  template <typename T>
  class matrix
  {
  public:
 
    matrix() = default;
 
    matrix(std::initializer_list<std::initializer_list<T>> rows)
      : rows_(rows.size())
    {
      std::size_t r {0};
      for (const std::initializer_list<T>& row : rows) {
        if (r == 0) {
          cols_ = row.size();
        }
        else if (row.size() != cols_) {
          throw dimension_error("matrix rows have unequal lengths");
        }
        data_.insert(data_.end(), row.begin(), row.end());
        ++r;
      }
    }
 
    explicit matrix(std::size_t rows,
                    std::size_t cols,
                    buffer<T>   elements)
      : rows_(rows), cols_(cols), data_(std::move(elements))
    {
      if (data_.size() != detail::checked_mul(rows_, cols_)) {
        throw dimension_error("matrix buffer size does not match dimensions");
      }
    }
 
    explicit matrix(std::size_t           rows,
                    std::size_t           cols,
                    const std::vector<T>& elements)
      : matrix(rows,
               cols,
               buffer<T>(elements.begin(), elements.end()))
    {}
 
    explicit matrix(std::size_t              rows,
                    std::size_t              cols,
                    std::initializer_list<T> elements)
      : matrix(rows,
               cols,
               buffer<T>(elements))
    {}
 
    [[nodiscard]] std::size_t rows() const noexcept
    {
      return rows_;
    }
 
    [[nodiscard]] std::size_t cols() const noexcept
    {
      return cols_;
    }
 
    [[nodiscard]] const T& operator()(std::size_t row,
                                      std::size_t col) const
    {
      return data_[(row * cols_) + col];
    }
 
    [[nodiscard]] bool operator==(const matrix& other) const
    {
      return rows_ == other.rows_ && cols_ == other.cols_ && data_ == other.data_;
    }
 
    [[nodiscard]] typename buffer<T>::const_iterator begin() const noexcept
    {
      return data_.begin();
    }
 
    [[nodiscard]] typename buffer<T>::const_iterator end() const noexcept
    {
      return data_.end();
    }
 
    [[nodiscard]] buffer<T> elements() &&
    {
      rows_ = 0;
      cols_ = 0;
      return std::move(data_);
    }
 
  private:
 
    std::size_t    rows_ {0}; 
    std::size_t    cols_ {0}; 
    buffer<T>      data_;     
  };
 
  template <typename T>
  [[nodiscard]] matrix<T> operator+(const matrix<T>& a,
                                    const matrix<T>& b)
  {
    if (a.rows() != b.rows() || a.cols() != b.cols()) {
      throw dimension_error("matrix shapes differ in '+'");
    }
    buffer<T> out(a.rows() * a.cols());
    std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::plus<T> {});
    return matrix<T> {a.rows(), a.cols(), std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> operator-(const matrix<T>& a,
                                    const matrix<T>& b)
  {
    if (a.rows() != b.rows() || a.cols() != b.cols()) {
      throw dimension_error("matrix shapes differ in '-'");
    }
    buffer<T> out(a.rows() * a.cols());
    std::transform(a.begin(), a.end(), b.begin(), out.begin(), std::minus<T> {});
    return matrix<T> {a.rows(), a.cols(), std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> operator*(const T&         scalar,
                                    const matrix<T>& m)
  {
    buffer<T> out(m.rows() * m.cols());
    std::transform(m.begin(), m.end(), out.begin(),
                   [&scalar](const T& element) { return scalar * element; });
    return matrix<T> {m.rows(), m.cols(), std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> operator*(const matrix<T>& m,
                                    const T&         scalar)
  {
    return scalar * m;
  }
 
  template <typename T>
  [[nodiscard]] vector<T> operator*(const matrix<T>& a,
                                    const vector<T>& v)
  {
    if (a.cols() != v.size()) {
      throw dimension_error("matrix–vector product: inner dimensions differ");
    }
    buffer<T> out(a.rows());
    if (a.cols() == 0) {
      std::fill(out.begin(), out.end(), T {}); // each element an empty sum
    }
    else {
      detail::matrix_times_vector(&*a.begin(), a.rows(), a.cols(), &v[0], out.data());
    }
    return vector<T> {std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> operator*(const matrix<T>& a,
                                    const matrix<T>& b)
  {
    if (a.cols() != b.rows()) {
      throw dimension_error("matrix–matrix product: inner dimensions differ");
    }
    // detail::gemm adds to `out`, every element in ascending order of the inner index, whatever the blocking and
    // the number of threads (see gemm.hpp): `out` starts at zero, written by the threads that will compute it.
    buffer<T> out(detail::checked_mul(a.rows(), b.cols())); // guard the cross-product alloc
    detail::parallel_zero(out.data(), out.size());
    if (!out.empty() && a.cols() != 0) {
      detail::gemm(a.rows(), b.cols(), a.cols(), &*a.begin(), a.cols(), &*b.begin(), b.cols(), out.data(), b.cols());
    }
    return matrix<T> {a.rows(), b.cols(), std::move(out)};
  }
 
  inline constexpr std::size_t transpose_tile {32};
 
  template <typename T>
  [[nodiscard]] matrix<T> transpose(const matrix<T>& m)
  {
    const std::size_t rows {m.rows()};
    const std::size_t cols {m.cols()};
    buffer<T>         out(rows * cols);
    detail::parallel_for((cols + transpose_tile - 1) / transpose_tile, transpose_tile * rows, [&](std::size_t block) {
                           const std::size_t j0 {block * transpose_tile};
                           const std::size_t j1 {std::min(j0 + transpose_tile, cols)};
                           for (std::size_t i0 {0}; i0 < rows; i0 += transpose_tile) {
                             const std::size_t i1 {std::min(i0 + transpose_tile, rows)};
                             for (std::size_t j {j0}; j < j1; ++j) {
                               for (std::size_t i {i0}; i < i1; ++i) {
                                 out[(j * rows) + i] = m(i, j);
                               }
                             }
                           }
                         });
    return matrix<T> {cols, rows, std::move(out)};
  }
 
  template <typename T>
  [[nodiscard]] bool approx_equal(const matrix<T>& a,
                                  const matrix<T>& b,
                                  double           tolerance = 1e-9)
  {
    if (a.rows() != b.rows() || a.cols() != b.cols()) {
      return false;
    }
    for (std::size_t i {0}; i < a.rows(); ++i) {
      for (std::size_t j {0}; j < a.cols(); ++j) {
        // Equal, or within the tolerance: a NaN is neither, so it is never close, and two equal
        // infinities are equal although their difference is NaN.
        if (!(a(i, j) == b(i, j)) && !(std::abs(a(i, j) - b(i, j)) <= tolerance)) {
          return false;
        }
      }
    }
    return true;
  }
} // namespace numerical
 
#endif // NUMERICAL_MATRIX_HPP