diff --git a/sor_convergence/sor_static_wavefront/Makefile b/sor_convergence/sor_static_wavefront/Makefile
index f772322..4cee729 100755
--- a/sor_convergence/sor_static_wavefront/Makefile
+++ b/sor_convergence/sor_static_wavefront/Makefile
@@ -1,9 +1,13 @@
-CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
-#CC = gcc
+#CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
+CC = gcc
 OMP_FLAG = -fopenmp
+<<<<<<< HEAD
+CFLAGS = -O3 -c ${OMP_FLAG}
+=======
 CFLAGS = -O1 -c ${OMP_FLAG}
 #CFLAGS = -O2 -c ${OMP_FLAG}
 #CFLAGS = -c ${OMP_FLAG}
+>>>>>>> 7046cb8ce94f73f6687fb2a3c9e75bd38a03fb26
 LFLAGS = -lm
 
 .SUFFIXES : .o .c
diff --git a/sor_convergence/sor_static_wavefront/sor_wavefront b/sor_convergence/sor_static_wavefront/sor_wavefront
new file mode 100755
index 0000000..a8e41a1
Binary files /dev/null and b/sor_convergence/sor_static_wavefront/sor_wavefront differ
diff --git a/sor_convergence/sor_static_wavefront/sor_wavefront.o b/sor_convergence/sor_static_wavefront/sor_wavefront.o
new file mode 100644
index 0000000..b714b27
Binary files /dev/null and b/sor_convergence/sor_static_wavefront/sor_wavefront.o differ
diff --git a/sor_original/Makefile b/sor_original/Makefile
new file mode 100755
index 0000000..731f767
--- /dev/null
+++ b/sor_original/Makefile
@@ -0,0 +1,17 @@
+#CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
+CC = gcc
+OMP_FLAG = -fopenmp
+#CFLAGS = -O2 -c ${OMP_FLAG}
+CFLAGS = -O3 -c ${OMP_FLAG}
+LFLAGS = -lm
+
+.SUFFIXES : .o .c
+
+.c.o:
+	${CC} ${CFLAGS} -o $@ $*.c
+
+sor: sor.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+clean:
+	rm *.o sor 
diff --git a/sor_original/sor b/sor_original/sor
new file mode 100755
index 0000000..5f03efc
Binary files /dev/null and b/sor_original/sor differ
diff --git a/sor_original/sor.c b/sor_original/sor.c
new file mode 100644
index 0000000..a3ef1ad
--- /dev/null
+++ b/sor_original/sor.c
@@ -0,0 +1,97 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+// ***  Solution of Laplace's Equation.
+// ***
+// ***  Uxx + Uyy = 0
+// ***  0 <= x <= pi, 0 <= y <= pi
+// ***  U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
+// ***
+// ***  then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
+// ***
+// ***  Should converge with
+// ***                 tol = 0.001 and M = 22  in  60 iterations.
+// ***   and with tol = 0.001 and M = 102 in 200 iterations.
+// ***   and with tol = 0.001 and M = 502 in 980 iterations.
+// *** 
+
+#define N 502
+#define MAX(a,b)  ( ( (a)>(b) ) ? (a) : (b) )
+
+double x[N][N], xnew[N][N], solution[N][N];
+
+double calcerror(double g[][N], int iter);
+
+int main(int argc, char *argv[]){
+	double tol=0.001, h, omega, error;
+        double pi = (double)4.0*atan((double)1.0);
+	int iter=0, i, j;
+
+	h = M_PI/(double)(N-1);
+
+	for(i=0; i<N; i++)
+		x[i][N-1] = sin((double)i*h);
+	
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N-1; j++)
+			x[i][j] = (double)j*h*x[i][N-1];
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N; j++)
+			solution[i][j] = sinh((double)j*h) * sin((double)i*h)/sinh(M_PI);
+	
+	omega = 2.0/(1.0+sin(M_PI/(double)(N-1)));
+   
+	for(i=0; i<N; i++){
+		for(j=0; j<N; j++){
+			xnew[0][j] = x[0][j];
+			xnew[N-1][j] = x[N-1][j];
+		}
+		xnew[i][0] = x[i][0];
+		xnew[i][N-1] = x[i][N-1];
+	}
+
+	error = calcerror(x, iter);
+
+	while(error >= tol){
+
+		for(i=1; i<N-1; i++)
+			for(j=1; j<N-1; j++){
+
+				xnew[i][j] = x[i][j]+0.25*omega*(xnew[i-1][j] + xnew[i][j-1] + x[i+1][j] + x[i][j+1] - (4*x[i][j]));
+			}
+		
+				
+
+		for(i=1; i<N-1; i++)
+			for(j=1; j<N-1; j++)
+				x[i][j] = xnew[i][j];
+
+		iter++;
+
+		if (fmod(iter, 20) == 0)
+		  error = calcerror(x, iter);
+		
+	}
+	
+	printf("Omega = %0.20f\n", omega);
+	printf("Convergence in %d iterations for %dx%d grid with tolerance %f.\n", iter, N, N, tol);
+
+
+	return 0;
+}
+
+double calcerror(double g[][N], int iter){
+
+	int i,j;
+	double error = 0.0;
+
+	for(i=1; i<N-1; i++)
+		for(j=1; j<N-1; j++)
+			error = MAX(error, fabs(solution[i][j] - x[i][j]));
+
+	printf("On iteration %d error= %f\n",iter, error);
+	return error;
+}
diff --git a/sor_original/sor.o b/sor_original/sor.o
new file mode 100644
index 0000000..23ae536
Binary files /dev/null and b/sor_original/sor.o differ
diff --git a/sor_redblack/Makefile b/sor_redblack/Makefile
new file mode 100755
index 0000000..731f767
--- /dev/null
+++ b/sor_redblack/Makefile
@@ -0,0 +1,17 @@
+#CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
+CC = gcc
+OMP_FLAG = -fopenmp
+#CFLAGS = -O2 -c ${OMP_FLAG}
+CFLAGS = -O3 -c ${OMP_FLAG}
+LFLAGS = -lm
+
+.SUFFIXES : .o .c
+
+.c.o:
+	${CC} ${CFLAGS} -o $@ $*.c
+
+sor: sor.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+clean:
+	rm *.o sor 
diff --git a/sor_redblack/sor b/sor_redblack/sor
new file mode 100755
index 0000000..2caa468
Binary files /dev/null and b/sor_redblack/sor differ
diff --git a/sor_redblack/sor.c b/sor_redblack/sor.c
new file mode 100644
index 0000000..4991bd2
--- /dev/null
+++ b/sor_redblack/sor.c
@@ -0,0 +1,122 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+
+// ***  Solution of Laplace's Equation.
+// ***
+// ***  Uxx + Uyy = 0
+// ***  0 <= x <= pi, 0 <= y <= pi
+// ***  U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
+// ***
+// ***  then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
+// ***
+// ***  Should converge with
+// ***                 tol = 0.001 and M = 22  in  60 iterations.
+// ***   and with tol = 0.001 and M = 102 in 200 iterations.
+// ***   and with tol = 0.001 and M = 502 in 980 iterations.
+// *** 
+
+#define N 502
+#define MAX(a,b)  ( ( (a)>(b) ) ? (a) : (b) )
+
+double x[N][N], xnew[N][N], solution[N][N];
+
+double calcerror(double g[][N], int iter);
+
+int main(int argc, char *argv[]){
+	double tol=0.001, h, omega, error;
+    double pi = (double)4.0*atan((double)1.0);
+	int iter=0, i, j;
+	double total_start;
+	double total_time = 0.0;
+	
+	total_start = omp_get_wtime();
+	h = M_PI/(double)(N-1);
+
+	for(i=0; i<N; i++)
+		x[i][N-1] = sin((double)i*h);
+	
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N-1; j++)
+			x[i][j] = (double)j*h*x[i][N-1];
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N; j++)
+			solution[i][j] = sinh((double)j*h) * sin((double)i*h)/sinh(M_PI);
+	
+	omega = 2.0/(1.0+sin(M_PI/(double)(N-1)));
+   
+	for(i=0; i<N; i++){
+		for(j=0; j<N; j++){
+			xnew[0][j] = x[0][j];
+			xnew[N-1][j] = x[N-1][j];
+		}
+		xnew[i][0] = x[i][0];
+		xnew[i][N-1] = x[i][N-1];
+	}
+
+	error = calcerror(x, iter);
+
+	while(error >= tol){
+		//red
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==0){
+					xnew[i][j] = x[i][j]+0.25*omega*(x[i-1][j] + x[i][j-1] + x[i+1][j] + x[i][j+1] - (4*x[i][j]));
+				}
+			}
+		}
+
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==0){
+					x[i][j] = xnew[i][j];
+				}
+			}
+		}
+
+		//black
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==1){
+					xnew[i][j] = x[i][j]+0.25*omega*(x[i-1][j] + x[i][j-1] + x[i+1][j] + x[i][j+1] - (4*x[i][j]));
+				}
+			}
+		}
+
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==1){
+					x[i][j] = xnew[i][j];
+				}
+			}
+		}
+
+		iter++;
+
+		if (fmod(iter, 20) == 0)
+		  error = calcerror(x, iter);
+		
+	}
+	total_time = omp_get_wtime() - total_start;
+	printf("Omega = %0.20f\n", omega);
+	printf("Convergence in %d iterations for %dx%d grid with tolerance %f.\n", iter, N, N, tol);
+	printf("Total time to convergence: %f seconds\n", total_time);
+
+	return 0;
+}
+
+double calcerror(double g[][N], int iter){
+
+	int i,j;
+	double error = 0.0;
+
+	for(i=1; i<N-1; i++)
+		for(j=1; j<N-1; j++)
+			error = MAX(error, fabs(solution[i][j] - x[i][j]));
+
+	printf("On iteration %d error= %f\n",iter, error);
+	return error;
+}
diff --git a/sor_redblack/sor.o b/sor_redblack/sor.o
new file mode 100644
index 0000000..0c52548
Binary files /dev/null and b/sor_redblack/sor.o differ
diff --git a/sor_redblack_OpenMP/Makefile b/sor_redblack_OpenMP/Makefile
new file mode 100755
index 0000000..731f767
--- /dev/null
+++ b/sor_redblack_OpenMP/Makefile
@@ -0,0 +1,17 @@
+#CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
+CC = gcc
+OMP_FLAG = -fopenmp
+#CFLAGS = -O2 -c ${OMP_FLAG}
+CFLAGS = -O3 -c ${OMP_FLAG}
+LFLAGS = -lm
+
+.SUFFIXES : .o .c
+
+.c.o:
+	${CC} ${CFLAGS} -o $@ $*.c
+
+sor: sor.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+clean:
+	rm *.o sor 
diff --git a/sor_redblack_OpenMP/sor b/sor_redblack_OpenMP/sor
new file mode 100755
index 0000000..407ce1e
Binary files /dev/null and b/sor_redblack_OpenMP/sor differ
diff --git a/sor_redblack_OpenMP/sor.c b/sor_redblack_OpenMP/sor.c
new file mode 100644
index 0000000..5ebb1f0
--- /dev/null
+++ b/sor_redblack_OpenMP/sor.c
@@ -0,0 +1,130 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+
+// ***  Solution of Laplace's Equation.
+// ***
+// ***  Uxx + Uyy = 0
+// ***  0 <= x <= pi, 0 <= y <= pi
+// ***  U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
+// ***
+// ***  then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
+// ***
+// ***  Should converge with
+// ***                 tol = 0.001 and M = 22  in  60 iterations.
+// ***   and with tol = 0.001 and M = 102 in 200 iterations.
+// ***   and with tol = 0.001 and M = 502 in 980 iterations.
+// *** 
+
+#define N 180 
+#define MAX(a,b)  ( ( (a)>(b) ) ? (a) : (b) )
+
+
+double x[N][N], xnew[N][N], solution[N][N];
+
+double calcerror(double g[][N], int iter);
+
+int main(int argc, char *argv[]){
+	double tol=0.001, h, omega, error;
+    double pi = (double)4.0*atan((double)1.0);
+	int iter=0, i, j;
+	double total_start;
+	double total_time = 0.0;
+	int numthreads;
+	
+	total_start = omp_get_wtime();
+	h = M_PI/(double)(N-1);
+
+	for(i=0; i<N; i++)
+		x[i][N-1] = sin((double)i*h);
+	
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N-1; j++)
+			x[i][j] = (double)j*h*x[i][N-1];
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N; j++)
+			solution[i][j] = sinh((double)j*h) * sin((double)i*h)/sinh(M_PI);
+	
+	omega = 2.0/(1.0+sin(M_PI/(double)(N-1)));
+   
+	for(i=0; i<N; i++){
+		for(j=0; j<N; j++){
+			xnew[0][j] = x[0][j];
+			xnew[N-1][j] = x[N-1][j];
+		}
+		xnew[i][0] = x[i][0];
+		xnew[i][N-1] = x[i][N-1];
+	}
+
+	error = calcerror(x, iter);
+
+	while(error >= tol){
+		#pragma omp parallel for			\
+  		schedule (static)	
+		//red
+		for(i=1; i<N-1; i++){
+		numthreads = omp_get_num_threads();	
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==0){
+					xnew[i][j] = x[i][j]+0.25*omega*(x[i-1][j] + x[i][j-1] + x[i+1][j] + x[i][j+1] - (4*x[i][j]));
+				}
+			}
+		}
+		#pragma omp parallel for			\
+  		schedule (static)
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==0){
+					x[i][j] = xnew[i][j];
+				}
+			}
+		}
+		#pragma omp parallel for			\
+  		schedule (static)
+		//black
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==1){
+					xnew[i][j] = x[i][j]+0.25*omega*(x[i-1][j] + x[i][j-1] + x[i+1][j] + x[i][j+1] - (4*x[i][j]));
+				}
+			}
+		}
+		#pragma omp parallel for			\
+  		schedule (static)
+		for(i=1; i<N-1; i++){
+			for(j=1; j<N-1; j++){
+				if((i+j)%2==1){
+					x[i][j] = xnew[i][j];
+				}
+			}
+		}
+
+		iter++;
+
+		if (fmod(iter, 20) == 0)
+		  error = calcerror(x, iter);
+		
+	}
+	total_time = omp_get_wtime() - total_start;
+	printf("Omega = %0.20f\n", omega);
+	printf("Convergence in %d iterations for %dx%d grid with tolerance %f.\n", iter, N, N, tol);
+	printf("Total time to convergence: %f seconds\n", total_time);
+	printf("Number of threads: %i\n", numthreads);
+	return 0;
+}
+
+double calcerror(double g[][N], int iter){
+
+	int i,j;
+	double error = 0.0;
+
+	for(i=1; i<N-1; i++)
+		for(j=1; j<N-1; j++)
+			error = MAX(error, fabs(solution[i][j] - x[i][j]));
+
+	printf("On iteration %d error= %f\n",iter, error);
+	return error;
+}
diff --git a/sor_redblack_OpenMP/sor.o b/sor_redblack_OpenMP/sor.o
new file mode 100644
index 0000000..0c503ea
Binary files /dev/null and b/sor_redblack_OpenMP/sor.o differ
diff --git a/sor_redblack_OpenMP/sor.sct b/sor_redblack_OpenMP/sor.sct
new file mode 100755
index 0000000..c03024e
--- /dev/null
+++ b/sor_redblack_OpenMP/sor.sct
@@ -0,0 +1,32 @@
+#!/bin/sh
+#
+# request Bourne shell as shell for job
+#$ -S /bin/sh
+#
+# use current working directory
+#$ -cwd
+#
+# join the output and error output in one file
+#$ -j y
+#
+# use the par environment for parallel jobs
+#$ -pe par 36
+#
+# set up affinity mask (defines cores on which to run threads)
+# default mask
+#export #GOMP_CPU_AFFINITY="0,1,2,4,5,6,8,9,10,12,13,14,16,17,18,20,21,22,24,25,26,28,29,30,32,33,34,36,37,38,40,41,42,44,45,46"
+
+# compact mask (4 cores on 1 proc)
+#export GOMP_CPU_AFFINITY="0,4,8,12"
+# scattered mask (1 core on each of 4 procs)
+#export GOMP_CPU_AFFINITY="0,1,2,3"
+
+#for n in `seq 1 7`
+#do
+for p in `seq 1 36`
+do
+export OMP_NUM_THREADS=$p
+echo Running on $p threads
+./sor
+done
+#done
diff --git a/sor_static_wavefront_redblack/Makefile b/sor_static_wavefront_redblack/Makefile
new file mode 100755
index 0000000..ec93c3c
--- /dev/null
+++ b/sor_static_wavefront_redblack/Makefile
@@ -0,0 +1,27 @@
+#CC = /opt/rh/devtoolset-2/root/usr/bin/gcc
+CC = gcc
+OMP_FLAG = -fopenmp
+#<<<<<<< HEAD
+CFLAGS = -O3 -c ${OMP_FLAG}
+#=======
+#CFLAGS = -O1 -c ${OMP_FLAG}
+#CFLAGS = -O2 -c ${OMP_FLAG}
+#CFLAGS = -c ${OMP_FLAG}
+LFLAGS = -lm
+
+.SUFFIXES : .o .c
+
+.c.o:
+	${CC} ${CFLAGS} -o $@ $*.c
+
+sor: sor.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+sor_wavefront_parallel: sor_wavefront_parallel.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+sor_wavefront_rb_parallel: sor_wavefront_rb_parallel.o
+	${CC} ${OMP_FLAG} -o $@ $@.o ${LFLAGS}
+
+clean:
+	rm *.o sor sor_wavefront_sequential sor_wavefront_parallel sor_wavefront_rb_parallel
diff --git a/sor_static_wavefront_redblack/README b/sor_static_wavefront_redblack/README
new file mode 100644
index 0000000..ef56a03
--- /dev/null
+++ b/sor_static_wavefront_redblack/README
@@ -0,0 +1,4 @@
+better_sor_wavefront_sequential.c finishes within 0.022 seconds
+sor_wavefront_sequential finishes within 1.2 seconds. 
+
+Compare these two and report. 
diff --git a/sor_static_wavefront_redblack/run_sor_par.sct b/sor_static_wavefront_redblack/run_sor_par.sct
new file mode 100755
index 0000000..2f513b3
--- /dev/null
+++ b/sor_static_wavefront_redblack/run_sor_par.sct
@@ -0,0 +1,23 @@
+#!/bin/sh
+#
+# request Bourne shell as shell for job
+#$ -S /bin/sh
+#
+# use current working directory
+#$ -cwd
+#
+# join the output and error output in one file
+#$ -j y
+#
+# use the par environment for parallel jobs
+#$ -pe par 36
+#
+# set up affinity mask (defines cores on which to run threads)
+export GOMP_CPU_AFFINITY="0,1,2,4,5,6,8,9,10,12,13,14,16,17,18,20,21,22,24,25,26,28,29,30,32,33,34,36,37,38,40,41,42,44,45,46"
+
+for p in 1 2 3 4 6 8 12 24 36 
+do
+	export OMP_NUM_THREADS=$p
+	echo "Running on $p threads"
+	./sor_wavefront_parallel
+done
diff --git a/sor_static_wavefront_redblack/sor_wavefront b/sor_static_wavefront_redblack/sor_wavefront
new file mode 100755
index 0000000..a8e41a1
Binary files /dev/null and b/sor_static_wavefront_redblack/sor_wavefront differ
diff --git a/sor_static_wavefront_redblack/sor_wavefront_parallel b/sor_static_wavefront_redblack/sor_wavefront_parallel
new file mode 100755
index 0000000..f151adf
Binary files /dev/null and b/sor_static_wavefront_redblack/sor_wavefront_parallel differ
diff --git a/sor_static_wavefront_redblack/sor_wavefront_parallel.c b/sor_static_wavefront_redblack/sor_wavefront_parallel.c
new file mode 100644
index 0000000..4aaa216
--- /dev/null
+++ b/sor_static_wavefront_redblack/sor_wavefront_parallel.c
@@ -0,0 +1,145 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+
+// ***  Solution of Laplace's Equation.
+// ***   Note that this experiment is using wavefront
+// ***  Uxx + Uyy = 0
+// ***  0 <= x <= pi, 0 <= y <= pi
+// ***  U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
+// ***
+// ***  then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
+// ***
+// ***  Should converge with
+// ***   tol = 0.001 and N = 22  in  60 iterations.
+// ***   and with tol = 0.001 and N = 102 in 200 iterations.
+// ***   and with tol = 0.001 and N = 502 in 980 iterations.
+// *** 
+
+#define N 502 // 50 is a good size to be used for debugging. The wavefront algorithm is VERY slow on sequential machines. 
+#define MAX(a,b)  ( ( (a)>(b) ) ? (a) : (b) )
+
+// To paralise the code, x[][], xnew[][] and solution[][] should not be global
+//double x[N][N], xnew[N][N], solution[N][N];
+
+double calcerror(double **g, int iter, double **s);
+void matrix_initialise(double **x_matrix, double **xnew_matrix, double **solution_matrix, double h);
+
+int main(int argc, char *argv[]){
+	double tol=0.001, h, omega, error;
+	double pi = (double)4.0*atan((double)1.0);
+	int iter=0, i, j;
+	double x[N][N], **x_ptr; 
+	double xnew[N][N], **xnew_ptr;
+	double solution[N][N], **solution_ptr;
+	int slength; // variable to skip the first and last two diagonal strips
+	int z; // variable to skip boundary elements
+
+	// variables for omp timmer
+	double calcerror_start; 
+	double calcerror_time = 0.0;
+
+	// calculate constant values, h and omega 
+	h = M_PI/(double)(N-1);
+	omega = 2.0/(1.0+sin(M_PI/(double)(N-1)));
+
+	// set up pointers to allow x[][], xnew[][] and solution[][] access in calcerror
+	x_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		x_ptr[i] = x[i];
+	}
+
+	xnew_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		xnew_ptr[i] = xnew[i];
+	}
+
+	solution_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		solution_ptr[i] = solution[i];
+	}
+   
+	// initialise x, xnew, y
+	matrix_initialise(x_ptr, xnew_ptr, solution_ptr, h);
+
+	// boundary conditions, copy the box boundaires from x to xnew array
+	for(i=0; i<N; i++){
+		for(j=0; j<N; j++){
+			xnew[0][j] = x[0][j];
+			xnew[N-1][j] = x[N-1][j];
+		}
+		xnew[i][0] = x[i][0];
+		xnew[i][N-1] = x[i][N-1];
+	}
+
+	// start time
+    calcerror_start = omp_get_wtime();
+
+	// calculate the initial error
+	error = calcerror(x_ptr, iter, solution_ptr);
+
+	while(error >= tol){
+		// wavefront parallelism
+		for (i=0; i < (2*N -1); i++){
+			int k = (i < N) ? 0 : (i - N + 1);
+			slength = (i - k -k + 1);
+			if (slength > 2) {
+				#pragma omp parallel for
+				for (j=k; j <= i-k; j++){
+					if (j==k || j == i-k) {
+						continue; // skip boundary conditions
+					}
+					xnew[j][i-j] = x[j][i-j]+0.25*omega*(xnew[j-1][i-j] + xnew[j][i-j-1] \
+									+ x[j+1][i-j] + x[j][i-j+1] - (4*x[j][i-j]));
+				}
+			}
+		}
+		#pragma omp parallel for
+		for(i=1; i<N-1; i++)
+			for(j=1; j<N-1; j++)
+				x[i][j] = xnew[i][j];
+
+		iter++;
+
+		if (fmod(iter, 20) == 0)
+			error = calcerror(x_ptr, iter, solution_ptr);
+	}
+	printf("Omega = %0.20f\n", omega);
+	printf("Convergence in %d iterations for %dx%d grid with tolerance %f.\n", iter, N, N, tol);
+
+	calcerror_time = calcerror_time + (omp_get_wtime() - calcerror_start);
+
+	printf("Total iteration to converge = %d \n", iter);
+	printf("Total time elapsed to converge = %f \n", calcerror_time);
+
+	return 0;
+}
+
+void matrix_initialise(double **x_matrix, double **xnew_matrix, double **solution_matrix, double h){
+	int i, j;
+	for(i=0; i<N; i++)
+		x_matrix[i][N-1] = sin((double)i*h);
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N-1; j++)
+			x_matrix[i][j] = (double)j*h*x_matrix[i][N-1];
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N; j++)
+			solution_matrix[i][j] = sinh((double)j*h) * sin((double)i*h)/sinh(M_PI);
+}
+
+double calcerror(double **g, int iter, double **s){
+	
+	// TO-DO: openmp reduce
+	int i,j;
+	double error = 0.0;
+
+	for(i=1; i<N-1; i++)
+		for(j=1; j<N-1; j++)
+			error = MAX(error, fabs(s[i][j] - g[i][j]));
+
+	printf("On iteration %d error= %f\n",iter, error);
+	return error;
+}
diff --git a/sor_static_wavefront_redblack/sor_wavefront_parallel.o b/sor_static_wavefront_redblack/sor_wavefront_parallel.o
new file mode 100644
index 0000000..7f5dd0d
Binary files /dev/null and b/sor_static_wavefront_redblack/sor_wavefront_parallel.o differ
diff --git a/sor_static_wavefront_redblack/sor_wavefront_rb_parallel b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel
new file mode 100755
index 0000000..b26b198
Binary files /dev/null and b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel differ
diff --git a/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.c b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.c
new file mode 100644
index 0000000..4425d6f
--- /dev/null
+++ b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.c
@@ -0,0 +1,159 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+
+// ***  Solution of Laplace's Equation.
+// ***   Note that this experiment is using wavefront
+// ***  Uxx + Uyy = 0
+// ***  0 <= x <= pi, 0 <= y <= pi
+// ***  U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
+// ***
+// ***  then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
+// ***
+// ***  Should converge with
+// ***   tol = 0.001 and N = 22  in  60 iterations.
+// ***   and with tol = 0.001 and N = 102 in 200 iterations.
+// ***   and with tol = 0.001 and N = 502 in 980 iterations.
+// *** 
+
+#define N 180 // 50 is a good size to be used for debugging. The wavefront algorithm is VERY slow on sequential machines. 
+#define MAX(a,b)  ( ( (a)>(b) ) ? (a) : (b) )
+
+// To paralise the code, x[][], xnew[][] and solution[][] should not be global
+//double x[N][N], xnew[N][N], solution[N][N];
+
+double calcerror(double **g, int iter, double **s);
+void matrix_initialise(double **x_matrix, double **xnew_matrix, double **solution_matrix, double h);
+
+int main(int argc, char *argv[]){
+	double tol=0.001, h, omega, error;
+	double pi = (double)4.0*atan((double)1.0);
+	int iter=0, i, j;
+	double x[N][N], **x_ptr; 
+	double xnew[N][N], **xnew_ptr;
+	double solution[N][N], **solution_ptr;
+	int slength; // variable to skip the first and last two diagonal strips
+	int z; // variable to skip boundary elements
+
+	// variables for omp timmer
+	double calcerror_start; 
+	double calcerror_time = 0.0;
+
+	// calculate constant values, h and omega 
+	h = M_PI/(double)(N-1);
+	omega = 2.0/(1.0+sin(M_PI/(double)(N-1)));
+
+	// set up pointers to allow x[][], xnew[][] and solution[][] access in calcerror
+	x_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		x_ptr[i] = x[i];
+	}
+
+	xnew_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		xnew_ptr[i] = xnew[i];
+	}
+
+	solution_ptr = (double **)malloc((unsigned) N * sizeof(double));
+	for (i=0; i < N; i++) {
+		solution_ptr[i] = solution[i];
+	}
+   
+	// initialise x, xnew, y
+	matrix_initialise(x_ptr, xnew_ptr, solution_ptr, h);
+
+	// boundary conditions, copy the box boundaires from x to xnew array
+	for(i=0; i<N; i++){
+		for(j=0; j<N; j++){
+			xnew[0][j] = x[0][j];
+			xnew[N-1][j] = x[N-1][j];
+		}
+		xnew[i][0] = x[i][0];
+		xnew[i][N-1] = x[i][N-1];
+	}
+
+	// start time
+    calcerror_start = omp_get_wtime();
+
+	// calculate the initial error
+	error = calcerror(x_ptr, iter, solution_ptr);
+
+	while(error >= tol){
+		// wavefront parallelism
+		for (i=0; i < (2*N -1); i++){
+			int k = (i < N) ? 0 : (i - N + 1);
+			slength = (i - k -k + 1);
+			if (slength > 2 && ((slength+2)%2==0)) {
+				#pragma omp parallel for
+				for (j=k; j <= i-k; j++){
+					if (j==k || j == i-k) {
+						continue; // skip boundary conditions
+					}
+					xnew[j][i-j] = x[j][i-j]+0.25*omega*(xnew[j-1][i-j] + xnew[j][i-j-1] \
+									+ x[j+1][i-j] + x[j][i-j+1] - (4*x[j][i-j]));
+				}
+				for(i=1; i<N-1; i++)
+					for(j=1; j<N-1; j++)
+						x[i][j] = xnew[i][j];
+			}
+			if (slength > 2 && ((slength+2)%2==1)) {
+				#pragma omp parallel for
+				for (j=k; j <= i-k; j++){
+					if (j==k || j == i-k) {
+						continue; // skip boundary conditions
+					}
+					xnew[j][i-j] = x[j][i-j]+0.25*omega*(xnew[j-1][i-j] + xnew[j][i-j-1] \
+									+ x[j+1][i-j] + x[j][i-j+1] - (4*x[j][i-j]));
+				}
+				for(i=1; i<N-1; i++)
+					for(j=1; j<N-1; j++)
+						x[i][j] = xnew[i][j];
+			}
+		}
+
+	
+
+		iter++;
+
+		if (fmod(iter, 20) == 0)
+			error = calcerror(x_ptr, iter, solution_ptr);
+	}
+	printf("Omega = %0.20f\n", omega);
+	printf("Convergence in %d iterations for %dx%d grid with tolerance %f.\n", iter, N, N, tol);
+
+	calcerror_time = calcerror_time + (omp_get_wtime() - calcerror_start);
+
+	printf("Total iteration to converge = %d \n", iter);
+	printf("Total time elapsed to converge = %f \n", calcerror_time);
+
+	return 0;
+}
+
+void matrix_initialise(double **x_matrix, double **xnew_matrix, double **solution_matrix, double h){
+	int i, j;
+	for(i=0; i<N; i++)
+		x_matrix[i][N-1] = sin((double)i*h);
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N-1; j++)
+			x_matrix[i][j] = (double)j*h*x_matrix[i][N-1];
+
+	for(i=0; i<N; i++)
+		for(j=0; j<N; j++)
+			solution_matrix[i][j] = sinh((double)j*h) * sin((double)i*h)/sinh(M_PI);
+}
+
+double calcerror(double **g, int iter, double **s){
+	
+	// TO-DO: openmp reduce
+	int i,j;
+	double error = 0.0;
+
+	for(i=1; i<N-1; i++)
+		for(j=1; j<N-1; j++)
+			error = MAX(error, fabs(s[i][j] - g[i][j]));
+
+	printf("On iteration %d error= %f\n",iter, error);
+	return error;
+}
diff --git a/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.o b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.o
new file mode 100644
index 0000000..1c680d6
Binary files /dev/null and b/sor_static_wavefront_redblack/sor_wavefront_rb_parallel.o differ