aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.clang-format7
-rw-r--r--.gitignore6
-rw-r--r--CMakeLists.txt64
-rw-r--r--README.md24
-rw-r--r--include/linalg_error.hpp19
-rw-r--r--include/matrix.hpp41
-rw-r--r--include/vector.hpp46
-rw-r--r--src/matrix.cpp72
-rw-r--r--src/vector.cpp92
-rw-r--r--tests/test_matrix.cpp96
-rw-r--r--tests/test_vector.cpp88
11 files changed, 550 insertions, 5 deletions
diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000..02cc57f
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,7 @@
+BasedOnStyle: LLVM
+IndentWidth: 4
+ColumnLimit: 100
+BreakBeforeBraces: Attach
+AllowShortFunctionsOnASingleLine: Empty
+PointerAlignment: Left
+SortIncludes: true
diff --git a/.gitignore b/.gitignore
index 0180838..c8fac71 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,3 @@
-.zig-cache/
-zig-out/
-*.o
+build/
+.vscode/
+.DS_Store \ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..d38773b
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,64 @@
+cmake_minimum_required(VERSION 3.20)
+
+project(linear_algebra_cpp VERSION 0.1.0 LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+option(LINEAR_ALGEBRA_BUILD_TESTS "Build unit tests" ON)
+option(LINEAR_ALGEBRA_BUILD_EXAMPLES "Build example programs" ON)
+
+add_library(linear_algebra
+ src/vector.cpp
+ src/matrix.cpp
+)
+
+add_library(linear_algebra::core ALIAS linear_algebra)
+
+target_include_directories(linear_algebra
+ PUBLIC
+ ${PROJECT_SOURCE_DIR}/include
+)
+
+target_compile_features(linear_algebra PUBLIC cxx_std_20)
+
+if(CMAKE_CXX_COMPILER_ID MATCHES "Clang|GNU")
+ target_compile_options(linear_algebra PRIVATE -Wall -Wextra -Wpedantic -Wconversion)
+elseif(MSVC)
+ target_compile_options(linear_algebra PRIVATE /W4 /permissive-)
+endif()
+
+if(LINEAR_ALGEBRA_BUILD_EXAMPLES)
+ add_executable(solve_linear_system examples/solve_linear_system.cpp)
+ target_link_libraries(solve_linear_system PRIVATE linear_algebra::core)
+endif()
+
+if(LINEAR_ALGEBRA_BUILD_TESTS)
+ include(FetchContent)
+
+ FetchContent_Declare(
+ Catch2
+ GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+ GIT_TAG v3.5.4
+ )
+
+ FetchContent_MakeAvailable(Catch2)
+
+ enable_testing()
+
+ add_executable(linear_algebra_tests
+ tests/test_vector.cpp
+ tests/test_matrix.cpp
+ )
+
+ target_link_libraries(linear_algebra_tests
+ PRIVATE
+ linear_algebra::core
+ Catch2::Catch2WithMain
+ )
+
+ include(Catch)
+ catch_discover_tests(linear_algebra_tests)
+endif()
diff --git a/README.md b/README.md
index 5622524..cb776c1 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,22 @@
-# zig-raytracer
-Lightweight raytracer in Zig
+# Numerical Linear Algebra
+
+This repo contains a small C++ dense numerical linear algebra library for double, with a companion experiments directory for evaluating performance.
+
+## Build
+
+```bash
+cmake -S . -B build
+cmake --build build
+```
+
+## Run tests
+
+```bash
+ctest --test-dir build --output-on-failure
+```
+
+## Run the example
+
+```bash
+./build/solve_linear_system
+``` \ No newline at end of file
diff --git a/include/linalg_error.hpp b/include/linalg_error.hpp
new file mode 100644
index 0000000..8258899
--- /dev/null
+++ b/include/linalg_error.hpp
@@ -0,0 +1,19 @@
+#pragma once
+
+#include <stdexcept>
+#include <string>
+
+namespace linalg {
+
+class LinAlgError : public std::runtime_error {
+public:
+ using std::runtime_error::runtime_error;
+};
+
+class DimensionMismatchError : public LinAlgError {
+public:
+ explicit DimensionMismatchError(const std::string& message)
+ : LinAlgError(message) {}
+};
+
+} // namespace linalg
diff --git a/include/matrix.hpp b/include/matrix.hpp
new file mode 100644
index 0000000..ba3e4ae
--- /dev/null
+++ b/include/matrix.hpp
@@ -0,0 +1,41 @@
+#pragma once
+
+#include <cstddef>
+#include <initializer_list>
+#include <vector>
+
+namespace linalg {
+
+class Matrix {
+public:
+ Matrix() = default;
+ Matrix(std::size_t rows, std::size_t cols);
+ Matrix(std::size_t rows, std::size_t cols, double value);
+ Matrix(std::initializer_list<std::initializer_list<double>> values);
+
+ [[nodiscard]] std::size_t rows() const noexcept;
+ [[nodiscard]] std::size_t cols() const noexcept;
+ [[nodiscard]] bool empty() const noexcept;
+
+ double& operator()(std::size_t i, std::size_t j);
+ const double& operator()(std::size_t i, std::size_t j) const;
+
+ void fill(double value);
+
+ double* data() noexcept;
+ const double* data() const noexcept;
+
+ static Matrix identity(std::size_t n);
+ static Matrix zeros(std::size_t rows, std::size_t cols);
+
+private:
+ [[nodiscard]] std::size_t index(std::size_t i, std::size_t j) const;
+ void check_bounds(std::size_t i, std::size_t j) const;
+
+ std::size_t rows_ = 0;
+ std::size_t cols_ = 0;
+ std::vector<double> data_;
+};
+
+} // namespace linalg
+
diff --git a/include/vector.hpp b/include/vector.hpp
new file mode 100644
index 0000000..e6bc0ce
--- /dev/null
+++ b/include/vector.hpp
@@ -0,0 +1,46 @@
+#pragma once
+
+#include <cstddef>
+#include <initializer_list>
+#include <vector>
+
+namespace linalg {
+
+class Vector {
+public:
+ Vector() = default;
+ explicit Vector(std::size_t n);
+ Vector(std::size_t n, double value);
+ Vector(std::initializer_list<double> values);
+
+ [[nodiscard]] std::size_t size() const noexcept;
+ [[nodiscard]] bool empty() const noexcept;
+
+ double& operator[](std::size_t i);
+ const double& operator[](std::size_t i) const;
+
+ void fill(double value);
+
+ double* data() noexcept;
+ const double* data() const noexcept;
+
+ auto begin() noexcept { return data_.begin(); }
+ auto end() noexcept { return data_.end(); }
+ auto begin() const noexcept { return data_.begin(); }
+ auto end() const noexcept { return data_.end(); }
+ auto cbegin() const noexcept { return data_.cbegin(); }
+ auto cend() const noexcept { return data_.cend(); }
+
+private:
+ void check_index(std::size_t i) const;
+
+ std::vector<double> data_;
+};
+
+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);
+double dot(const Vector& lhs, const Vector& rhs);
+
+} // namespace linalg
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
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
new file mode 100644
index 0000000..b5fb60a
--- /dev/null
+++ b/tests/test_matrix.cpp
@@ -0,0 +1,96 @@
+#include "matrix.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <type_traits>
+#include <utility>
+
+using linalg::Matrix;
+
+TEST_CASE("Matrix constructors initialize dimensions and values", "[matrix]") {
+ const Matrix empty;
+ CHECK(empty.rows() == 0);
+ CHECK(empty.cols() == 0);
+ CHECK(empty.empty());
+
+ const Matrix zeroed(2, 3);
+ REQUIRE(zeroed.rows() == 2);
+ REQUIRE(zeroed.cols() == 3);
+ for (std::size_t i = 0; i < zeroed.rows(); ++i) {
+ for (std::size_t j = 0; j < zeroed.cols(); ++j) {
+ CHECK(zeroed(i, j) == 0.0);
+ }
+ }
+
+ const Matrix filled(2, 2, 1.5);
+ CHECK(filled(0, 0) == 1.5);
+ CHECK(filled(0, 1) == 1.5);
+ CHECK(filled(1, 0) == 1.5);
+ CHECK(filled(1, 1) == 1.5);
+}
+
+TEST_CASE("Matrix supports checked element access and fill", "[matrix]") {
+ Matrix a(2, 3);
+ a(0, 0) = 1.0;
+ a(0, 1) = 2.0;
+ a(1, 2) = 5.0;
+
+ CHECK(a(0, 0) == 1.0);
+ CHECK(a(0, 1) == 2.0);
+ CHECK(a(1, 2) == 5.0);
+
+ a.fill(-3.0);
+ for (std::size_t i = 0; i < a.rows(); ++i) {
+ for (std::size_t j = 0; j < a.cols(); ++j) {
+ CHECK(a(i, j) == -3.0);
+ }
+ }
+
+ CHECK_THROWS_AS(a(2, 0), std::out_of_range);
+ CHECK_THROWS_AS(a(0, 3), std::out_of_range);
+}
+
+TEST_CASE("Matrix identity and zero factories work", "[matrix]") {
+ const Matrix identity = Matrix::identity(3);
+ REQUIRE(identity.rows() == 3);
+ REQUIRE(identity.cols() == 3);
+
+ for (std::size_t i = 0; i < identity.rows(); ++i) {
+ for (std::size_t j = 0; j < identity.cols(); ++j) {
+ const double expected = (i == j) ? 1.0 : 0.0;
+ CHECK(identity(i, j) == expected);
+ }
+ }
+
+ const Matrix zeros = Matrix::zeros(2, 4);
+ CHECK(zeros.rows() == 2);
+ CHECK(zeros.cols() == 4);
+ for (std::size_t i = 0; i < zeros.rows(); ++i) {
+ for (std::size_t j = 0; j < zeros.cols(); ++j) {
+ CHECK(zeros(i, j) == 0.0);
+ }
+ }
+}
+
+TEST_CASE("Matrix copy, move, and initializer list keeps row-major layout", "[matrix]") {
+ static_assert(std::is_nothrow_move_constructible_v<Matrix>);
+
+ Matrix original{{1.0, 2.0}, {3.0, 4.0}};
+ Matrix copied = original;
+ copied(0, 0) = 10.0;
+
+ CHECK(original(0, 0) == 1.0);
+ CHECK(copied(0, 0) == 10.0);
+
+ Matrix moved = std::move(original);
+ REQUIRE(moved.rows() == 2);
+ REQUIRE(moved.cols() == 2);
+ CHECK(moved(0, 0) == 1.0);
+ CHECK(moved(0, 1) == 2.0);
+ CHECK(moved(1, 0) == 3.0);
+ CHECK(moved(1, 1) == 4.0);
+}
+
+TEST_CASE("Matrix initializer list rejects unequal row lengths", "[matrix]") {
+ CHECK_THROWS_AS((Matrix{{1.0, 2.0}, {3.0}}), std::invalid_argument);
+}
diff --git a/tests/test_vector.cpp b/tests/test_vector.cpp
new file mode 100644
index 0000000..1d7b91a
--- /dev/null
+++ b/tests/test_vector.cpp
@@ -0,0 +1,88 @@
+#include "linalg_error.hpp"
+#include "vector.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <type_traits>
+#include <utility>
+
+using linalg::DimensionMismatchError;
+using linalg::Vector;
+
+TEST_CASE("Vector constructors initialize size and values", "[vector]") {
+ const Vector empty;
+ REQUIRE(empty.size() == 0);
+ REQUIRE(empty.empty());
+
+ const Vector zeroed(4);
+ REQUIRE(zeroed.size() == 4);
+ for (std::size_t i = 0; i < zeroed.size(); ++i) {
+ CHECK(zeroed[i] == 0.0);
+ }
+
+ const Vector filled(3, 2.5);
+ REQUIRE(filled.size() == 3);
+ CHECK(filled[0] == 2.5);
+ CHECK(filled[1] == 2.5);
+ CHECK(filled[2] == 2.5);
+}
+
+TEST_CASE("Vector supports checked element access and fill", "[vector]") {
+ Vector v{1.0, 2.0, 3.0};
+ REQUIRE(v.size() == 3);
+
+ v[1] = 7.0;
+ CHECK(v[0] == 1.0);
+ CHECK(v[1] == 7.0);
+ CHECK(v[2] == 3.0);
+
+ v.fill(-2.0);
+ CHECK(v[0] == -2.0);
+ CHECK(v[1] == -2.0);
+ CHECK(v[2] == -2.0);
+
+ CHECK_THROWS_AS(v[3], std::out_of_range);
+}
+
+TEST_CASE("Vector copy and move preserve values", "[vector]") {
+ static_assert(std::is_nothrow_move_constructible_v<Vector>);
+
+ Vector original{4.0, -1.0, 8.0};
+ Vector copied = original;
+ copied[0] = 10.0;
+
+ CHECK(original[0] == 4.0);
+ CHECK(copied[0] == 10.0);
+
+ Vector moved = std::move(original);
+ REQUIRE(moved.size() == 3);
+ CHECK(moved[0] == 4.0);
+ CHECK(moved[1] == -1.0);
+ CHECK(moved[2] == 8.0);
+}
+
+TEST_CASE("Vector arithmetic enforces shape compatibility", "[vector]") {
+ const Vector a{1.0, 2.0, 3.0};
+ const Vector b{4.0, 5.0, 6.0};
+
+ const Vector sum = a + b;
+ CHECK(sum[0] == 5.0);
+ CHECK(sum[1] == 7.0);
+ CHECK(sum[2] == 9.0);
+
+ const Vector diff = b - a;
+ CHECK(diff[0] == 3.0);
+ CHECK(diff[1] == 3.0);
+ CHECK(diff[2] == 3.0);
+
+ const Vector scaled = 0.5 * a;
+ CHECK(scaled[0] == 0.5);
+ CHECK(scaled[1] == 1.0);
+ CHECK(scaled[2] == 1.5);
+
+ CHECK(linalg::dot(a, b) == 32.0);
+
+ const Vector short_vec{1.0, 2.0};
+ CHECK_THROWS_AS(a + short_vec, DimensionMismatchError);
+ CHECK_THROWS_AS(linalg::dot(a, short_vec), DimensionMismatchError);
+}