Add basic ops

10 files changed, 217 insertions, 1 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index d38773b..1c854de 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,6 +13,7 @@ option(LINEAR_ALGEBRA_BUILD_EXAMPLES "Build example programs" ON)
 add_library(linear_algebra
     src/vector.cpp
     src/matrix.cpp
+    src/norms.cpp
 )
 
 add_library(linear_algebra::core ALIAS linear_algebra)
diff --git a/README.md b/README.md
index cb776c1..6f49640 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # Numerical Linear Algebra
 
-This repo contains a small C++ dense numerical linear algebra library for double, with a companion experiments directory for evaluating performance.
+This repo contains a small C++ dense numerical linear algebra library for `double`, with a companion experiments directory for evaluating performance.
 
 ## Build
 
diff --git a/include/matrix.hpp b/include/matrix.hpp
index ba3e4ae..d78f814 100644
--- a/include/matrix.hpp
+++ b/include/matrix.hpp
@@ -3,6 +3,7 @@
 #include <cstddef>
 #include <initializer_list>
 #include <vector>
+#include "vector.hpp"
 
 namespace linalg {
 
@@ -37,5 +38,11 @@ private:
     std::vector<double> data_;
 };
 
+Matrix transpose(const Matrix& matrix);
+Matrix operator+(const Matrix& lhs, const Matrix& rhs);
+Matrix operator-(const Matrix& lhs, const Matrix& rhs);
+Vector operator*(const Matrix& matrix, const Vector& vector);
+Matrix operator*(const Matrix& lhs, const Matrix& rhs);
+
 }  // namespace linalg
 
diff --git a/include/norms.hpp b/include/norms.hpp
new file mode 100644
index 0000000..85fad41
--- /dev/null
+++ b/include/norms.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "vector.hpp"
+
+namespace linalg {
+
+double norm2(const Vector& vector);
+
+}
diff --git a/include/vector.hpp b/include/vector.hpp
index e6bc0ce..3e797d8 100644
--- a/include/vector.hpp
+++ b/include/vector.hpp
@@ -41,6 +41,7 @@ Vector operator+(const Vector& lhs, const Vector& rhs);
 Vector operator-(const Vector& lhs, const Vector& rhs);
 Vector operator*(const Vector& v, double scalar);
 Vector operator*(double scalar, const Vector& v);
+Vector operator/(const Vector& v, double scalar);
 double dot(const Vector& lhs, const Vector& rhs);
 
 }  // namespace linalg
diff --git a/src/matrix.cpp b/src/matrix.cpp
index dd19e02..c6888a4 100644
--- a/src/matrix.cpp
+++ b/src/matrix.cpp
@@ -1,10 +1,24 @@
 #include "matrix.hpp"
+#include "linalg_error.hpp"
 
 #include <algorithm>
+#include <sstream>
 #include <stdexcept>
 
 namespace linalg {
 
+namespace {
+
+void check_same_shape(const Matrix& lhs, const Matrix& rhs, const char* operation) {
+    if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) {
+        std::ostringstream oss;
+        oss << operation << " requires equal matrix shapes, got " << lhs.rows() << "x" << lhs.cols() << " and " << rhs.rows() << "x" << rhs.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+}
+
+}  // namespace
+
 Matrix::Matrix(std::size_t rows, std::size_t cols)
     : rows_(rows), cols_(cols), data_(rows * cols) {}
 
@@ -69,4 +83,77 @@ void Matrix::check_bounds(std::size_t i, std::size_t j) const {
     }
 }
 
+Matrix transpose(const Matrix& matrix) {
+    Matrix result(matrix.cols(), matrix.rows());
+    for (std::size_t i = 0; i < matrix.rows(); ++i) {
+        for (std::size_t j = 0; j < matrix.cols(); ++j) {
+            result(j, i) = matrix(i, j);
+        }
+    }
+    return result;
+}
+
+Matrix operator+(const Matrix& lhs, const Matrix& rhs) {
+    check_same_shape(lhs, rhs, "Matrix addition");
+
+    Matrix result(lhs.rows(), lhs.cols());
+    for (std::size_t i = 0; i < lhs.rows(); ++i) {
+        for (std::size_t j = 0; j < lhs.cols(); ++j) {
+            result(i, j) = lhs(i, j) + rhs(i, j);
+        }
+    }
+    return result;
+}
+
+Matrix operator-(const Matrix& lhs, const Matrix& rhs) {
+    check_same_shape(lhs, rhs, "Matrix subtraction");
+
+    Matrix result(lhs.rows(), lhs.cols());
+    for (std::size_t i = 0; i < lhs.rows(); ++i) {
+        for (std::size_t j = 0; j < lhs.cols(); ++j) {
+            result(i, j) = lhs(i, j) - rhs(i, j);
+        }
+    }
+    return result;
+}
+
+Vector operator*(const Matrix& matrix, const Vector& vector) {
+    if (matrix.cols() != vector.size()) {
+        std::ostringstream oss;
+        oss << "Matrix-vector multiplication requires matrix columns to match vector size, got "
+            << matrix.cols() << " and " << vector.size();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    Vector result(matrix.rows());
+    for (std::size_t i = 0; i < matrix.rows(); ++i) {
+        double sum = 0.0;
+        for (std::size_t j = 0; j < matrix.cols(); ++j) {
+            sum += matrix(i, j) * vector[j];
+        }
+        result[i] = sum;
+    }
+    return result;
+}
+
+Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
+    if (lhs.cols() != rhs.rows()) {
+        std::ostringstream oss;
+        oss << "Matrix multiplication requires lhs.cols() == rhs.rows(), got " << lhs.cols()
+            << " and " << rhs.rows();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    Matrix result(lhs.rows(), rhs.cols());
+    for (std::size_t i = 0; i < lhs.rows(); ++i) {
+        for (std::size_t k = 0; k < lhs.cols(); ++k) {
+            const double lhs_ik = lhs(i, k);
+            for (std::size_t j = 0; j < rhs.cols(); ++j) {
+                result(i, j) += lhs_ik * rhs(k, j);
+            }
+        }
+    }
+    return result;
+}
+
 }  // namespace linalg
diff --git a/src/norms.cpp b/src/norms.cpp
new file mode 100644
index 0000000..cafee69
--- /dev/null
+++ b/src/norms.cpp
@@ -0,0 +1,9 @@
+#include "norms.hpp"
+
+#include <cmath>
+
+namespace linalg {
+
+double norm2(const Vector& vector) { return std::sqrt(dot(vector, vector)); }
+
+}  // namespace linalg
diff --git a/src/vector.cpp b/src/vector.cpp
index b8b1c6f..8c4a8c4 100644
--- a/src/vector.cpp
+++ b/src/vector.cpp
@@ -83,6 +83,18 @@ Vector operator*(const Vector& v, double scalar) {
 
 Vector operator*(double scalar, const Vector& v) { return v * scalar; }
 
+Vector operator/(const Vector& v, double scalar) {
+    if (scalar == 0.0) {
+        throw std::invalid_argument("Vector scalar division requires a nonzero scalar");
+    }
+
+    Vector result(v.size());
+    for (std::size_t i = 0; i < v.size(); ++i) {
+        result[i] = v[i] / scalar;
+    }
+    return result;
+}
+
 double dot(const Vector& lhs, const Vector& rhs) {
     check_same_size(lhs, rhs, "Dot product");
 
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
index b5fb60a..4cfc783 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -1,11 +1,15 @@
+#include "linalg_error.hpp"
 #include "matrix.hpp"
 
+#include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
 
 #include <type_traits>
 #include <utility>
 
 using linalg::Matrix;
+using linalg::Vector;
+using linalg::DimensionMismatchError;
 
 TEST_CASE("Matrix constructors initialize dimensions and values", "[matrix]") {
     const Matrix empty;
@@ -94,3 +98,79 @@ TEST_CASE("Matrix copy, move, and initializer list keeps row-major layout", "[ma
 TEST_CASE("Matrix initializer list rejects unequal row lengths", "[matrix]") {
     CHECK_THROWS_AS((Matrix{{1.0, 2.0}, {3.0}}), std::invalid_argument);
 }
+
+TEST_CASE("Matrix transpose swaps rows and columns", "[matrix]") {
+    const Matrix a{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}};
+    const Matrix at = linalg::transpose(a);
+
+    REQUIRE(at.rows() == 3);
+    REQUIRE(at.cols() == 2);
+    CHECK(at(0, 0) == 1.0);
+    CHECK(at(1, 0) == 2.0);
+    CHECK(at(2, 0) == 3.0);
+    CHECK(at(0, 1) == 4.0);
+    CHECK(at(1, 1) == 5.0);
+    CHECK(at(2, 1) == 6.0);
+}
+
+TEST_CASE("Matrix addition and subtraction enforce equal shapes", "[matrix]") {
+    const Matrix a{{1.0, 2.0}, {3.0, 4.0}};
+    const Matrix b{{0.5, -1.0}, {2.0, 1.5}};
+
+    const Matrix sum = a + b;
+    CHECK(sum(0, 0) == Catch::Approx(1.5));
+    CHECK(sum(0, 1) == Catch::Approx(1.0));
+    CHECK(sum(1, 0) == Catch::Approx(5.0));
+    CHECK(sum(1, 1) == Catch::Approx(5.5));
+
+    const Matrix diff = a - b;
+    CHECK(diff(0, 0) == Catch::Approx(0.5));
+    CHECK(diff(0, 1) == Catch::Approx(3.0));
+    CHECK(diff(1, 0) == Catch::Approx(1.0));
+    CHECK(diff(1, 1) == Catch::Approx(2.5));
+
+    const Matrix wrong_shape(3, 1);
+    CHECK_THROWS_AS(a + wrong_shape, DimensionMismatchError);
+    CHECK_THROWS_AS(a - wrong_shape, DimensionMismatchError);
+}
+
+TEST_CASE("Matrix-vector multiply checks dimensions", "[matrix]") {
+    const Matrix a{{1.0, 2.0, 3.0}, {0.0, -1.0, 4.0}};
+    const Vector x{2.0, -1.0, 0.5};
+
+    const Vector y = a * x;
+    REQUIRE(y.size() == 2);
+    CHECK(y[0] == Catch::Approx(1.5));
+    CHECK(y[1] == Catch::Approx(3.0));
+
+    const Vector wrong_size{1.0, 2.0};
+    CHECK_THROWS_AS(a * wrong_size, DimensionMismatchError);
+}
+
+TEST_CASE("Matrix-matrix multiply handles identity and shape checks", "[matrix]") {
+    const Matrix a{{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
+    const Matrix identity = Matrix::identity(2);
+    const Matrix product = a * identity;
+
+    REQUIRE(product.rows() == a.rows());
+    REQUIRE(product.cols() == a.cols());
+    for (std::size_t i = 0; i < a.rows(); ++i) {
+        for (std::size_t j = 0; j < a.cols(); ++j) {
+            CHECK(product(i, j) == Catch::Approx(a(i, j)));
+        }
+    }
+
+    const Matrix lhs{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}};
+    const Matrix rhs{{7.0, 8.0}, {9.0, 10.0}, {11.0, 12.0}};
+    const Matrix dense_product = lhs * rhs;
+
+    REQUIRE(dense_product.rows() == 2);
+    REQUIRE(dense_product.cols() == 2);
+    CHECK(dense_product(0, 0) == Catch::Approx(58.0));
+    CHECK(dense_product(0, 1) == Catch::Approx(64.0));
+    CHECK(dense_product(1, 0) == Catch::Approx(139.0));
+    CHECK(dense_product(1, 1) == Catch::Approx(154.0));
+
+    const Matrix incompatible(4, 1);
+    CHECK_THROWS_AS(lhs * incompatible, DimensionMismatchError);
+}
diff --git a/tests/test_vector.cpp b/tests/test_vector.cpp
index 1d7b91a..364ce31 100644
--- a/tests/test_vector.cpp
+++ b/tests/test_vector.cpp
@@ -1,6 +1,8 @@
 #include "linalg_error.hpp"
+#include "norms.hpp"
 #include "vector.hpp"
 
+#include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
 
 #include <type_traits>
@@ -86,3 +88,11 @@ TEST_CASE("Vector arithmetic enforces shape compatibility", "[vector]") {
     CHECK_THROWS_AS(a + short_vec, DimensionMismatchError);
     CHECK_THROWS_AS(linalg::dot(a, short_vec), DimensionMismatchError);
 }
+
+TEST_CASE("Vector 2-norm matches manually computed values", "[vector][norms]") {
+    const Vector v{3.0, 4.0};
+    CHECK(linalg::norm2(v) == Catch::Approx(5.0));
+
+    const Vector zero(5);
+    CHECK(linalg::norm2(zero) == Catch::Approx(0.0));
+}