Vector.hpp

#pragma once
#include "Matrix.hpp"
 
namespace LinAlg {
 
//==============================================================================
template <typename T = double>
class Vector : public Matrix<T> {
public:
  Vector() : Matrix<T>() {}
 
  Vector(std::size_t dimension) : Matrix<T>(dimension, 1) {}
 
  Vector(std::initializer_list<T> l) : Matrix<T>(l.size(), 1, l) {}
 
  Vector(std::vector<T> &&v)
      : Matrix<T>(v.size(), 1, std::forward<std::vector<T>>(v)) {}
 
  Vector(const std::vector<T> &v) : Matrix<T>(v.size(), 1, v) {}
 
  Vector(const Matrix<T> &&m) : Matrix<T>(std::move(m)) {
    assert(m.cols() == 1 && "Can only convert Matrix to Vector if matrix has 1 "
                            "column. Traspose first?");
  }
 
  Vector(const Matrix<T> &m) : Matrix<T>(m) {
    assert(m.cols() == 1 && "Can only convert Matrix to Vector if matrix has 1 "
                            "column. Traspose first?");
  }
 
  //============================================================================
 
  T operator[](std::size_t i) const { return this->data()[i]; }
  T &operator[](std::size_t i) { return this->data()[i]; }
  T &at(std::size_t i) {
    assert(i < this->size());
    return this->data()[i];
  }
  T at(std::size_t i) const {
    assert(i < this->size());
    return this->data()[i];
  }
  T &operator()(std::size_t i) { return at(i); }
  T operator()(std::size_t i) const { return at(i); }
 
  //============================================================================
  [[nodiscard]] Vector<T> conj() const;
  [[nodiscard]] auto real() const;
  [[nodiscard]] auto imag() const;
  [[nodiscard]] auto complex() const;
 
  Vector<T> transpose() const = delete;
 
  //============================================================================
  Vector<T> &operator+=(const Vector<T> &rhs);
  Vector<T> &operator-=(const Vector<T> rhs);
  Vector<T> &operator*=(const T x);
  Vector<T> &operator/=(const T x);
 
  [[nodiscard]] friend Vector<T> operator+(Vector<T> lhs,
                                           const Matrix<T> &rhs) {
    return (lhs += rhs);
  }
  [[nodiscard]] friend Vector<T> operator-(Vector<T> lhs,
                                           const Matrix<T> &rhs) {
    return (lhs -= rhs);
  }
  [[nodiscard]] friend Vector<T> operator*(const T x, Vector<T> rhs) {
    return (rhs *= x);
  }
  [[nodiscard]] friend Vector<T> operator*(Vector<T> lhs, const T x) {
    return (lhs *= x);
  }
  [[nodiscard]] friend Vector<T> operator/(Vector<T> lhs, const T x) {
    return (lhs /= x);
  }
 
  //============================================================================
  [[nodiscard]] friend Vector<T> operator*(const Matrix<T> &a,
                                           const Vector<T> &b) {
    // https://www.gnu.org/software/gsl/doc/html/blas.html
    assert(a.cols() == b.rows());
    Vector<T> product(b.rows());
    const auto a_gsl = a.as_gsl_view();
    const auto b_gsl = b.as_gsl_view();
    auto product_gsl = product.as_gsl_view();
    if constexpr (std::is_same_v<T, double>) {
      gsl_blas_dgemv(CblasNoTrans, 1.0, &a_gsl.matrix, &b_gsl.vector, 0.0,
                     &product_gsl.vector);
    } else if constexpr (std::is_same_v<T, float>) {
      gsl_blas_sgemv(CblasNoTrans, 1.0f, &a_gsl.matrix, &b_gsl.vector, 0.0f,
                     &product_gsl.vector);
    } else if constexpr (std::is_same_v<T, std::complex<double>>) {
      gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, &a_gsl.matrix,
                     &b_gsl.vector, GSL_COMPLEX_ZERO, &product_gsl.vector);
    } else if constexpr (std::is_same_v<T, std::complex<float>>) {
      const gsl_complex_float one{1.0f, 0.0f};
      const gsl_complex_float zero{0.0f, 0.0f};
      gsl_blas_cgemv(CblasNoTrans, one, &a_gsl.matrix, &b_gsl.vector, zero,
                     &product_gsl.vector);
    }
 
    return product;
  }
 
  [[nodiscard]] friend T operator*(const Vector<T> &a, const Vector<T> &b) {
    // https://www.gnu.org/software/gsl/doc/html/blas.html
    assert(a.rows() == b.rows());
    const auto a_gsl = a.as_gsl_view();
    const auto b_gsl = b.as_gsl_view();
    T product = T(0);
    if constexpr (std::is_same_v<T, double>) {
      gsl_blas_ddot(&a_gsl.vector, &b_gsl.vector, &product);
    } else if constexpr (std::is_same_v<T, float>) {
      gsl_blas_sdot(&a_gsl.vector, &b_gsl.vector, &product);
    } else if constexpr (std::is_same_v<T, std::complex<double>>) {
      gsl_blas_zdotu(&a_gsl.vector, &b_gsl.vector,
                     reinterpret_cast<gsl_complex *>(&product));
    } else if constexpr (std::is_same_v<T, std::complex<float>>) {
      gsl_blas_cdotu(&a_gsl.vector, &b_gsl.vector,
                     reinterpret_cast<gsl_complex_float *>(&product));
    }
 
    return product;
  }
 
  [[nodiscard]] friend Matrix<T> outer_product(const Vector<T> &a,
                                               const Vector<T> &b) {
    Matrix<T> op(a.rows(), b.rows());
    for (std::size_t i = 0; i < op.rows(); ++i) {
      for (std::size_t j = 0; j < op.cols(); ++j) {
        op[i][j] = a[i] * b[j];
      }
    }
    return op;
  }
 
  //============================================================================
  [[nodiscard]] auto as_gsl_view();
 
  [[nodiscard]] auto as_gsl_view() const;
};
 
} // namespace LinAlg
 
#include "Vector.ipp"