From 7d763b5778cfbf81ee420558d13e8acbb66d860a Mon Sep 17 00:00:00 2001 From: Jasper Date: Sun, 31 Aug 2025 15:47:31 +0200 Subject: Initial commit --- .gitignore | 3 + build.sh | 5 + debug.c | 12 +++ matrix.c | 339 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ matrix.h | 49 +++++++++ tests.c | 24 +++++ tests.h | 7 ++ utils.c | 39 +++++++ utils.h | 9 ++ 9 files changed, 487 insertions(+) create mode 100644 .gitignore create mode 100755 build.sh create mode 100644 debug.c create mode 100644 matrix.c create mode 100644 matrix.h create mode 100644 tests.c create mode 100644 tests.h create mode 100644 utils.c create mode 100644 utils.h diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2eb88f3 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +misc/ +bin/* +tags diff --git a/build.sh b/build.sh new file mode 100755 index 0000000..6ec6b63 --- /dev/null +++ b/build.sh @@ -0,0 +1,5 @@ +if [ ! -d "bin" ]; then + mkdir bin +fi + +gcc -Wall -Wextra -pedantic -lm -o bin/linalg_debug debug.c tests.c utils.c matrix.c diff --git a/debug.c b/debug.c new file mode 100644 index 0000000..df2071b --- /dev/null +++ b/debug.c @@ -0,0 +1,12 @@ +#include +#include +#include + +#include "tests.h" +#include "matrix.h" + +int main() +{ + srand(time(NULL)); + run_tests(); +} diff --git a/matrix.c b/matrix.c new file mode 100644 index 0000000..46302f7 --- /dev/null +++ b/matrix.c @@ -0,0 +1,339 @@ +#include +#include +#include +#include +#include +#include + +#include "matrix.h" +#include "utils.h" + +Matrix *matrix_alloc(const size_t m, const size_t n) +{ + double *xs = malloc(m * n * sizeof(double)); + assert(xs != NULL && "Allocation of matrix elements failed"); + Matrix *mat = malloc(sizeof(Matrix)); + assert(mat != NULL && "Allocation of matrix failed"); + mat->xs = xs; + mat->m = m; + mat->n = n; + return mat; +} + +void matrix_free(Matrix *mat) +{ + free(mat->xs); + free(mat); +} + +void matrix_freen(Matrix **mats) +{ + Matrix *mat, **p; + p = mats; + while ((mat = *(p++))) + matrix_free(mat); + free(mats); +} + +void matrix_print(const Matrix *mat) +{ + for (size_t i = 0; i < mat->m; ++i) + { + for (size_t j = 0; j < mat->n; ++j) + printf("%lf ", matrix_at(mat, i, j)); + printf("\n"); + } +} + +/* Expected format of a 3x2 matrix [3 2;1 2;0 0] + * NOTE: No spaces around ';' + * NOTE: No proper error handling is done */ +Matrix *matrix_from_str(char *str) +{ + assert(str != NULL && "Passed null string"); + char *row, *tmp_row, *entry; + + str = str_delete_at(str, 0); + str = str_delete_at(str, strlen(str) - 1); + + size_t rows = str_count_occ(str, ';') + 1; + size_t cols = (str_count_occ(str, ' ') + rows) / rows; + + Matrix *mat_tmp = matrix_const(rows, cols, 0); + + for (size_t i = 0; (row = strsep(&str, ";")); ++i) + { + tmp_row = strdup(row); + for (size_t j = 0; (entry = strsep(&tmp_row, " ")); ++j) + matrix_at(mat_tmp, i, j)= atof(entry); + } + + free(tmp_row); + return mat_tmp; +} + +Matrix *matrix_from_arr(double arr[], size_t m, size_t n) +{ + Matrix *mat = matrix_alloc(m, n); + + for (size_t i = 0; i < m * n; ++i) + mat->xs[i] = arr[i]; + + return mat; +} + +char *matrix_to_str(const Matrix *mat) +{ + assert(mat != NULL); + return "TODO"; +} + +Matrix *matrix_id(const size_t n) +{ + Matrix *tmp = matrix_const(n, n, 0); + + for (size_t i = 0; i < n; ++i) + matrix_at(tmp, i, i) = 1; + + return tmp; +} + +Matrix *matrix_const(const size_t m, const size_t n, const double x) +{ + Matrix *mat = matrix_alloc(m, n); + for (size_t i = 0; i < mat->m; ++i) + { + for (size_t j = 0; j < mat->n; ++j) + matrix_at(mat, i, j) = x; + } + return mat; +} + +Matrix *matrix_copy (const Matrix *mat) +{ + Matrix *tmp = matrix_alloc(mat->m, mat->n); + + for (size_t i = 0; i < mat->m; ++i) + { + for (size_t j = 0; j < mat->n; ++j) + matrix_at(tmp, i, j) = matrix_at(mat, i, j); + } + return tmp; +} + +double matrix_trace(const Matrix *mat) +{ + assert(matrix_is_square(mat) && "Trace is not defined for non-square matrices"); + double sum = 0; + + for (size_t i = 0; i < mat->n; ++i) + sum += matrix_at(mat, i, i); + + return sum; +} + +Matrix *matrix_transpose(const Matrix *mat) +{ + assert(mat != NULL && "Matrix points to null"); + assert(mat->xs != NULL && "Matrix elements point to null"); + + Matrix *t = matrix_alloc(mat->n, mat->m); + for (size_t i = 0; i < mat->n; ++i) + { + for (size_t j = 0; j < mat->m; ++j) + matrix_at(t, i, j) = matrix_at(mat, j, i); + } + return t; +} + +Matrix *matrix_add(const Matrix *A, const Matrix *B) +{ + assert(A->m == B->m && A->n == B->n && "Dimension mismatch"); + Matrix *C = matrix_alloc(A->m, A->n); + + for (size_t i = 0; i < A->m; ++i) + { + for (size_t j = 0; j < A->n; ++j) + matrix_at(C, i, j) = matrix_at(A, i, j) + matrix_at(B, i, j); + } + + return C; +} + +Matrix *matrix_scale(const double x, const Matrix *A) +{ + Matrix *tmp = matrix_alloc(A->m, A->n); + + for (size_t i = 0; i < A->m; ++i) + { + for (size_t j = 0; j < A->n; ++j) + matrix_at(tmp, i, j) = x * matrix_at(A, i, j); + } + + return tmp; +} + +Matrix *matrix_sub(const Matrix *A, const Matrix *B) +{ + return matrix_add(A, matrix_scale(-1, B)); +} + +Matrix *matrix_mult(const Matrix *A, const Matrix *B) +{ + assert(A->n == B->m && "Dimension mismatch"); + Matrix *C = matrix_const(A->m, B->n, 0); + + for (size_t i = 0; i < C->m; ++i) + { + for (size_t j = 0; j < C->n; ++j) + { + for (size_t k = 0; k < A->n; ++k) + matrix_at(C, i, j) += matrix_at(A, i, k) * matrix_at(B, k, j); + } + } + + return C; +} + +Matrix *matrix_rand(const size_t m, const size_t n, const int bound_l, const int bound_u, MatrixType type) +{ + assert(bound_l <= bound_u && "Lower bound is expected to be less or equal than upper bound"); + + if (bound_l == bound_u) + { + return matrix_const(m, n, bound_l); + } + + Matrix *mat = matrix_const(m, n, 0); + + for (size_t i = 0; i < mat->m; ++i) + { + for (size_t j = 0; j < mat->n; ++j) + { + switch (type) + { + case MATRIX_DIAG: + if (i != j) continue; + break; + case MATRIX_TRIAG_UPPER: + if (i > j) continue; + break; + case MATRIX_TRIAG_LOWER: + if (i < j) continue; + break; + case MATRIX_TRIAG_SUPPER: + if (i >= j) continue; + break; + case MATRIX_TRIAG_SLOWER: + if (i <= j) continue; + break; + case MATRIX_NONE: break; + } + matrix_at(mat, i, j) = int_rand(bound_l, bound_u); + } + } + + return mat; +} + +int matrix_is_square(const Matrix *mat) +{ + return mat->m == mat->n; +} + +int matrix_eq(const Matrix *A, const Matrix *B, const double tol) +{ + assert(A->m == B->m && A->n == B->n && "Dimension mismatch"); + for (size_t i = 0; i < A->m; ++i) + { + for (size_t j = 0; j < A->n; ++j) + if (matrix_at(A, i, j) - matrix_at(B, i, j) > tol) return FALSE; + } + return TRUE; +} + +double matrix_norm_frob(const Matrix *mat) +{ + return sqrt(matrix_trace(matrix_mult(matrix_transpose(mat), mat))); +} + +Matrix **matrix_LR(const Matrix *A, const Matrix *b) +{ + assert(matrix_is_square(A) && "LR decomposition only works for square matrices"); + assert(A->n == b->m && "Dimension mismatch"); + assert(b->n == 1 && "Vector b expected"); + + Matrix *L = matrix_id(A->n); + Matrix *A_prev = matrix_copy(A); + Matrix *A_curr = matrix_const(A->n, A->n, 0); + + for (size_t k = 1; k < A->n; ++k) + { + for (size_t i = 0; i < k; ++i) + { + for (size_t j = 0; j < A->n; ++j) + matrix_at(A_curr, i, j) = matrix_at(A_prev, i, j); + + } + + assert(matrix_at(A_prev, k - 1, k - 1) != 0 && "Pivoting strategy needed"); + + for (size_t i = k; i < A->n; ++i) + { + matrix_at(L, i, k - 1) = matrix_at(A_prev, i, k - 1) / matrix_at(A_prev, k - 1, k - 1); + for (size_t j = k; j < A->n; ++j) + { + matrix_at(A_curr, i, j) = matrix_at(A_prev, i, j) - matrix_at(L, i, k - 1) * matrix_at(A_prev, k - 1, j); + } + } + + A_prev = matrix_copy(A_curr); + if (k != A->n - 1) + A_curr = matrix_const(A->n, A->n, 0); + } + + matrix_free(A_prev); + Matrix **result = malloc(sizeof(Matrix) * 2); + result[0] = L; + result[1] = A_curr; + result[2] = NULL; + return result; +} + +Matrix *matrix_forwardel(const Matrix *L, const Matrix *b) +{ + assert(L->n == b->m && "Dimension mismatch"); + assert(b->n == 1 && "Vector b expected"); + + double sum; + Matrix *y = matrix_const(L->n, 1, 0); + + for (size_t i = 0; i < y->m; ++i) + { + sum = matrix_at(b, i, 0); + for (size_t j = 0; j < i; ++j) + sum -= matrix_at(L, i, j) * matrix_at(y, j, 0); + matrix_at(y, i, 0) = sum / matrix_at(L, i, i); + } + + return y; +} + +Matrix *matrix_backsubst(const Matrix *R, const Matrix *y) +{ + assert(R->n == y->m && "Dimension mismatch"); + assert(y->n == 1 && "Vector y expected"); + Matrix *x = matrix_const(R->n, 1, 0); + + double sum; + + for (int i = x->m; i >= 0; --i) + { + sum = matrix_at(y, i, 0); + for (size_t j = i + 1; j < x->m; ++j) + sum -= matrix_at(R, i, j) * matrix_at(x, j, 0); + matrix_at(x, i, 0) = sum / matrix_at(R, i, i); + } + + return x; +} diff --git a/matrix.h b/matrix.h new file mode 100644 index 0000000..e066777 --- /dev/null +++ b/matrix.h @@ -0,0 +1,49 @@ +#ifndef MATRIX_H +#define MATRIX_H + +#include + +#define matrix_at(mat, i, j) ((mat)->xs[(i) * (mat)->n + (j)]) +#define TRUE 1 +#define FALSE 0 + +typedef enum _MatrixType { + MATRIX_DIAG, + MATRIX_TRIAG_UPPER, + MATRIX_TRIAG_LOWER, + MATRIX_TRIAG_SUPPER, /* S for strict */ + MATRIX_TRIAG_SLOWER, + MATRIX_NONE +} MatrixType; + +typedef struct _Matrix { + double *xs; + size_t m; + size_t n; +} Matrix; + +Matrix *matrix_alloc(const size_t m, const size_t n); +void matrix_free(Matrix *mat); +void matrix_freen(Matrix **mats); +void matrix_print(const Matrix *mat); +Matrix *matrix_from_str(char *str); +Matrix *matrix_from_arr(double arr[], size_t m, size_t n); +char *matrix_to_str(const Matrix *mat); +Matrix *matrix_id(const size_t n); +Matrix *matrix_const(const size_t m, const size_t n, const double x); +Matrix *matrix_copy (const Matrix *mat); +double matrix_trace(const Matrix *mat); +Matrix *matrix_transpose(const Matrix *mat); +Matrix *matrix_add(const Matrix *A, const Matrix *B); +Matrix *matrix_scale(const double x, const Matrix *A); +Matrix *matrix_sub(const Matrix *A, const Matrix *B); +Matrix *matrix_mult(const Matrix *A, const Matrix *B); +Matrix *matrix_rand(const size_t m, const size_t n, const int bound_l, const int bound_u, MatrixType type); +int matrix_is_square(const Matrix *mat); +int matrix_eq(const Matrix *A, const Matrix *B, const double tol); +double matrix_norm_frob(const Matrix *mat); +Matrix **matrix_LR(const Matrix *A, const Matrix *b); +Matrix *matrix_forwardel(const Matrix *L, const Matrix *b); +Matrix *matrix_backsubst(const Matrix *R, const Matrix *y); + +#endif diff --git a/tests.c b/tests.c new file mode 100644 index 0000000..9b2a4dc --- /dev/null +++ b/tests.c @@ -0,0 +1,24 @@ +#include "tests.h" + +#include + +#include "matrix.h" + +void test_LR() +{ + Matrix *A = matrix_from_str("[2 3 1;4 1 -5;-1 2 3]"); + Matrix *b = matrix_from_str("[-1;2;3]"); + + Matrix **LR = matrix_LR(A, b); + Matrix *L = LR[0]; + Matrix *R = LR[1]; + Matrix *y = matrix_forwardel(L, b); + Matrix *x = matrix_backsubst(R, y); + + assert(matrix_eq(matrix_mult(A, x), b, 0.01)); +} + +void run_tests() +{ + test_LR(); +} diff --git a/tests.h b/tests.h new file mode 100644 index 0000000..65e82d5 --- /dev/null +++ b/tests.h @@ -0,0 +1,7 @@ +#ifndef TEST_H +#define TEST_H + +void test_LR(); +void run_tests(); + +#endif diff --git a/utils.c b/utils.c new file mode 100644 index 0000000..dff3b3f --- /dev/null +++ b/utils.c @@ -0,0 +1,39 @@ +#include +#include +#include + +char *str_delete_at(const char *str, const int pos) +{ + char *tmp = malloc(sizeof(char) * (strlen(str) + 1)); + assert(tmp != NULL); + strcpy(tmp, str); + + int len = strlen(str); + if (pos >= len) + return tmp; + + for (int i = 0; i < len; ++i) + { + if (i >= pos) + tmp[i] = tmp[i + 1]; + } + + return tmp; +} + +size_t str_count_occ(const char *str, const char ch) +{ + size_t count = 0; + for (int i = 0; str[i] != '\0'; ++i) + { + if (str[i] == ch) + count += 1; + } + return count; +} + +/* Return random integer between lower bound bound_l (inclusive) and upper bound bound_u (exclusive) */ +int int_rand(const int bound_l, const int bound_u) +{ + return bound_l + (rand() % bound_u); +} diff --git a/utils.h b/utils.h new file mode 100644 index 0000000..8503c1b --- /dev/null +++ b/utils.h @@ -0,0 +1,9 @@ +#ifndef UTIL_H +#define UTIL_H + +char *str_delete_at(const char *str, const int pos); +size_t str_count_occ(const char *str, const char ch); + +int int_rand(const int bound_l, const int bound_u); + +#endif -- cgit v1.2.3