Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 0 additions & 15 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,12 @@ CXXFLAGS = -std=c++2b -Wall -Wextra -pedantic
SRC_MAIN = tests/main.cpp
SRC_FMAT = tests/fmat.cpp
BIN_MAIN = tests/main.o
BIN_FMAT = tests/fmat.o

build:
$(CXX) $(CXXFLAGS) $(SRC_MAIN) -o $(BIN_MAIN)

build-fmat:
$(CXX) $(CXXFLAGS) $(SRC_FMAT) -o $(BIN_FMAT)

test: build
./$(BIN_MAIN)

test-fmat: build-fmat
./$(BIN_FMAT)

test-all: test test-fmat

test-asan:
$(CXX) $(CXXFLAGS) -fsanitize=address,undefined -fno-omit-frame-pointer -g $(SRC_MAIN) -o $(BIN_MAIN)
./$(BIN_MAIN)
$(CXX) $(CXXFLAGS) -fsanitize=address,undefined -fno-omit-frame-pointer -g $(SRC_FMAT) -o $(BIN_FMAT)
./$(BIN_FMAT)

clean:
rm -f $(BIN_MAIN) $(BIN_FMAT)
60 changes: 9 additions & 51 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,58 +1,16 @@
## JSF Compliant Linear Algebra Engine
## Features
### Core Matrix & Vector Mathematics
* **N-Dimensional Vector Spaces:** Full implementation of multi-dimensional vector operations, including Euclidean norms, dot products, cross products, and vector projections.

### Matrix & Vector Algebra
* **Comprehensive Matrix Algebra:** Complete support for basic arithmetic (addition, subtraction, scalar multiplication) and advanced matrix operations (matrix-matrix multiplication, transposition, determinants, matrix inversion, and trace evaluation).
* **Structured Matrix Types:** Optimized storage and algorithms for specialized matrix representations, including dense, sparse, symmetric, band, tridiagonal, and identity structures.

### Numerical Methods & Factorizations
* **Systems of Linear Equations:** Direct and iterative solvers for linear systems $Ax = b$, optimized for computational efficiency and error propagation mitigation.
* **Matrix Reduction:** Stable Gaussian Elimination with partial and complete pivoting mechanisms to prevent catastrophic cancellation and numerical instability.
* **Matrix Factorizations:** Advanced matrix decompositions for structural analysis and efficient backend calculation:
* **LU Decomposition:** (with pivoting, $PA = LU$) for fast subsequent solves.
* **QR Factorization:** via Householder reflections or Gram-Schmidt orthonormalization for least-squares approximations.
* **Cholesky Decomposition:** for symmetric positive-definite systems.
* **Structured Matrix Types:** Optimized storage and algorithms for specialized matrix representations.

## Project Structure
```text
linear-algebra-engine/
├── .github/ # CI/CD pipelines and compilation workflows
├── build/ # Compiled binaries, static archives, and object targets
├── docs/ # API specifications, mathematical proofs, and JSF matrices
├── include/ # Public API headers (.hpp / .h)
│ ├── core/ # Fundamental vector and matrix definitions
│ │ ├── matrix.hpp # Matrix class template and structural definition
│ │ └── vector.hpp # N-dimensional vector template and inline arithmetic
│ ├── memory/ # Safety-critical allocation layer
│ │ ├── arena_allocator.hpp # Monolithic bump allocator implementation
│ │ └── memory_block.hpp # Statically sized arena block templates
│ ├── methods/ # Numerical solver implementations
│ │ ├── decomposition.hpp # LU, QR, and Cholesky factorization interfaces
│ │ ├── elimination.hpp # Gaussian elimination and pivoting logic
│ │ └── solvers.hpp # Direct and iterative system solvers (Ax = b)
│ └── structures/ # Specialized matrix structural layouts
│ ├── sparse_matrix.hpp # Compressed Sparse Row (CSR) implementation
│ └── tridiagonal.hpp # Highly optimized tridiagonal matrix storage
├── src/ # Source implementation files (.cpp / .c)
│ ├── core/
│ │ ├── matrix.cpp
│ │ └── vector.cpp
│ ├── memory/
│ │ └── arena_allocator.cpp
│ └── methods/
│ ├── decomposition.cpp
│ ├── elimination.cpp
│ └── solvers.cpp
├── demos/ # Reference implementations and modeling sandboxes
│ ├── ode_solver/
│ │ ├── main.cpp # Main executable for RK4 mass-spring system
│ │ └── state_space.cpp # Linear state-space vector differential maps
│ └── pde_heat_diffusion/
│ └── heat_solver.cpp # Finite difference matrix mesh simulation
├── tests/ # Unit testing suite and validation frameworks
│ ├── test_matrix_ops.cpp # Validation profiles for algebraic compliance
│ ├── test_solvers.cpp # Accuracy verification against ill-conditioned systems
│ └── test_memory.cpp # Boundary testing for bump allocator allocations
├── CMakeLists.txt # Main project build orchestration script
├── LICENSE # Licensing documentation
└── README.md # Project documentation overview
linear-algebra/
├── lib/
├── src/
├── tests/
└── docs/
```
41 changes: 0 additions & 41 deletions lib/arena.h

This file was deleted.

128 changes: 88 additions & 40 deletions lib/matrix.h
Original file line number Diff line number Diff line change
@@ -1,99 +1,148 @@
#include <iostream>
#include <vector>

using Dimension = std::pair<size_t, size_t>;

template<typename T>
using Array2D = std::vector<std::vector<T>>;

template<typename T>
using Array = std::vector<T>;

template<typename T>
class Matrix;

template<typename T>
Dimension matrix_size(const Matrix<T>& mat);
template<typename T>
void matrix_fill(Matrix<T>& mat, const T& data);
template<typename T>
void matrix_eq(Matrix<T>& mat, const Matrix<T>& eq_mat);
template<typename T>
void matrix_add(Matrix<T>& mat, const Matrix<T>& add_mat);
template<typename T>
void matrix_subtr(Matrix<T>& mat, const Matrix<T>& subtr_mat);
template<typename T>
void matrix_mul(Matrix<T>& mat, const Matrix<T>& mul_mat);
template<typename T>
void matrix_set(Matrix<T>& mat, const size_t& row, const size_t& col, const T& data);


template<typename T>
class Matrix {
public:
Matrix(size_t numRows, size_t numCols, const T& data);
public:

Matrix(size_t numRows, size_t numCols, const T& data);
Matrix(size_t numRows, size_t numCols, const std::vector<std::vector<T>>& array_2d);
Matrix(const Matrix<T>& mat);
~Matrix();

friend std::pair<size_t, size_t> & matrix_size(const Matrix<T>& mat);

friend void matrix_fill(Matrix<T>& mat, const T& data);
friend void matrix_eq(Matrix<T>& mat, const Matrix<T>& eq_mat);
friend void matrix_add(Matrix<T>& mat, const Matrix<T>& add_mat);
friend void matrix_subtr(Matrix<T> &mat, const Matrix<T>& subtr_mat);
friend void matrix_mul(Matrix<T>& mat, const Matrix<T>& mul_mat);
std::vector<T>& operator[](size_t row) { return data[row]; }
const std::vector<T>& operator[](size_t row) const { return data[row]; }

friend T& matrix_get(const Matrix<T>& mat, const size_t& row, const size_t& col);
friend void matrix_set(Matrix<T>& mat, const size_t& row, const size_t& col, const T& data);
friend Dimension matrix_size<>(const Matrix<T>& mat);

friend void matrix_fill<>(Matrix<T>& mat, const T& data);
friend void matrix_eq<>(Matrix<T>& mat, const Matrix<T>& eq_mat);
friend void matrix_add<>(Matrix<T>& mat, const Matrix<T>& add_mat);
friend void matrix_subtr<>(Matrix<T>& mat, const Matrix<T>& subtr_mat);
friend void matrix_mul<>(Matrix<T>& mat, const Matrix<T>& mul_mat);
friend void matrix_set<>(Matrix<T>& mat, const size_t& row, const size_t& col, const T& data);

private:
size_t numRows;
size_t numCols;
std::vector<T> data;
std::vector<std::vector<T>> data;
};


template<typename T>
Matrix<T>::Matrix(size_t numRows, size_t numCols, const T& data) :
numRows(numRows), numCols(numCols), data(std::vector<std::vector<T>>(numRows, std::vector<T>(numCols, data))) {}

template<typename T>
Matrix<T>::Matrix(const Matrix<T>& mat) {
std::pair<size_t, size_t> sz_mat = matrix_size(mat);
Dimension sz_mat = matrix_size(mat);
numRows = sz_mat.first;
numCols = sz_mat.second;

data = std::vector<std::vector<T>>(numRows, std::vector<T>(numCols));

for(size_t row = 0; row <= numRows-1; row++) {
for(size_t col = 0; col <= numCols-1; col++) {
data[row][col] = mat[row][col];
for(size_t row = 0; row < numRows; row++) {
for(size_t col = 0; col < numCols; col++) {
data[row][col] = mat[row][col];
}
}
}

template<typename T>
Matrix<T>::~Matrix() = default;

template<typename T>
Dimension matrix_size(const Matrix<T>& mat) {
return { mat.numRows, mat.numCols };
}

template<typename T>
void matrix_set(Matrix<T>& mat, const size_t& row, const size_t& col, const T& value) {
mat.data[row][col] = value;
}

template<typename T>
Matrix<T>::Matrix(size_t numRows, size_t numCols, const std::vector<std::vector<T>>& vec_data) :
numRows(numRows), numCols(numCols), data(vec_data) {}

template<typename T>
void matrix_fill(Matrix<T> &mat, const T& val) {
std::pair<size_t, size_t> sz_mat = matrix_size(mat);
Dimension sz_mat = matrix_size(mat);

const size_t& numRows = sz_mat.first;
const size_t& numCols = sz_mat.second;

for(size_t row = 0; row <= numRows-1; row++) {
for(size_t col = 0; col <= numCols-1; col++) {
mat.data[row][col] = val;
for(size_t row = 0; row < numRows; row++) {
for(size_t col = 0; col < numCols; col++) {
matrix_set(mat, row, col, val);
}
}
}

template<typename T>
void matrix_add(Matrix<T> &mat, const Matrix<T>& add_mat) {
const std::pair<size_t, size_t> & sz_mat = matrix_size(mat);
const std::pair<size_t, size_t> & sz_add_mat = matrix_size(add_mat);
const Dimension & sz_mat = matrix_size(mat);
const Dimension & sz_add_mat = matrix_size(add_mat);

if (sz_mat.first != sz_add_mat.first) return;
if (sz_mat.second != sz_add_mat.second) return;

for(size_t row = 0; row <= sz_mat.first-1; row++) {
for(size_t col = 0; col <= sz_mat.second-1; col++) {
mat.data[row][col] += (add_mat[row][col]);
for(size_t row = 0; row < sz_mat.first; row++) {
for(size_t col = 0; col < sz_mat.second; col++) {
mat[row][col] += add_mat[row][col];
}
}
}

template<typename T>
void matrix_subtr(Matrix<T> &mat, const Matrix<T>& subtr_mat) {
const std::pair<size_t, size_t> & sz_mat = matrix_size(mat);
const std::pair<size_t, size_t> & sz_subtr_mat = matrix_size(subtr_mat);
const Dimension & sz_mat = matrix_size(mat);
const Dimension & sz_subtr_mat = matrix_size(subtr_mat);

if (sz_mat.first != sz_subtr_mat.first) return;
if (sz_mat.second != sz_subtr_mat.second) return;

for(size_t row = 0; row <= sz_mat.first-1; row++) {
for(size_t col = 0; col <= sz_mat.second-1; col++) {
mat.data[row][col] -= (subtr_mat[row][col]);
for(size_t row = 0; row < sz_mat.first; row++) {
for(size_t col = 0; col < sz_mat.second; col++) {
mat[row][col] -= subtr_mat[row][col];
}
}
}


template<typename T>
void matrix_mul(Matrix<T> &mat, const Matrix<T>& mul_mat) {
const Matrix<T>& ref(mat);
const Matrix<T> lhs(mat);

const std::pair<size_t, size_t> & sz_mat = matrix_size(mat);
const std::pair<size_t, size_t> & sz_mul_mat = matrix_size(mul_mat);
const Dimension & sz_mat = matrix_size(mat);
const Dimension & sz_mul_mat = matrix_size(mul_mat);

if(sz_mat.second != sz_mul_mat.first) return;

Expand All @@ -102,14 +151,13 @@ void matrix_mul(Matrix<T> &mat, const Matrix<T>& mul_mat) {

mat.data = std::vector<std::vector<T>>(mat.numRows, std::vector<T>(mat.numCols));

for(size_t row = 0; row <= mat.numRows-1; row++) {
for(size_t col = 0; col <= mat.numCols-1; col++) {
mat.data[row][col] = 0;
for(size_t i = 0; i <= sz_mat.first; i++) {
for(size_t j = 0; j <= sz_mat.second; j++) {
mat.data[row][col] += (matrix_get(ref, i, j) * matrix_get(mul_mat, j, i));
}
for(size_t row = 0; row < mat.numRows; row++) {
for(size_t col = 0; col < mat.numCols; col++) {
T sum = T{};
for(size_t k = 0; k < sz_mat.second; k++) {
sum += (lhs[row][k] * mul_mat[k][col]);
}
mat[row][col] = sum;
}
}
}
Loading