diff --git a/Cargo.toml b/Cargo.toml
index 4e02e20..0eced3e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -11,5 +11,6 @@ edition = "2018"
 
 [dependencies]
 rand = "0.8.3"
+rayon = "1.5.0"
 hwloc2 = "2.2.0"
 lazy_static = "1.4.0"
diff --git a/README.md b/README.md
index f45c840..07083a8 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,8 @@ library written in Rust. For more details, please see the project report.
 
 RustMP requires the following dependencies to run:
 - `rand-v0.8.3`: Random number generation used for benchmarking examples
+- `rayon-v1.5.0`: Rayon parallel iterator library used in RustMP/Rayon..
+comparision tests
 - `hwloc2-v2.2.0`: C HWLoc wrapper used for NUMA aware process pinning
 - `lazy_static-v1.4.0`: Lazy static macro used for delayed singleton init
 
diff --git a/omp/Makefile b/omp/Makefile
new file mode 100644
index 0000000..57373c9
--- /dev/null
+++ b/omp/Makefile
@@ -0,0 +1,55 @@
+#
+# File: Makefile
+# Author: Jack Yu (yyu57)
+# Simple generic makefile for my C projects
+# Run 'make build' to compile binary
+#
+
+CC ?= gcc
+CFLAGS := -std=gnu17 -Werror -Wall -O3
+CFLAGS_SEQ := $(CFLAGS)
+CFLAGS_OMP := $(CFLAGS) -fopenmp -DPAR
+CFLAGS_OMP_SAFE := $(CFLAGS_OMP) -DCRIT
+
+BIN := matmul
+BIN_SEQ := seqc_$(BIN)
+BIN_OMP := omp_$(BIN)
+BIN_OMP_SAFE := ompsafe_$(BIN)
+SRCDIR := src
+BLDDIR := build
+SRCS := $(shell find $(SRCDIR) -name '*.c')
+
+# === Recipes ===
+.DEFAULT_GOAL := all
+
+all: $(BIN_SEQ) $(BIN_OMP) $(BIN_OMP_SAFE)
+
+$(BIN_SEQ): $(BLDDIR)/$(BIN_SEQ)
+	@cp $(BLDDIR)/$(BIN_SEQ) $(BIN_SEQ)
+
+$(BIN_OMP): $(BLDDIR)/$(BIN_OMP)
+	@cp $(BLDDIR)/$(BIN_OMP) $(BIN_OMP)
+
+$(BIN_OMP_SAFE): $(BLDDIR)/$(BIN_OMP_SAFE)
+	@cp $(BLDDIR)/$(BIN_OMP_SAFE) $(BIN_OMP_SAFE)
+
+$(BLDDIR)/$(BIN_SEQ): $(BLDDIR) $(SRCS)
+	$(CC) $(CFLAGS_SEQ) $(SRCS) -o $@
+
+$(BLDDIR)/$(BIN_OMP): $(BLDDIR) $(SRCS)
+	$(CC) $(CFLAGS_OMP) $(SRCS) -o $@
+
+$(BLDDIR)/$(BIN_OMP_SAFE): $(BLDDIR) $(SRCS)
+	$(CC) $(CFLAGS_OMP_SAFE) $(SRCS) -o $@
+
+build: all
+
+$(BLDDIR):
+	@mkdir -p $(BLDDIR)
+
+clean:
+	@rm -rf $(BLDDIR)
+	@rm -f $(BIN_SEQ) $(BIN_OMP) $(BIN_OMP_SAFE)
+
+.PHONY: build run clean all
+
diff --git a/omp/src/main.c b/omp/src/main.c
new file mode 100644
index 0000000..a8e5487
--- /dev/null
+++ b/omp/src/main.c
@@ -0,0 +1,100 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/sysinfo.h>
+#include <time.h>
+#include <limits.h>
+#include <unistd.h>
+
+void
+print_time_diff (struct timespec *start, struct timespec *end)
+{
+  time_t sec;
+  long nsec;
+  sec = end->tv_sec - start->tv_sec;
+  nsec = end->tv_nsec - start->tv_nsec;
+
+  if (nsec < 0)
+    {
+      sec--;
+      nsec += 1000000000;
+    }
+
+  printf ("Elapsed time: %ld.%09ld sec\n", sec, nsec);
+}
+
+void
+generate_matrices (double ***matrix, double ***result, size_t nsize)
+{
+    size_t i, j;
+    *matrix = (double **) malloc (nsize * sizeof (double *));
+    *result = (double **) malloc (nsize * sizeof (double *));
+    srand (time (NULL));
+    for (i = 0; i < nsize; i++) {
+        (*matrix)[i] = (double *) malloc (nsize * sizeof (double));
+        (*result)[i] = (double *) malloc (nsize * sizeof (double));
+        for (j = 0; j < nsize; j++) {
+            (*matrix)[i][j] = 255.0f * ((double) rand () - ((double) RAND_MAX/2.0f));
+            (*result)[i][j] = 0;
+        }
+    }
+}
+
+void
+free_all (double **matrix, double **result, size_t nsize)
+{
+    size_t i;
+    for (i = 0; i < nsize; i++) {
+        free (matrix[nsize]);
+        free (result[nsize]);
+    }
+    free (matrix);
+    free (result);
+}
+
+
+int main (int argc, char *argv[])
+{
+    struct timespec start, end;
+    size_t i, j, k, nsize;
+    long raw_size;
+    double sum;
+    double **matrix, **result;
+    if (argc != 2) {
+        printf ("Usage: %s <nsize>\n", argv[0]);
+        exit(-1);
+    }
+    raw_size = atol (argv[1]);
+    if (raw_size < 1) {
+        printf ("Usage: %s <nsize>\n", argv[0]);
+        exit(-1);
+    }
+    nsize = raw_size;
+    generate_matrices (&matrix, &result, nsize);
+    clock_gettime (CLOCK_MONOTONIC, &start);
+#ifdef PAR
+    #pragma omp parallel for default(shared) private(j, k, sum)
+#endif
+    for (i = 0; i < nsize; i++) {
+        for (j = 0; j < nsize; j++) {
+            sum = 0;
+            for (k = 0; k < nsize; k++) {
+                sum += matrix[i][k] * matrix[k][j];
+            }
+#ifdef PAR
+#ifdef CRIT
+            #pragma omp critical
+            {
+#endif
+#endif
+                result[i][j] = sum;
+#ifdef PAR
+#ifdef CRIT
+            }
+#endif
+#endif
+        }
+    }
+    clock_gettime (CLOCK_MONOTONIC, &end);
+    print_time_diff (&start, &end);
+    free_all (matrix, result, nsize);
+}
diff --git a/src/bin/rayon_matmul.rs b/src/bin/rayon_matmul.rs
new file mode 100644
index 0000000..7968e8d
--- /dev/null
+++ b/src/bin/rayon_matmul.rs
@@ -0,0 +1,41 @@
+use std::env;
+use std::cmp::max;
+use rand::random;
+use std::time::Instant;
+use rayon::prelude::*;
+
+fn gen_matrix(nsize: usize) -> Vec<Vec<f64>> {
+    let mut ret = Vec::with_capacity(nsize);
+    for _ in 0..nsize {
+        let mut row = Vec::with_capacity(nsize);
+        for _ in 0..nsize {
+            row.push((random::<f64>() - 0.5) * 255.0);
+        }
+        ret.push(row);
+    }
+    ret
+}
+
+fn main() {
+    let args: Vec<String> = env::args().collect();
+    if args.len() != 2 {
+        eprintln!("Usage: {} <msize>", args[0]);
+        return;
+    }
+    let nsize = max(args[1].parse::<usize>().expect("Usage: matrix_mul <msize>"), 1);
+    let matrix = gen_matrix(nsize);
+    let timer = Instant::now();
+    let _result = (0..nsize).into_par_iter().map(|i| {
+        let mut res_row = Vec::with_capacity(nsize);
+        for j in 0..nsize {
+            let mut sum: f64 = 0.0;
+            for k in 0..nsize {
+                sum += matrix[i][k] * matrix[k][j];
+            }
+            res_row.push(sum);
+        }
+        res_row
+    }).collect::<Vec<Vec<f64>>>();
+    let interval = timer.elapsed();
+    println!("Elapsed time: {:?}", interval);
+}
diff --git a/src/bin/rmp_matmul.rs b/src/bin/rmp_matmul.rs
new file mode 100644
index 0000000..a927868
--- /dev/null
+++ b/src/bin/rmp_matmul.rs
@@ -0,0 +1,58 @@
+use std::env;
+use std::cmp::max;
+use rand::random;
+use std::time::Instant;
+
+use rustmp::{par_for, critical};
+
+fn gen_matrix(nsize: usize) -> Vec<Vec<f64>> {
+    let mut ret = Vec::with_capacity(nsize);
+    for _ in 0..nsize {
+        let mut row = Vec::with_capacity(nsize);
+        for _ in 0..nsize {
+            row.push((random::<f64>() - 0.5) * 255.0);
+        }
+        ret.push(row);
+    }
+    ret
+}
+
+fn gen_empty(nsize: usize) -> Vec<Vec<f64>> {
+    let mut ret = Vec::with_capacity(nsize);
+    let mut row = Vec::with_capacity(nsize);
+    for _ in 0..nsize {
+        row.push(0 as f64);
+    }
+    for _ in 0..nsize {
+        ret.push(row.clone());
+    }
+    ret
+}
+
+fn main() {
+    let args: Vec<String> = env::args().collect();
+    if args.len() != 2 {
+        eprintln!("Usage: {} <msize>", args[0]);
+        return;
+    }
+    let nsize = max(args[1].parse::<usize>().expect("Usage: matrix_mul <msize>"), 1);
+    let matrix = gen_matrix(nsize);
+    let result = gen_empty(nsize);
+    let timer = Instant::now();
+    par_for! {
+        for i in 0..nsize, locked result, read matrix, {
+            for j in 0..nsize {
+                let mut sum = 0.0;
+                for k in 0..nsize {
+                    sum += matrix[i][k] * matrix[k][j];
+                }
+                critical! {
+                    readwrite result;
+                    result[i][j] = sum;
+                }
+            }
+        }
+    }
+    let interval = timer.elapsed();
+    println!("Elapsed time: {:?}", interval);
+}
diff --git a/src/bin/seqrs_matmul.rs b/src/bin/seqrs_matmul.rs
new file mode 100644
index 0000000..4c6f975
--- /dev/null
+++ b/src/bin/seqrs_matmul.rs
@@ -0,0 +1,51 @@
+use std::env;
+use std::cmp::max;
+use rand::random;
+use std::time::Instant;
+
+fn gen_matrix(nsize: usize) -> Vec<Vec<f64>> {
+    let mut ret = Vec::with_capacity(nsize);
+    for _ in 0..nsize {
+        let mut row = Vec::with_capacity(nsize);
+        for _ in 0..nsize {
+            row.push((random::<f64>() - 0.5) * 255.0);
+        }
+        ret.push(row);
+    }
+    ret
+}
+
+fn gen_empty(nsize: usize) -> Vec<Vec<f64>> {
+    let mut ret = Vec::with_capacity(nsize);
+    let mut row = Vec::with_capacity(nsize);
+    for _ in 0..nsize {
+        row.push(0 as f64);
+    }
+    for _ in 0..nsize {
+        ret.push(row.clone());
+    }
+    ret
+}
+
+fn main() {
+    let args: Vec<String> = env::args().collect();
+    if args.len() != 2 {
+        eprintln!("Usage: {} <msize>", args[0]);
+        return;
+    }
+    let nsize = max(args[1].parse::<usize>().expect("Usage: matrix_mul <msize>"), 1);
+    let matrix = gen_matrix(nsize);
+    let mut result = gen_empty(nsize);
+    let timer = Instant::now();
+    for i in 0..nsize {
+        for j in 0..nsize {
+            let mut sum = 0.0;
+            for k in 0..nsize {
+                sum += matrix[i][k] * matrix[k][j];
+            }
+            result[i][j] = sum;
+        }
+    }
+    let interval = timer.elapsed();
+    println!("Elapsed time: {:?}", interval);
+}