matrix.h

Go to the documentation of this file.

 
#pragma once
 
#include <array>
#include <cassert>
#include <vector>
 
namespace qualpal {
 
template<typename T>
class Matrix
{
public:
  Matrix(std::size_t rows, std::size_t cols)
    : rows(rows)
    , cols(cols)
    , data(rows * cols)
  {
  }
 
  Matrix()
    : rows(0)
    , cols(0)
    , data()
  {
  }
 
  Matrix(std::size_t rows, std::size_t cols, const std::vector<T>& data)
    : rows(rows)
    , cols(cols)
    , data(data)
  {
    assert(rows * cols == data.size());
  }
 
  T& operator()(std::size_t row, std::size_t col)
  {
    return data[col * rows + row];
  }
 
  const T& operator()(std::size_t row, std::size_t col) const
  {
    return data[col * rows + row];
  }
 
  Matrix<T> transpose() const
  {
    Matrix<T> result(cols, rows);
    for (std::size_t i = 0; i < rows; ++i) {
      for (std::size_t j = 0; j < cols; ++j) {
        result(j, i) = (*this)(i, j);
      }
    }
    return result;
  }
 
  std::size_t ncol() const { return cols; }
 
  std::size_t nrow() const { return rows; }
 
private:
  std::size_t rows, cols;
  std::vector<T> data;
};
 
template<typename T, std::size_t rows, std::size_t cols>
class FixedMatrix
{
public:
  constexpr FixedMatrix()
    : data{}
  {
  }
 
  constexpr FixedMatrix(std::initializer_list<std::initializer_list<T>> list)
    : data{}
  {
    auto it = list.begin();
    for (std::size_t i = 0; i < rows && it != list.end(); ++i, ++it) {
      auto col_it = it->begin();
      for (std::size_t j = 0; j < cols && col_it != it->end(); ++j, ++col_it) {
        data[i * cols + j] = *col_it;
      }
    }
  }
 
  constexpr FixedMatrix<T, cols, rows> t() const
  {
    FixedMatrix<T, cols, rows> result;
    for (std::size_t i = 0; i < rows; ++i) {
      for (std::size_t j = 0; j < cols; ++j) {
        result(j, i) = (*this)(i, j);
      }
    }
    return result;
  }
 
  constexpr T& operator()(std::size_t row, std::size_t col)
  {
    return data[row * cols + col];
  }
 
  constexpr const T& operator()(std::size_t row, std::size_t col) const
  {
    return data[row * cols + col];
  }
 
  constexpr FixedMatrix<T, rows, cols> operator*(const T& scalar) const
  {
    FixedMatrix<T, rows, cols> result;
    for (std::size_t i = 0; i < rows * cols; ++i) {
      result.data[i] = this->data[i] * scalar;
    }
    return result;
  }
 
  constexpr std::array<T, rows> operator*(const std::array<T, cols>& vec) const
  {
    std::array<T, rows> result{};
    for (std::size_t i = 0; i < rows; ++i) {
      for (std::size_t j = 0; j < cols; ++j) {
        result[i] += (*this)(i, j) * vec[j];
      }
    }
    return result;
  }
 
  template<std::size_t other_cols>
  constexpr FixedMatrix<T, rows, other_cols> operator*(
    const FixedMatrix<T, cols, other_cols>& other) const
  {
    FixedMatrix<T, rows, other_cols> result;
    result.zeros();
    for (std::size_t i = 0; i < rows; ++i) {
      for (std::size_t j = 0; j < other_cols; ++j) {
        for (std::size_t k = 0; k < cols; ++k) {
          result(i, j) += (*this)(i, k) * other(k, j);
        }
      }
    }
    return result;
  }
 
  constexpr Matrix<T> operator*(const Matrix<T>& other) const
  {
    assert(cols == other.nrow());
    Matrix<T> result(rows, other.ncol());
    for (std::size_t i = 0; i < rows; ++i) {
      for (std::size_t j = 0; j < other.ncol(); ++j) {
        for (std::size_t k = 0; k < cols; ++k) {
          result(i, j) += (*this)(i, k) * other(k, j);
        }
      }
    }
    return result;
  }
 
  constexpr FixedMatrix<T, rows, cols> operator+(
    const FixedMatrix<T, rows, cols>& other) const
  {
    FixedMatrix<T, rows, cols> result;
    for (std::size_t i = 0; i < rows * cols; ++i) {
      result.data[i] = this->data[i] + other.data[i];
    }
    return result;
  }
 
  constexpr FixedMatrix<T, rows, cols> operator-(
    const FixedMatrix<T, rows, cols>& other) const
  {
    FixedMatrix<T, rows, cols> result;
    for (std::size_t i = 0; i < rows * cols; ++i) {
      result.data[i] = this->data[i] - other.data[i];
    }
    return result;
  }
 
  constexpr void fill(const T& value) { data.fill(value); }
 
  constexpr void zeros() { data.fill(0); }
 
private:
  std::array<T, rows * cols> data;
};
 
} // namespace qualpal