From 6e2469414bd9247ef0493ab4440edb7ba306b260 Mon Sep 17 00:00:00 2001
From: Tom Deakin <tom.deakin@bristol.ac.uk>
Date: Thu, 12 Mar 2020 15:45:48 +0000
Subject: [PATCH] Add host OpenMP version

---
 heat.c | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/heat.c b/heat.c
index 406ed1d..4117a6c 100644
--- a/heat.c
+++ b/heat.c
@@ -159,6 +159,7 @@ int main(int argc, char *argv[]) {
   printf("Error (L2norm): %E\n", norm);
   printf("Solve time (s): %lf\n", toc-tic);
   printf("Total time (s): %lf\n", stop-start);
+  printf("Bandwidth (GB/s): %lf\n", 1.0E-9*2.0*n*n*nsteps*sizeof(double)/(toc-tic));
   printf(LINE);
 
   // Free the memory
@@ -172,13 +173,14 @@ int main(int argc, char *argv[]) {
 void initial_value(const int n, const double dx, const double length, double * restrict u) {
 
   double y = dx;
+#pragma omp parallel for firstprivate(y)
   for (int j = 0; j < n; ++j) {
+    double y = (double)(j+1)*dx; // Physical y position
     double x = dx; // Physical x position
     for (int i = 0; i < n; ++i) {
       u[i+j*n] = sin(PI * x / length) * sin(PI * y / length);
       x += dx;
     }
-    y += dx; // Physical y position
   }
 
 }
@@ -187,6 +189,7 @@ void initial_value(const int n, const double dx, const double length, double * r
 // Zero the array u
 void zero(const int n, double * restrict u) {
 
+#pragma omp parallel for
   for (int j = 0; j < n; ++j) {
     for (int i = 0; i < n; ++i) {
       u[i+j*n] = 0.0;
@@ -204,6 +207,7 @@ void solve(const int n, const double alpha, const double dx, const double dt, co
   const double r2 = 1.0 - 4.0*r;
 
   // Loop over the nxn grid
+#pragma omp parallel for
   for (int j = 0; j < n; ++j) {
     for (int i = 0; i < n; ++i) {