diff options
| -rw-r--r-- | CMakeLists.txt | 1 | ||||
| -rw-r--r-- | README.md | 2 | ||||
| -rw-r--r-- | include/matrix.hpp | 7 | ||||
| -rw-r--r-- | include/norms.hpp | 9 | ||||
| -rw-r--r-- | include/vector.hpp | 1 | ||||
| -rw-r--r-- | src/matrix.cpp | 87 | ||||
| -rw-r--r-- | src/norms.cpp | 9 | ||||
| -rw-r--r-- | src/vector.cpp | 12 | ||||
| -rw-r--r-- | tests/test_matrix.cpp | 80 | ||||
| -rw-r--r-- | tests/test_vector.cpp | 10 |
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) @@ -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)); +} |