Initial commit

1 files changed, 339 insertions, 0 deletions

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;
+}