diff --git a/Makefile b/Makefile
index c322603..18aa6ee 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,13 @@
 # Compiler and flags
 CC = gcc
-CFLAGS = -I nn/include -I tests -Wall -Wextra -Werror -Wpedantic -Wstrict-prototypes -Wold-style-definition -g $(CU_CFLAGS)
+# OpenMP Flag
+ifdef USE_OPENMP
+    OPENMP_FLAG = -fopenmp
+else
+    OPENMP_FLAG = 
+endif
+
+CFLAGS = -I nn/include -I tests -Wall -Wextra -Werror -Wpedantic -Wstrict-prototypes -Wold-style-definition -g $(CU_CFLAGS) $(OPENMP_FLAG)
 
 # Source files
 SRCS = $(shell find nn/src -name '*.c' -not -path 'nn/src/main.c')
diff --git a/model_summary.txt b/model_summary.txt
new file mode 100644
index 0000000..d978c8b
--- /dev/null
+++ b/model_summary.txt
@@ -0,0 +1,26 @@
+==================================
+      Neural Network Summary      
+==================================
+Number of layers: 2
+----------------------------------
+Layer 1:
+  Weights matrix: 2 x 4
+  Bias matrix:    1 x 4
+  Activation:     RELU
+----------------------------------
+Layer 2:
+  Weights matrix: 4 x 1
+  Bias matrix:    1 x 1
+  Activation:     SIGMOID
+==================================
+
+==================================
+         Model Predictions        
+==================================
+Input -> Expected | Predicted (Rounded)
+----------------------------------
+Input: (0, 0) -> Expected: 0 | Predicted: 0.5000 (Rounded: 1)
+Input: (0, 1) -> Expected: 1 | Predicted: 0.9829 (Rounded: 1)
+Input: (1, 0) -> Expected: 1 | Predicted: 0.9924 (Rounded: 1)
+Input: (1, 1) -> Expected: 0 | Predicted: 0.9989 (Rounded: 1)
+==================================
diff --git a/nn/include/feedforward.h b/nn/include/feedforward.h
index 6e62f65..6488b18 100644
--- a/nn/include/feedforward.h
+++ b/nn/include/feedforward.h
@@ -31,4 +31,4 @@ void free_network(NeuralNetwork* nn);
  * @return Output activation of the last layer (batch_size x output_features).
  *         Caller owns and must free.
  */
-Matrix* feedforward(NeuralNetwork* nn, const Matrix* input);
+Matrix* feedforward(const NeuralNetwork* nn, const Matrix* input);
diff --git a/nn/include/neural_network.h b/nn/include/neural_network.h
index 4092442..91d232f 100644
--- a/nn/include/neural_network.h
+++ b/nn/include/neural_network.h
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <stdio.h>
+
 #include "activation.h"
 #include "cache.h"
 #include "linalg.h"
diff --git a/nn/include/summary.h b/nn/include/summary.h
new file mode 100644
index 0000000..3986b99
--- /dev/null
+++ b/nn/include/summary.h
@@ -0,0 +1,13 @@
+#ifndef NN_SUMMARY_H
+#define NN_SUMMARY_H
+
+#include <stdio.h>
+
+#include "neural_network.h"
+
+void fprint_network_summary(FILE* stream, const NeuralNetwork* nn);
+void flog_training_progress(FILE* stream, int epoch, int epochs, double loss);
+void fprint_model_predictions(FILE* stream, const NeuralNetwork* nn,
+                              const Matrix* x_test, const Matrix* y_test);
+
+#endif  // NN_SUMMARY_H
diff --git a/nn/src/examples/xor.c b/nn/src/examples/xor.c
index 825860a..65bcda4 100644
--- a/nn/src/examples/xor.c
+++ b/nn/src/examples/xor.c
@@ -17,6 +17,7 @@
 #include "linalg.h"
 #include "loss.h"
 #include "neural_network.h"
+#include "summary.h"
 #include "utils.h"
 
 /**
@@ -76,7 +77,16 @@ int main() {
   }
 
   // 3. Print Network Summary
-  print_network_summary(nn);
+  FILE* summary_file = fopen("model_summary.txt", "w");
+  if (summary_file == NULL) {
+    LOG_ERROR("Failed to open model_summary.txt for writing.");
+    // Handle error, but continue for now to avoid stopping the whole process
+  } else {
+    fprint_network_summary(summary_file, nn);
+    fprint_model_predictions(summary_file, nn, x_train,
+                             y_train);  // Add predictions
+    fclose(summary_file);
+  }
 
   // 4. Training Parameters
   double learning_rate = 0.1;
@@ -85,6 +95,11 @@ int main() {
   printf("Training XOR network with %d epochs, learning rate %.2f\n", epochs,
          learning_rate);
 
+  FILE* log_file = fopen("training_log.txt", "w");
+  if (log_file == NULL) {
+    LOG_ERROR("Failed to open training_log.txt for writing.");
+  }
+
   // 5. Training Loop
   for (int epoch = 0; epoch < epochs; epoch++) {
     double total_loss = 0;
@@ -124,7 +139,9 @@ int main() {
     }
     free_matrix(y_hat);
 
-    log_training_progress(epoch, epochs, total_loss);
+    if (log_file != NULL) {
+      flog_training_progress(log_file, epoch, epochs, total_loss);
+    }
   }
 
   printf("\nTraining complete. Testing network...\n");
@@ -153,6 +170,9 @@ int main() {
   free_matrix(x_train);
   free_matrix(y_train);
   free_matrix(predictions);
+
+  if (log_file != NULL) fclose(log_file);
+
   free_network(nn);
 
   printf("\nMemory freed. XOR example finished.\n");
diff --git a/nn/src/linalg/operations.c b/nn/src/linalg/operations.c
index 8c66b71..02986f3 100644
--- a/nn/src/linalg/operations.c
+++ b/nn/src/linalg/operations.c
@@ -15,6 +15,10 @@
 #include "linalg.h"
 #include "utils.h"
 
+#ifdef USE_OPENMP
+#include <omp.h>
+#endif
+
 //============================
 // Functions for Matrix Operations
 //============================
@@ -46,6 +50,9 @@ Matrix* add_matrix(Matrix* a, Matrix* b) {
   Matrix* result = create_matrix(a->rows, a->cols);
   size_t total_elements = a->rows * a->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = a->matrix_data[i] + b->matrix_data[i];
   }
@@ -69,6 +76,9 @@ Matrix* subtract_matrix(Matrix* a, Matrix* b) {
   Matrix* result = create_matrix(a->rows, a->cols);
   size_t total_elements = a->rows * a->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = a->matrix_data[i] - b->matrix_data[i];
   }
@@ -94,6 +104,9 @@ Matrix* multiply_matrix(Matrix* a, Matrix* b) {
   Matrix* result = create_matrix(a->rows, a->cols);
   size_t total_elements = a->rows * a->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = a->matrix_data[i] * b->matrix_data[i];
   }
@@ -115,6 +128,9 @@ Matrix* apply_onto_matrix(double (*func)(double), Matrix* m) {
   Matrix* result = create_matrix(m->rows, m->cols);
   size_t total_elements = m->rows * m->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = func(m->matrix_data[i]);
   }
@@ -135,6 +151,9 @@ Matrix* add_scalar_to_matrix(Matrix* m, double n) {
   Matrix* result = create_matrix(m->rows, m->cols);
   size_t total_elements = m->rows * m->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = m->matrix_data[i] + n;
   }
@@ -149,29 +168,32 @@ Matrix* add_scalar_to_matrix(Matrix* m, double n) {
  * @param b The second matrix.
  * @return A new matrix containing the result of the dot product a * b.
  */
-Matrix* dot_matrix(Matrix* a, Matrix* b) {
-  ASSERT(a != NULL && b != NULL, "Input matrices cannot be NULL.");
-  ASSERT(a->cols == b->rows,
-         "The number of columns in the first matrix must equal the number of "
-         "rows in the second matrix for dot product.");
-
-  LOG_INFO("Performing dot product on a %zux%zu and a %zux%zu matrix.", a->rows,
-           a->cols, b->rows, b->cols);
-  Matrix* result = create_matrix(a->rows, b->cols);
-
-  for (size_t i = 0; i < a->rows; i++) {
-    for (size_t j = 0; j < b->cols; j++) {
-      double sum = 0;
-      for (size_t k = 0; k < a->cols; k++) {
-        sum +=
-            a->matrix_data[i * a->cols + k] * b->matrix_data[k * b->cols + j];
+Matrix* dot_matrix(Matrix* m1, Matrix* m2) {
+  ASSERT(m1 != NULL && m2 != NULL, "Input matrices cannot be NULL.");
+  ASSERT(m1->cols == m2->rows,
+         "Matrices dimensions are incompatible for dot product.");
+
+  LOG_INFO("Performing dot product of %zux%zu and %zux%zu matrices.", m1->rows,
+           m1->cols, m2->rows, m2->cols);
+  Matrix* result = create_matrix(m1->rows, m2->cols);
+  // Initialize result matrix with zeros
+  memset(result->matrix_data, 0, m1->rows * m2->cols * sizeof(double));
+
+#ifdef USE_OPENMP
+#pragma omp parallel for collapse(2)
+#endif
+  for (size_t i = 0; i < m1->rows; i++) {
+    for (size_t j = 0; j < m2->cols; j++) {
+      for (size_t k = 0; k < m1->cols; k++) {
+        result->matrix_data[i * result->cols + j] +=
+            m1->matrix_data[i * m1->cols + k] *
+            m2->matrix_data[k * m2->cols + j];
       }
-      result->matrix_data[i * result->cols + j] = sum;
     }
   }
 
-  LOG_INFO("Matrix dot product complete. Resulting matrix is %zux%zu.",
-           result->rows, result->cols);
+  LOG_INFO("Dot product complete. Resulting matrix is %zux%zu.", result->rows,
+           result->cols);
   return result;
 }
 
@@ -185,6 +207,9 @@ Matrix* transpose_matrix(Matrix* m) {
   LOG_INFO("Transposing a %zux%zu matrix.", m->rows, m->cols);
   Matrix* result = create_matrix(m->cols, m->rows);
 
+#ifdef USE_OPENMP
+#pragma omp parallel for collapse(2)
+#endif
   for (size_t i = 0; i < m->rows; i++) {
     for (size_t j = 0; j < m->cols; j++) {
       result->matrix_data[j * result->cols + i] =
@@ -209,6 +234,9 @@ Matrix* scale_matrix(double n, Matrix* m) {
   Matrix* result = create_matrix(m->rows, m->cols);
   size_t total_elements = m->rows * m->cols;
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t i = 0; i < total_elements; i++) {
     result->matrix_data[i] = m->matrix_data[i] * n;
   }
@@ -233,6 +261,9 @@ Matrix* add_bias_to_matrix(Matrix* m, Matrix* bias) {
   Matrix* result = create_matrix(m->rows, m->cols);
   ASSERT(result != NULL, "Failed to create matrix for bias addition.");
 
+#ifdef USE_OPENMP
+#pragma omp parallel for collapse(2)
+#endif
   for (size_t i = 0; i < m->rows; i++) {
     for (size_t j = 0; j < m->cols; j++) {
       result->matrix_data[i * m->cols + j] =
@@ -255,6 +286,9 @@ Matrix* sum_matrix_columns(Matrix* m) {
   Matrix* result = create_matrix(1, m->cols);
   ASSERT(result != NULL, "Failed to create matrix for column summation.");
 
+#ifdef USE_OPENMP
+#pragma omp parallel for
+#endif
   for (size_t j = 0; j < m->cols; j++) {
     double sum = 0;
     for (size_t i = 0; i < m->rows; i++) {
diff --git a/nn/src/neural_network/feedforward.c b/nn/src/neural_network/feedforward.c
index 7cd105c..5063782 100644
--- a/nn/src/neural_network/feedforward.c
+++ b/nn/src/neural_network/feedforward.c
@@ -81,7 +81,7 @@ void free_network(NeuralNetwork* nn) {
  * @return A new matrix containing the output of the last layer of the network.
  * The caller is responsible for freeing this matrix.
  */
-Matrix* feedforward(NeuralNetwork* nn, const Matrix* input) {
+Matrix* feedforward(const NeuralNetwork* nn, const Matrix* input) {
   ASSERT(nn != NULL, "Neural Network pointer cannot be NULL.");
   ASSERT(input != NULL, "Input matrix cannot be NULL.");
   ASSERT(input->cols == nn->layers[0]->weights->rows,
diff --git a/nn/src/neural_network/summary.c b/nn/src/neural_network/summary.c
index 550b860..2fea2c6 100644
--- a/nn/src/neural_network/summary.c
+++ b/nn/src/neural_network/summary.c
@@ -1,5 +1,7 @@
+#include <math.h>
 #include <stdio.h>
 
+#include "feedforward.h"
 #include "neural_network.h"
 #include "utils.h"
 
@@ -32,4 +34,40 @@ void flog_training_progress(FILE* stream, int epoch, int epochs, double loss) {
   if (epoch % 100 == 0 || epoch == epochs - 1) {
     fprintf(stream, "Epoch %d/%d, Loss: %f\n", epoch, epochs, loss);
   }
+}
+
+void fprint_model_predictions(FILE* stream, const NeuralNetwork* nn,
+                              const Matrix* x_test, const Matrix* y_test) {
+  if (nn == NULL || x_test == NULL || y_test == NULL) {
+    fprintf(stream,
+            "Cannot print predictions: NeuralNetwork, x_test, or y_test is "
+            "NULL.\n");
+    return;
+  }
+
+  Matrix* predictions = feedforward(nn, x_test);
+  if (predictions == NULL) {
+    fprintf(stream, "Failed to generate predictions.\n");
+    return;
+  }
+
+  fprintf(stream, "\n==================================\n");
+  fprintf(stream, "         Model Predictions        \n");
+  fprintf(stream, "==================================\n");
+  fprintf(stream, "Input -> Expected | Predicted (Rounded)\n");
+  fprintf(stream, "----------------------------------\n");
+
+  for (size_t i = 0; i < x_test->rows; i++) {
+    fprintf(stream, "Input: (");
+    for (size_t j = 0; j < x_test->cols; j++) {
+      fprintf(stream, "%.0f%s", x_test->matrix_data[i * x_test->cols + j],
+              (j == x_test->cols - 1) ? "" : ", ");
+    }
+    fprintf(stream, ") -> Expected: %.0f | Predicted: %.4f (Rounded: %.0f)\n",
+            y_test->matrix_data[i], predictions->matrix_data[i],
+            round(predictions->matrix_data[i]));
+  }
+  fprintf(stream, "==================================\n");
+
+  free_matrix(predictions);
 }
\ No newline at end of file
diff --git a/training_log.txt b/training_log.txt
new file mode 100644
index 0000000..7961c69
--- /dev/null
+++ b/training_log.txt
@@ -0,0 +1,21 @@
+Epoch 0/2000, Loss: 0.312058
+Epoch 100/2000, Loss: 0.257256
+Epoch 200/2000, Loss: 0.019960
+Epoch 300/2000, Loss: 0.006621
+Epoch 400/2000, Loss: 0.003768
+Epoch 500/2000, Loss: 0.002665
+Epoch 600/2000, Loss: 0.002038
+Epoch 700/2000, Loss: 0.001627
+Epoch 800/2000, Loss: 0.001347
+Epoch 900/2000, Loss: 0.001140
+Epoch 1000/2000, Loss: 0.000984
+Epoch 1100/2000, Loss: 0.000863
+Epoch 1200/2000, Loss: 0.000772
+Epoch 1300/2000, Loss: 0.000702
+Epoch 1400/2000, Loss: 0.000642
+Epoch 1500/2000, Loss: 0.000591
+Epoch 1600/2000, Loss: 0.000547
+Epoch 1700/2000, Loss: 0.000509
+Epoch 1800/2000, Loss: 0.000475
+Epoch 1900/2000, Loss: 0.000445
+Epoch 1999/2000, Loss: 0.000419