Initial commit

9 files changed, 487 insertions, 0 deletions

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 <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#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 <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+
+#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 <stddef.h>
+
+#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 <assert.h>
+
+#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 <string.h>
+#include <assert.h>
+#include <stdlib.h>
+
+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