2 files changed, 164 insertions, 0 deletions
diff --git a/src/matrix.cpp b/src/matrix.cpp
new file mode 100644
index 0000000..dd19e02
--- /dev/null
+++ b/src/matrix.cpp
@@ -0,0 +1,72 @@
+#include "matrix.hpp"
+
+#include <algorithm>
+#include <stdexcept>
+
+namespace linalg {
+
+Matrix::Matrix(std::size_t rows, std::size_t cols)
+    : rows_(rows), cols_(cols), data_(rows * cols) {}
+
+Matrix::Matrix(std::size_t rows, std::size_t cols, double value)
+    : rows_(rows), cols_(cols), data_(rows * cols, value) {}
+
+Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) : rows_(values.size()) {
+    if (rows_ == 0) {
+        cols_ = 0;
+        return;
+    }
+
+    cols_ = values.begin()->size();
+    data_.reserve(rows_ * cols_);
+
+    for (const auto& row : values) {
+        if (row.size() != cols_) {
+            throw std::invalid_argument("Matrix initializer rows must have equal length");
+        }
+        data_.insert(data_.end(), row.begin(), row.end());
+    }
+}
+
+
+std::size_t Matrix::rows() const noexcept { return rows_; }
+
+std::size_t Matrix::cols() const noexcept { return cols_; }
+
+bool Matrix::empty() const noexcept { return data_.empty(); }
+
+double& Matrix::operator()(std::size_t i, std::size_t j) {
+    check_bounds(i, j);
+    return data_[index(i, j)];
+}
+
+const double& Matrix::operator()(std::size_t i, std::size_t j) const {
+    check_bounds(i, j);
+    return data_[index(i, j)];
+}
+
+void Matrix::fill(double value) { std::fill(data_.begin(), data_.end(), value); }
+
+double* Matrix::data() noexcept { return data_.data(); }
+
+const double* Matrix::data() const noexcept { return data_.data(); }
+
+Matrix Matrix::identity(std::size_t n) {
+    Matrix result(n, n);
+    for (std::size_t i = 0; i < n; ++i) {
+        result(i, i) = 1.0;
+    }
+    return result;
+}
+
+Matrix Matrix::zeros(std::size_t rows, std::size_t cols) { return Matrix(rows, cols, 0.0); }
+
+std::size_t Matrix::index(std::size_t i, std::size_t j) const { return i * cols_ + j; }
+
+void Matrix::check_bounds(std::size_t i, std::size_t j) const {
+    if (i >= rows_ || j >= cols_) {
+        throw std::out_of_range("Matrix index out of bounds");
+    }
+}
+
+}  // namespace linalg
diff --git a/src/vector.cpp b/src/vector.cpp
new file mode 100644
index 0000000..b8b1c6f
--- /dev/null
+++ b/src/vector.cpp
@@ -0,0 +1,92 @@
+#include "vector.hpp"
+#include "linalg_error.hpp"
+
+#include <algorithm>
+#include <numeric>
+#include <sstream>
+#include <stdexcept>
+
+namespace linalg {
+
+namespace {
+
+void check_same_size(const Vector& lhs, const Vector& rhs, const char* operation) {
+    if (lhs.size() != rhs.size()) {
+        std::ostringstream oss;
+        oss << operation << " requires equal vector sizes, got " << lhs.size() << " and "
+            << rhs.size();
+        throw DimensionMismatchError(oss.str());
+    }
+}
+
+}  // namespace
+
+Vector::Vector(std::size_t n) : data_(n) {}
+
+Vector::Vector(std::size_t n, double value) : data_(n, value) {}
+
+Vector::Vector(std::initializer_list<double> values) : data_(values) {}
+
+std::size_t Vector::size() const noexcept { return data_.size(); }
+
+bool Vector::empty() const noexcept { return data_.empty(); }
+
+double& Vector::operator[](std::size_t i) {
+    check_index(i);
+    return data_[i];
+}
+
+const double& Vector::operator[](std::size_t i) const {
+    check_index(i);
+    return data_[i];
+}
+
+void Vector::fill(double value) { std::fill(data_.begin(), data_.end(), value); }
+
+double* Vector::data() noexcept { return data_.data(); }
+
+const double* Vector::data() const noexcept { return data_.data(); }
+
+void Vector::check_index(std::size_t i) const {
+    if (i >= data_.size()) {
+        throw std::out_of_range("Vector index out of range");
+    }
+}
+
+Vector operator+(const Vector& lhs, const Vector& rhs) {
+    check_same_size(lhs, rhs, "Vector addition");
+
+    Vector result(lhs.size());
+    for (std::size_t i = 0; i < lhs.size(); ++i) {
+        result[i] = lhs[i] + rhs[i];
+    }
+    return result;
+}
+
+Vector operator-(const Vector& lhs, const Vector& rhs) {
+    check_same_size(lhs, rhs, "Vector subtraction");
+
+    Vector result(lhs.size());
+    for (std::size_t i = 0; i < lhs.size(); ++i) {
+        result[i] = lhs[i] - rhs[i];
+    }
+    return result;
+}
+
+Vector operator*(const Vector& v, double scalar) {
+    Vector result(v.size());
+    for (std::size_t i = 0; i < v.size(); ++i) {
+        result[i] = v[i] * scalar;
+    }
+    return result;
+}
+
+Vector operator*(double scalar, const Vector& v) { return v * scalar; }
+
+double dot(const Vector& lhs, const Vector& rhs) {
+    check_same_size(lhs, rhs, "Dot product");
+
+    return std::inner_product(lhs.begin(), lhs.end(), rhs.begin(), 0.0);
+}
+
+}  // namespace linalg