FESADev/src/fesa/math/matrix.cpp

#include "fesa/math/matrix.h"

#include <mkl.h>

#include <limits>
#include <stdexcept>
#include <utility>

#include "math/dense_blas_internal.h"

namespace fesa {
namespace {

/// @brief Rejects shape overflow before logical dimensions diverge from
/// storage.
std::size_t CheckedStorageSize(const std::size_t rows,
                               const std::size_t columns) {
  if (columns != 0 &&
      rows > (std::numeric_limits<std::size_t>::max)() / columns) {
    throw std::length_error{
        "Dense matrix dimensions exceed the storage size range."};
  }
  return rows * columns;
}

/// @brief Preserves the Matrix exception contract around private size
/// conversion.
MKL_INT MklSizeOrThrow(const std::size_t size) {
  const auto converted = dense_blas_internal::ToMklSize(size);
  if (!converted.HasValue()) {
    throw std::length_error{
        "Dense matrix dimension exceeds the MKL integer range."};
  }
  return converted.Value();
}

}  // namespace

Matrix::Matrix(const std::size_t rows, const std::size_t columns,
               const double value)
    : rows_(rows),
      columns_(columns),
      values_(CheckedStorageSize(rows, columns), value) {}

Matrix::Matrix(const Matrix& other)
    : rows_(other.rows_),
      columns_(other.columns_),
      values_(other.values_.size()) {
  dense_blas_internal::CopyValues(other.values_.data(), other.values_.size(),
                                  values_.data());
}

Matrix::Matrix(Matrix&& other) noexcept
    : rows_(other.rows_),
      columns_(other.columns_),
      values_(std::move(other.values_)) {
  other.rows_ = 0;
  other.columns_ = 0;
  other.values_.clear();
}

Matrix& Matrix::operator=(const Matrix& other) {
  if (this != &other) {
    std::vector<double> copied(other.values_.size());
    dense_blas_internal::CopyValues(other.values_.data(), other.values_.size(),
                                    copied.data());
    rows_ = other.rows_;
    columns_ = other.columns_;
    values_.swap(copied);
  }
  return *this;
}

Matrix& Matrix::operator=(Matrix&& other) noexcept {
  if (this != &other) {
    rows_ = other.rows_;
    columns_ = other.columns_;
    values_ = std::move(other.values_);
    other.rows_ = 0;
    other.columns_ = 0;
    other.values_.clear();
  }
  return *this;
}

std::size_t Matrix::Rows() const noexcept { return rows_; }

std::size_t Matrix::Columns() const noexcept { return columns_; }

double& Matrix::operator()(const std::size_t row, const std::size_t column) {
  if (row >= rows_ || column >= columns_) {
    throw std::out_of_range{"Matrix index is outside its dimensions."};
  }
  return values_[row * columns_ + column];
}

const double& Matrix::operator()(const std::size_t row,
                                 const std::size_t column) const {
  if (row >= rows_ || column >= columns_) {
    throw std::out_of_range{"Matrix index is outside its dimensions."};
  }
  return values_[row * columns_ + column];
}

Vector Matrix::Multiply(const Vector& rhs) const {
  if (columns_ != rhs.Size()) {
    throw std::invalid_argument{
        "Matrix-vector multiplication has incompatible dimensions."};
  }

  Vector result{rows_};
  if (rows_ == 0 || columns_ == 0) {
    return result;
  }

  // The owned layout is row-major, so the leading dimension is the column
  // count for the adapter call and remains invisible to public consumers.
  cblas_dgemv(CblasRowMajor, CblasNoTrans, MklSizeOrThrow(rows_),
              MklSizeOrThrow(columns_), 1.0, values_.data(),
              MklSizeOrThrow(columns_), rhs.Data(), 1, 0.0, result.Data(), 1);
  return result;
}

Matrix Matrix::Multiply(const Matrix& rhs) const {
  if (columns_ != rhs.rows_) {
    throw std::invalid_argument{
        "Matrix multiplication has incompatible dimensions."};
  }

  Matrix result{rows_, rhs.columns_};
  if (rows_ == 0 || columns_ == 0 || rhs.columns_ == 0) {
    return result;
  }

  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, MklSizeOrThrow(rows_),
              MklSizeOrThrow(rhs.columns_), MklSizeOrThrow(columns_), 1.0,
              values_.data(), MklSizeOrThrow(columns_), rhs.values_.data(),
              MklSizeOrThrow(rhs.columns_), 0.0, result.values_.data(),
              MklSizeOrThrow(rhs.columns_));
  return result;
}

}  // namespace fesa