feat(linear-static-3d-euler-beam): step 9 - dense-math-adapters

src/fesa/math/matrix.cpp

@@ -0,0 +1,150 @@
+#include "fesa/math/matrix.hpp"
+
+#include <mkl.h>
+
+#include <limits>
+#include <stdexcept>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+MKL_INT toMklSize(const std::size_t size) {
+    if (size > static_cast<std::size_t>((std::numeric_limits<MKL_INT>::max)())) {
+        throw std::length_error{"Dense matrix dimension exceeds the MKL integer range."};
+    }
+    return static_cast<MKL_INT>(size);
+}
+
+void copyValues(const std::vector<double>& source, std::vector<double>& destination) {
+    if (source.empty()) {
+        return;
+    }
+
+    cblas_dcopy(toMklSize(source.size()), source.data(), 1, destination.data(), 1);
+}
+
+} // namespace
+
+Matrix::Matrix(
+    const std::size_t rows,
+    const std::size_t columns,
+    const double value)
+    : rows_(rows), columns_(columns), values_(rows * columns, value) {}
+
+Matrix::Matrix(const Matrix& other)
+    : rows_(other.rows_), columns_(other.columns_), values_(other.values_.size()) {
+    copyValues(other.values_, values_);
+}
+
+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());
+        copyValues(other.values_, copied);
+        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,
+        toMklSize(rows_),
+        toMklSize(columns_),
+        1.0,
+        values_.data(),
+        toMklSize(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,
+        toMklSize(rows_),
+        toMklSize(rhs.columns_),
+        toMklSize(columns_),
+        1.0,
+        values_.data(),
+        toMklSize(columns_),
+        rhs.values_.data(),
+        toMklSize(rhs.columns_),
+        0.0,
+        result.values_.data(),
+        toMklSize(rhs.columns_));
+    return result;
+}
+
+} // namespace fesa