Matrix.hpp

#pragma once
#include "qip/Vector.hpp" // for std::vector overloads
#include <array>
#include <cassert>
#include <complex>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <iostream>
#include <type_traits>
#include <utility>
#include <vector>
 
template <typename T>
struct is_complex : std::false_type {};
template <typename T>
struct is_complex<std::complex<T>> : std::true_type {};
template <typename T>
constexpr bool is_complex_v = is_complex<T>::value;
 
//==============================================================================
namespace LinAlg {
 
template <typename T>
class View;
 
//==============================================================================
template <typename T = double>
class Matrix {
 
protected:
  std::size_t m_rows;
  std::size_t m_cols;
  std::vector<T> m_data{};
 
public:
  Matrix() : m_rows(0), m_cols(0) {}
 
  Matrix(std::size_t rows, std::size_t cols)
      : m_rows(rows), m_cols(cols), m_data(rows * cols) {}
 
  // excplicit, since don't alow flot->int converions
  explicit Matrix(std::size_t dimension)
      : m_rows(dimension), m_cols(dimension), m_data(dimension * dimension) {}
 
  Matrix(std::initializer_list<std::initializer_list<T>> ll)
      : m_rows(ll.size()), m_cols(ll.begin()->size()), m_data{} {
    // way to avoid copy?
    m_data.reserve(m_rows * m_cols);
    for (auto &l : ll) {
      m_data.insert(m_data.end(), l.begin(), l.end());
    }
  }
 
  Matrix(std::size_t rows, std::size_t cols, std::initializer_list<T> l)
      : m_rows(rows), m_cols(cols), m_data{l} {
    assert(m_data.size() == rows * cols &&
           "initializer_list must be rows*cols");
  }
 
  Matrix(std::size_t rows, std::size_t cols, std::vector<T> &&v)
      : m_rows(rows), m_cols(cols), m_data{std::forward<std::vector<T>>(v)} {
    assert(m_data.size() == rows * cols &&
           "initializer_list must be rows*cols");
  }
  Matrix(std::size_t rows, std::size_t cols, const std::vector<T> &v)
      : m_rows(rows), m_cols(cols), m_data{v} {
    assert(m_data.size() == rows * cols &&
           "initializer_list must be rows*cols");
  }
 
  //============================================================================
 
  std::size_t rows() const { return m_rows; }
  std::size_t cols() const { return m_cols; }
  std::size_t size() const { return m_data.size(); }
 
  T *data() { return m_data.data(); }
  const T *data() const { return m_data.data(); }
 
  //============================================================================
 
  const T *operator[](std::size_t i) const { return &(m_data[i * m_cols]); }
  T *operator[](std::size_t i) { return &(m_data[i * m_cols]); }
 
  T &at(std::size_t row_i, std::size_t col_j) {
    assert(row_i < m_rows && col_j < m_cols);
    return m_data[row_i * m_cols + col_j];
  }
  T at(std::size_t row_i, std::size_t col_j) const {
    assert(row_i < m_rows && col_j < m_cols);
    return m_data[row_i * m_cols + col_j];
  }
  T &operator()(std::size_t i, std::size_t j) { return at(i, j); }
  T operator()(std::size_t i, std::size_t j) const { return at(i, j); }
 
  //============================================================================
 
  auto begin() { return m_data.begin(); }
  auto cbegin() const { return m_data.cbegin(); }
  auto end() { return m_data.end(); }
  auto cend() const { return m_data.cend(); }
 
  // [[deprecated]]
  const T *row(std::size_t row) const {
    return m_data.data() + long(row * m_cols);
  }
 
  [[nodiscard]] View<T> row_view(std::size_t row) {
    return View<T>(this->data(), row * m_cols, m_cols, 1ul);
  }
  [[nodiscard]] View<const T> row_view(std::size_t row) const {
    return View<const T>(this->data(), row * m_cols, m_cols, 1ul);
  }
  [[nodiscard]] View<T> column_view(std::size_t col) {
    return View<T>(this->data(), col, m_rows, m_rows);
  }
  [[nodiscard]] View<const T> column_view(std::size_t col) const {
    return View<const T>(this->data(), col, m_rows, m_rows);
  }
 
  //============================================================================
  [[nodiscard]] auto as_gsl_view();
 
  [[nodiscard]] auto as_gsl_view() const;
 
  //============================================================================
  // Basic matrix operations:
  //============================================================================
  //============================================================================
 
  [[nodiscard]] T determinant() const;
 
  Matrix<T> &invert_in_place();
 
  [[nodiscard]] Matrix<T> inverse() const;
  [[nodiscard]] Matrix<T> transpose() const;
 
  //============================================================================
 
  Matrix<T> &make_identity();
  Matrix<T> &zero();
  // //! M -> M + aI, for I=identity (adds a to diag elements), in place
  // Matrix<T> &plusIdent(T a = T(1));
  //============================================================================
 
  [[nodiscard]] Matrix<T> conj() const;
  [[nodiscard]] auto real() const;
  [[nodiscard]] auto imag() const;
  [[nodiscard]] auto complex() const;
 
  Matrix<T> &conj_in_place();
 
  //============================================================================
  Matrix<T> &mult_elements_by(const Matrix<T> &a);
 
  [[nodiscard]] friend Matrix<T> mult_elements(Matrix<T> a,
                                               const Matrix<T> &b) {
    return a.mult_elements_by(b);
  }
 
  //============================================================================
  // Operator overloads: +,-, scalar */
  Matrix<T> &operator+=(const Matrix<T> &rhs);
  Matrix<T> &operator-=(const Matrix<T> &rhs);
  Matrix<T> &operator*=(const T x);
  Matrix<T> &operator/=(const T x);
 
  //============================================================================
  // nb: these are defined inline here to avoid ambiguous overload?
  [[nodiscard]] friend Matrix<T> operator+(Matrix<T> lhs,
                                           const Matrix<T> &rhs) {
    return (lhs += rhs);
  }
  [[nodiscard]] friend Matrix<T> operator-(Matrix<T> lhs,
                                           const Matrix<T> &rhs) {
    return (lhs -= rhs);
  }
  [[nodiscard]] friend Matrix<T> operator*(const T x, Matrix<T> rhs) {
    return (rhs *= x);
  }
  [[nodiscard]] friend Matrix<T> operator*(Matrix<T> lhs, const T x) {
    return (lhs *= x);
  }
  [[nodiscard]] friend Matrix<T> operator/(Matrix<T> lhs, const T x) {
    return (lhs /= x);
  }
 
  //============================================================================
  Matrix<T> &operator+=(T aI);
  Matrix<T> &operator-=(T aI);
 
  [[nodiscard]] friend Matrix<T> operator+(Matrix<T> M, T aI) {
    return (M += aI);
  }
  [[nodiscard]] friend Matrix<T> operator-(Matrix<T> M, T aI) {
    return (M -= aI);
  }
 
  //============================================================================
  template <typename U>
  friend Matrix<U> operator*(const Matrix<U> &a, const Matrix<U> &b);
 
  //============================================================================
  template <typename U>
  friend std::ostream &operator<<(std::ostream &os, const Matrix<U> &a);
};
 
//==============================================================================
//==============================================================================
//==============================================================================
 
//==============================================================================
//==============================================================================
//==============================================================================
 
template <typename T>
constexpr auto myEps();
template <typename T>
bool equal(const Matrix<T> &lhs, const Matrix<T> &rhs, T eps = myEps<T>());
 
} // namespace LinAlg
 
#include "Matrix.ipp"