From 0f617454d39d4bf1d91ed427a5bae09325140f5e Mon Sep 17 00:00:00 2001
From: Al <albarrentine@gmail.com>
Date: Mon, 31 Aug 2015 14:57:11 -0400
Subject: [PATCH] [math] Dense matrices

---
 src/matrix.c | 122 +++++++++++++++++++++++++++++++++++++++++++++++++++
 src/matrix.h |  34 ++++++++++++++
 2 files changed, 156 insertions(+)
 create mode 100644 src/matrix.c
 create mode 100644 src/matrix.h

diff --git a/src/matrix.c b/src/matrix.c
new file mode 100644
index 00000000..d7a4de11
--- /dev/null
+++ b/src/matrix.c
@@ -0,0 +1,122 @@
+#include "matrix.h"
+
+
+matrix_t *matrix_new(size_t m, size_t n) {
+    matrix_t *matrix = malloc(sizeof(matrix_t));
+
+    if (matrix == NULL) {
+        return NULL;
+    }
+
+    matrix->m = m;
+    matrix->n = n;
+
+    matrix->values = malloc(sizeof(double) * m * n);
+    if (matrix->values == NULL) {
+        free(matrix);
+        return NULL;
+    }
+
+    return matrix;
+
+}
+
+matrix_t *matrix_copy(matrix_t *self) {
+    matrix_t *cpy = matrix_new(self->m, self->n);
+    size_t num_values = self->m * self->n;
+    memcpy(cpy->values, self->values, num_values * sizeof(double));
+    return cpy;
+}
+
+inline void matrix_init_values(matrix_t *self, double *values) {
+    size_t num_values = self->m * self->n;
+    memcpy(self->values, values, num_values * sizeof(double));
+}
+
+inline void matrix_set(matrix_t *self, double value) {
+    double_array_set(self->values, self->m * self->n, value);
+}
+
+inline void matrix_set_row(matrix_t *self, size_t index, double *row) {
+    size_t offset = index * self->m;
+    double *values = self->values;
+    size_t n = self->n;
+    memcpy(values + offset, row, n * sizeof(double));
+}
+
+inline void matrix_set_scalar(matrix_t *self, size_t row_index, size_t col_index, double value) {
+    size_t offset = row_index * self->m + col_index;
+    double *values = self->values;
+    values[offset] = value;
+}
+
+inline matrix_t *matrix_new_value(size_t m, size_t n, double value) {
+    matrix_t *matrix = matrix_new(m, n);
+    matrix_set(matrix, value);
+    return matrix;
+}
+
+inline matrix_t *matrix_new_zeros(size_t m, size_t n) {
+    return matrix_new_value(m, n, 0.0);
+}
+
+inline matrix_t *matrix_new_ones(size_t m, size_t n) {
+    return matrix_new_value(m, n, 1.0);
+}
+
+matrix_t *matrix_new_values(size_t m, size_t n, double *values) {
+    matrix_t *matrix = matrix_new(m, n);
+    memcpy(matrix->values, values, m * n * sizeof(double));
+    return matrix;
+}
+
+inline void matrix_scale(matrix_t *self, double value) {
+    double_array_mul(self->values, self->m * self->n, value);
+}
+
+void matrix_dot_vector(matrix_t *self, double *vec, double *result) {
+    double *values = self->values;
+    size_t n = self->n;
+    for (int i = 0; i < self->m; i++) {
+        for (int j = 0; j < n; j++) {
+            result[i] += values[n * i + j] * vec[j];
+        }
+    }
+}
+
+int matrix_dot_matrix(matrix_t *m1, matrix_t *m2, matrix_t *result) {
+    if (m1->n != m2->m || m1->m != result->m || m2->n != result->n) {
+        return -1;
+    }
+
+    size_t m1_rows = m1->m;
+    size_t m1_cols = m1->n;
+    size_t m2_rows = m2->m;
+    size_t m2_cols = m2->n;
+
+    double *m1_values = m1->values;
+    double *m2_values = m2->values;
+    double *result_values = result->values;
+
+    for (int i = 0; i < m1_rows; i++) {
+        for (int j = 0; j < m2_cols; j++) {
+            size_t result_index = m2_cols * i + j;
+            result_values[result_index] = 0.0;
+            for (int k = 0; k < m2_rows; k++) {
+                result_values[result_index] += m1_values[m1_cols * i + k] * m2_values[m2_cols * k + j];
+            }
+        }
+    }
+
+    return 0;
+}
+
+void matrix_destroy(matrix_t *self) {
+    if (self == NULL) return;
+
+    if (self->values != NULL) {
+        free(self->values);
+    }
+
+    free(self);
+}
diff --git a/src/matrix.h b/src/matrix.h
new file mode 100644
index 00000000..93f46f45
--- /dev/null
+++ b/src/matrix.h
@@ -0,0 +1,34 @@
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <stdlib.h>
+#include <stdint.h>
+
+#include "collections.h"
+
+typedef struct matrix {
+    size_t m;
+    size_t n;
+    double *values;
+} matrix_t;
+
+matrix_t *matrix_new(size_t m, size_t n);
+matrix_t *matrix_new_value(size_t m, size_t n, double value);
+matrix_t *matrix_new_values(size_t m, size_t n, double *values);
+matrix_t *matrix_new_zeros(size_t m, size_t n);
+matrix_t *matrix_new_ones(size_t m, size_t n);
+
+matrix_t *matrix_copy(matrix_t *self);
+
+void matrix_init_values(matrix_t *self, double *values);
+void matrix_set(matrix_t *self, double value);
+void matrix_set_row(matrix_t *self, size_t index, double *row);
+void matrix_set_scalar(matrix_t *self, size_t row_index, size_t col_index, double value);
+
+void matrix_scale(matrix_t *self, double value);
+void matrix_dot_vector(matrix_t *self, double *vec, double *result);
+int matrix_dot_matrix(matrix_t *m1, matrix_t *m2, matrix_t *result);
+
+void matrix_destroy(matrix_t *self);
+
+#endif
\ No newline at end of file