From 7d763b5778cfbf81ee420558d13e8acbb66d860a Mon Sep 17 00:00:00 2001 From: Jasper Date: Sun, 31 Aug 2025 15:47:31 +0200 Subject: Initial commit --- matrix.c | 339 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 339 insertions(+) create mode 100644 matrix.c (limited to 'matrix.c') 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; +} -- cgit v1.2.3