fixed redblack files

redblack.cpp

@@ -0,0 +1,113 @@
+#include <iostream>
+#include <cmath>
+#include "redblack.h"
+using namespace std;
+#define TOL 1.0e-6
+
+/* Red-Black 2D algorithm without parallelization
+
+    In: T, b, n, m, itermax (optional)
+        n,m are # rows / colummns
+*/
+void redBlackSerial(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
+        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];
+            }
+        }
+        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
+        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];
+            }
+        }
+        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; }
+
+double computeNorm(double * x, int N) {
+    double norm = 0.0;
+    for (int i=0; i<N; i++) { norm = norm + (x[i]*x[i]); }
+    return sqrt(norm);
+}
+
+double computeResidual(double * T, double * b, int n, int m, double alpha) {
+    int i, j;
+    int index, index_left, index_right, index_up, index_down;
+
+    double residual = 0.0;
+    for (i=1; i<n-1; i++) {
+        for (j=1; j<m-1; j++) {
+            double r;
+            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);
+
+            r = (-1.0*alpha*T[index]) +
+                ((alpha/4.0) *
+                (T[index_left] + T[index_right] +
+                 T[index_up]   + T[index_down] )) + b[index];
+            residual += r*r;
+        }
+    }
+    return sqrt(residual);
+
+}

redblack.h

@@ -0,0 +1,4 @@
+void redBlackSerial(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

@@ -0,0 +1,45 @@
+#include "redblack.h"
+#include <iostream>
+#include <cmath>
+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];
+    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 { T[getIndex(i, j, m)] = 0.0; }
+            b[getIndex(i, j, m)] = 0.0;
+        }
+    }
+    b[getIndex(n/2, m/2, m)] = 1.0;
+
+    double cc = 99.0;
+    int iter = 0;
+    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 ++;
+
+    }
+
+    FILE * output = fopen("T.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)]);
+        }
+    }
+
+
+
+    return 0;
+}