parallel red-black attempt #1

redblack.cpp

@@ -1,6 +1,7 @@
 #include <iostream>
 #include <cmath>
 #include "redblack.h"
+#include <omp.h>
 using namespace std;
 #define TOL 1.0e-6
 
@@ -78,6 +79,88 @@ void redBlackSerial(double * T, double * b,
         iter ++;
     }
 }
+/* Red-Black 2D algorithm with parallelization
+
+    In: T, b, n, m, itermax (optional)
+        n,m are # rows / colummns
+*/
+void redBlackParallel(double * T, double * b,
+                    int n, int m, double alpha, int itermax) {
+    int i, j;
+    int index, index_left, index_right, index_up, index_down;
+
+    int iter = 0;
+    while (iter < itermax) {
+        // update half the board in parallel
+        #pragma omp parallel for collapse(2) \
+                private(i, j) shared(T, b, n, m, alpha)
+        for (i=1; i<n-1; i+=2) {
+            for (j=1; j<m-1; j+=2) {
+                index       = getIndex(i, j, m);
+                index_left  = getIndex(i-1, j, m);
+                index_right = getIndex(i+1, j, m);
+                index_up    = getIndex(i, j+1, m);
+                index_down  = getIndex(i, j-1, m);
+
+                T[index] = (1.0 - alpha) * T[index] +
+                           ((alpha/4.0) *
+                           (T[index_left] + T[index_right] +
+                            T[index_up]   + T[index_down] )) + b[index];
+            }
+        }
+        #pragma omp parallel for collapse(2) \
+                private(i, j) shared(T, b, n, m, alpha)
+        for (i=2; i<n-1; i+=2) {
+            for (j=2; j<m-1; j+=2) {
+                index       = getIndex(i, j, m);
+                index_left  = getIndex(i-1, j, m);
+                index_right = getIndex(i+1, j, m);
+                index_up    = getIndex(i, j+1, m);
+                index_down  = getIndex(i, j-1, m);
+
+                T[index] = (1.0 - alpha) * T[index] +
+                           ((alpha/4.0) *
+                           (T[index_left] + T[index_right] +
+                            T[index_up]   + T[index_down] )) + b[index];
+            }
+        }
+
+        // update second half of the board
+        #pragma omp parallel for collapse(2) \
+                private(i, j) shared(T, b, n, m, alpha)
+        for (i=1; i<n-1; i+=2) {
+            for(j=2; j<m-1; j+=2) {
+                index       = getIndex(i, j, m);
+                index_left  = getIndex(i-1, j, m);
+                index_right = getIndex(i+1, j, m);
+                index_up    = getIndex(i, j+1, m);
+                index_down  = getIndex(i, j-1, m);
+
+                T[index] = (1.0 - alpha) * T[index] +
+                           ((alpha/4.0) *
+                           (T[index_left] + T[index_right] +
+                            T[index_up]   + T[index_down] )) + b[index];
+            }
+        }
+        #pragma omp parallel for collapse(2) \
+                private(i, j) shared(T, b, n, m, alpha)
+        for (i=2; i<n-1; i+=2) {
+            for (j=1; j<m-1; j+=2) {
+                index       = getIndex(i, j, m);
+                index_left  = getIndex(i-1, j, m);
+                index_right = getIndex(i+1, j, m);
+                index_up    = getIndex(i, j+1, m);
+                index_down  = getIndex(i, j-1, m);
+
+                T[index] = (1.0 - alpha) * T[index] +
+                           ((alpha/4.0) *
+                           (T[index_left] + T[index_right] +
+                            T[index_up]   + T[index_down] )) + b[index];
+            }
+        }
+        iter ++;
+    }
+}
 
 int getIndex(int x, int y, int width) { return x*width + y; }
 

redblack.h

@@ -1,4 +1,6 @@
 void redBlackSerial(double * T, double * b,
                    int n, int m, double alpha, int itermax);
+void redBlackParallel(double * T, double * b,
+                       int n, int m, double alpha, int itermax);
 int getIndex(int x, int y, int width);
 double computeResidual(double * T, double * b, int n, int m, double alpha);

redblack_example.cpp

@@ -1,43 +1,78 @@
 #include "redblack.h"
 #include <iostream>
 #include <cmath>
+#include <ctime>
 using namespace std;
 
-
-int main() {
-    const int n = 100;
-    const int m = 100;
-    double alpha = 0.9;
-    double * T = new double[n*m];
-    double * b = new double[n*m];
+void initArrays(double * T, double * b, int n, int m) {
     for (int i=0; i<n; i++) {
         for (int j=0; j<m; j++) {
             if (i==0 || i==n-1) { T[getIndex(i, j, m)] = 0.0; }
-			else if (j==0 || j==m-1) { T[getIndex(i, j, m)] = 20.0; }
+            else if (j==0 || j==m-1) { T[getIndex(i, j, m)] = 20.0; }
             else { T[getIndex(i, j, m)] = 0.0; }
             b[getIndex(i, j, m)] = 0.0;
         }
     }
     b[getIndex(n/2, m/2, m)] = 1.0;
+}
+
+int main() {
+    const int n = 100;
+    const int m = 100;
+    double alpha = 0.9;
+    double * T = new double[n*m];
+    double * b = new double[n*m];
+    initArrays(T, b, n, m);
 
     double cc = 99.0;
     int iter = 0;
+    cout << "Starting without parallel!\n";
+    clock_t start;
+    double duration;
+    start = clock();
     while (cc > 1e-5) {
         redBlackSerial(T, b, n, m, alpha, 1);
         cc = computeResidual(T, b, n, m, alpha);
         if (iter%1000==0) {
             cout << "Iteration: " << iter << ", Res: " << cc << "\n";
         }
         iter ++;
+    }
+    duration = (clock() - start) / (double) CLOCKS_PER_SEC;
+    cout << "Non-parallel: " << duration << " seconds.\n";
+
+    FILE * output = fopen("T_serial.txt", "w");
+    for (int i=0; i<n; i++) {
+        for (int j=0; j<m; j++) {
+            fprintf(output, "%i %i %0.16g\n", i, j, T[getIndex(i, j, m)]);
+        }
+    }
+    fclose(output);
 
+    initArrays(T, b, n, m);
+
+    cc = 99.0;
+    iter = 0;
+    cout << "Starting with parallel!\n";
+    start = clock();
+    while (cc > 1e-5) {
+        redBlackParallel(T, b, n, m, alpha, 1);
+        cc = computeResidual(T, b, n, m, alpha);
+        if (iter%1000==0) {
+            cout << "Iteration: " << iter << ", Res: " << cc << "\n";
+        }
+        iter ++;
     }
+    duration = (clock() - start) / (double) CLOCKS_PER_SEC;
+    cout << "Parallel: " << duration << " seconds.\n";
 
-    FILE * output = fopen("T.txt", "w");
+    output = fopen("T_parallel.txt", "w");
     for (int i=0; i<n; i++) {
         for (int j=0; j<m; j++) {
             fprintf(output, "%i %i %0.16g\n", i, j, T[getIndex(i, j, m)]);
         }
     }
+    fclose(output);