diff --git a/nn/include/linalg.h b/nn/include/linalg.h
index 1d1f6c6..078edd1 100644
--- a/nn/include/linalg.h
+++ b/nn/include/linalg.h
@@ -18,16 +18,15 @@ typedef struct _Matrix {
 // Functions for Matrix IO
 //=====================
 
-Matrix* read_matrix(char* fileName);
+Matrix* read_matrix(const char* filename);
 Matrix* create_matrix(int rows, int cols);
 Matrix* copy_matrix(const Matrix* m);
 Matrix* flatten_matrix(Matrix* m, int axis);
 void fill_matrix(Matrix* m, double n);
-void randomize_matrix(Matrix* m);
+void randomize_matrix(Matrix* m, double n);
 void free_matrix(Matrix* m);
-void write_matrix(Matrix* m, char* filename);
+void write_matrix(Matrix* m, const char* filename);
 void print_matrix(Matrix* m);
-void save_matrix(Matrix* m);
 int matrix_argmax(Matrix* m);
 
 //============================
@@ -37,11 +36,8 @@ Matrix* identity_matrix(int n);
 Matrix* add_matrix(Matrix* a, Matrix* b);
 Matrix* subtract_matrix(Matrix* a, Matrix* b);
 Matrix* multiply_matrix(Matrix* a, Matrix* b);
-Matrix* apply_onto_matrix(double* (*func)(double),
-                          Matrix* m);  // Apply our own stuff onto the matrix
-                                       // items, useful for activation
+Matrix* apply_onto_matrix(double* (*func)(double), Matrix* m);
 Matrix* add_scalar_to_matrix(Matrix* m, double n);
 Matrix* dot_matrix(Matrix* a, Matrix* b);
 Matrix* transpose_matrix(Matrix* m);
-Matrix* scale_matrix(double n,
-                     Matrix* m);  // matrix can be scaled to non integer values
+Matrix* scale_matrix(double n, Matrix* m);
diff --git a/nn/src/linalg/io.c b/nn/src/linalg/io.c
new file mode 100644
index 0000000..258f881
--- /dev/null
+++ b/nn/src/linalg/io.c
@@ -0,0 +1,209 @@
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "linalg.h"
+#include "utils.h"
+
+Matrix* read_matrix(const char* filename) {
+  LOG_INFO("Attempting to load matrix from file: %s", filename);
+
+  FILE* file = fopen(filename, "r");
+  ASSERT(file != NULL, "Failed to open file for matrix loading.");
+
+  char entry[1024];
+  int rows = 0, cols = 0;
+
+  if (fgets(entry, sizeof(entry), file) == NULL) {
+    LOG_ERROR("Could not read rows from file: %s", filename);
+    fclose(file);
+    return NULL;
+  }
+  rows = atoi(entry);
+
+  if (fgets(entry, sizeof(entry), file) == NULL) {
+    LOG_ERROR("Could not read columns from file: %s", filename);
+    fclose(file);
+    return NULL;
+  }
+  cols = atoi(entry);
+
+  ASSERT(rows > 0 && cols > 0, "Invalid matrix dimensions read from file.");
+
+  Matrix* m = create_matrix(rows, cols);
+
+  for (int i = 0; i < rows; i++) {
+    if (fgets(entry, sizeof(entry), file) == NULL) {
+      LOG_ERROR("Unexpected end of file while reading matrix data.");
+      free_matrix(m);
+      fclose(file);
+      return NULL;
+    }
+
+    char* line_ptr = entry;
+    for (int j = 0; j < cols; j++) {
+      m->matrix_data[i * cols + j] = strtod(line_ptr, &line_ptr);
+    }
+  }
+
+  LOG_INFO("Successfully loaded a %dx%d matrix from %s.", m->rows, m->cols,
+           filename);
+  fclose(file);
+  return m;
+}
+
+Matrix* create_matrix(int rows, int cols) {
+  LOG_INFO("Creating a new matrix of size %dx%d.", rows, cols);
+
+  Matrix* matrix = (Matrix*)malloc(sizeof(Matrix));
+  CHECK_MALLOC(matrix, "Failed to allocate memory for Matrix struct.");
+
+  matrix->matrix_data = (double*)malloc(rows * cols * sizeof(double));
+  CHECK_MALLOC(matrix->matrix_data,
+               "Failed to allocate memory for matrix data.");
+
+  matrix->rows = rows;
+  matrix->cols = cols;
+
+  LOG_INFO("Matrix created successfully at address %p.", matrix);
+
+  return matrix;
+}
+
+Matrix* copy_matrix(const Matrix* m) {
+  ASSERT(m != NULL, "Input matrix for copy is NULL.");
+
+  LOG_INFO("Copying a %dx%d matrix.", m->rows, m->cols);
+  Matrix* new_matrix = create_matrix(m->rows, m->cols);
+
+  size_t total_bytes = m->rows * m->cols * sizeof(double);
+  memcpy(new_matrix->matrix_data, m->matrix_data, total_bytes);
+
+  LOG_INFO("Matrix copy complete.");
+
+  return new_matrix;
+}
+
+// Helper function for flattening
+Matrix* flatten_column_wise(const Matrix* m) {
+  Matrix* new_matrix = create_matrix(m->rows * m->cols, 1);
+
+  int k = 0;
+  for (int j = 0; j < m->cols; j++) {
+    for (int i = 0; i < m->rows; i++) {
+      new_matrix->matrix_data[k] = m->matrix_data[i * m->cols + j];
+      k++;
+    }
+  }
+  return new_matrix;
+}
+
+Matrix* flatten_matrix(Matrix* m, int axis) {
+  ASSERT(m != NULL, "Input matrix to flatten is NULL.");
+  ASSERT(axis == 0 || axis == 1,
+         "Axis must be 0 (row-wise) or 1 (column-wise).");
+
+  if (axis == 0) {
+    LOG_INFO(
+        "Flattening matrix row-wise. No operation needed as data is "
+        "already "
+        "contiguous.");
+    m->cols = m->rows * m->cols;
+    m->rows = 1;
+    return m;
+  } else {
+    LOG_INFO("Flattening matrix column-wise. A new matrix will be created.");
+    return flatten_column_wise(m);
+  }
+}
+
+void fill_matrix(Matrix* m, double n) {
+  ASSERT(m != NULL, "Input matrix for fill_matrix is NULL.");
+  LOG_INFO("Filling a %dx%d matrix with the value %.2f.", m->rows, m->cols, n);
+
+  for (int i = 0; i < m->rows; i++) {
+    for (int j = 0; j < m->cols; j++) {
+      m->matrix_data[i * m->cols + j] = n;
+    }
+  }
+}
+
+void randomize_matrix(Matrix* m, double n) {
+  LOG_INFO("Randomizing a %dx%d matrix.", m->rows, m->cols);
+  // Apparently a 1/n or 1/n^2 scaling leads to a vanishing gradient problem
+  double min = -1.0 / sqrt(n);
+  double max = 1.0 / sqrt(n);
+  double range = max - min;
+
+  for (int i = 0; i < m->rows; i++) {
+    for (int j = 0; j < m->cols; j++) {
+      double random_value = (double)rand() / (double)RAND_MAX;
+      m->matrix_data[i * m->cols + j] = min + random_value * range;
+    }
+  }
+  LOG_INFO("Matrix randomized successfully.");
+}
+
+void free_matrix(Matrix* m) {
+  LOG_INFO("Freeing matrix at address %p.", m);
+  if (m == NULL) {
+    LOG_WARN("Attempted to free a NULL pointer.");
+    return;
+  }
+  if (m->matrix_data != NULL) {
+    free(m->matrix_data);
+  }
+  free(m);
+  LOG_INFO("Matrix freed successfully.");
+}
+
+void print_matrix(Matrix* m) {
+  ASSERT(m != NULL, "Input matrix for print is NULL.");
+  LOG_INFO("Printing matrix of size %dx%d.", m->rows, m->cols);
+  for (int i = 0; i < m->rows; i++) {
+    for (int j = 0; j < m->cols; j++) {
+      printf("%.3f ", m->matrix_data[i * m->cols + j]);
+    }
+    printf("\n");
+  }
+}
+
+void write_matrix(Matrix* m, const char* filename) {
+  ASSERT(m != NULL, "Input matrix for save is NULL.");
+  LOG_INFO("Saving a %dx%d matrix to file: %s", m->rows, m->cols, filename);
+
+  FILE* file = fopen(filename, "w");
+  ASSERT(file != NULL, "Failed to open file for saving matrix.");
+
+  fprintf(file, "%d\n", m->rows);
+  fprintf(file, "%d\n", m->cols);
+
+  for (int i = 0; i < m->rows; i++) {
+    for (int j = 0; j < m->cols; j++) {
+      fprintf(file, "%.3f ", m->matrix_data[i * m->cols + j]);
+    }
+    fprintf(file, "\n");
+  }
+
+  fclose(file);
+  LOG_INFO("Matrix saved successfully.");
+}
+
+int matrix_argmax(Matrix* m) {
+  ASSERT(m != NULL, "Input matrix for argmax is NULL.");
+  ASSERT(m->cols == 1, "Input must be a column vector (Mx1).");
+
+  double maxValue = INT_MIN;
+  int maxIndex = 0;
+
+  for (int i = 0; i < m->rows; i++) {
+    if (m->matrix_data[i] > maxValue) {
+      maxIndex = i;
+      maxValue = m->matrix_data[i];
+    }
+  }
+  LOG_INFO("Max value found at index %d.", maxIndex);
+  return maxIndex;
+}
diff --git a/nn/src/linalg/operations.c b/nn/src/linalg/operations.c
new file mode 100644
index 0000000..e69de29