Commits

Ben Bass  committed 69d11d5

add pthread support for using multi-core speedup

  • Participants
  • Parent commits 3151ce6

Comments (0)

Files changed (3)

-all: out.ppm
+all: out.tiff
 
+out.tiff: out.ppm
+	ppm2tiff out.ppm out.tiff
 
 out.ppm: make_ppm mandel.dat
 	./make_ppm
   
 mandel.dat: mandelmap
-	./mandelmap
+	time ./mandelmap
 
 make_ppm: make_ppm.cc mandelmap.h
 	g++ make_ppm.cc -o make_ppm
   
 mandelmap: mandelmap.cc mandelmap.cc mandelmap.h
-	g++ mandelmap.cc -fast -o mandelmap -lc
+	g++ mandelmap.cc -Werror -O3 -o mandelmap -lc -lpthread
 	
 clean:
-	rm -f mandelmap make_ppm out.ppm
+	rm -f mandelmap make_ppm out.ppm out.tiff
 
 superclean: clean
 	rm -f mandel_*.dat

File mandelmap.cc

-#include <stdio.h>
-#include <math.h> // for NAN
+// MandelMap - copyright Ben Bass 2010-2011
+
+#include <cstdio>
+#include <math.h>  // for NAN. Not in the C++ library...
+#include <cstddef>
 #include "mandelmap.h"
 
 #include <sys/types.h>
 #include <unistd.h>
 #include <sys/mman.h>
 
-int mpoint(pinfo* p);
+#include <pthread.h>
 
-int inside_points = 0;
+long mpoint(pinfo* p);
+
+long inside_points = 0;
 
 // these are the dimensions of the fractal we are plotting.
 const FLOAT MIN_X = -2.0;
 const FLOAT MAX_Y = 1.5;
 const FLOAT BAILOUT = 4.0;
 
-inline void mpoint(FLOAT r,
-                   FLOAT i,
-                   pinfo* p,
-                   int old_max_iter=0)
+inline long mpoint(FLOAT r,
+                  FLOAT i,
+                  pinfo* p,
+                  int old_max_iter=0)
 {
     if (isnan(p->x)) {
-        return;
+        // we've already processed this entry -
+        // and it's not part of the Mset
+        return 0;
     }
     FLOAT x2, y2;
     FLOAT x = p->x;
         // carry on later.
         p->x = x;
         p->y = y;
-        inside_points ++;
+        return 1;
     }
     else
     {
         // point (so we can do smooth renderings)
         p->x = NAN;
         p->itercount = (old_max_iter + (MAX_ITER - k));
+        return 0;
     }
 }
 
+int INIT_MAX_ITER;
+pinfo* FPTR_START;
+FLOAT STEP_SIZE;
+pthread_mutex_t acc_lock;
+
+void* worker_start(void* arg)
+{
+    int t_num = *reinterpret_cast<int*>(arg);
+    FLOAT t_start = t_num*STEP_SIZE;
+
+    long local_inside = 0;
+
+    // get the number of points to skip per row
+    ptrdiff_t row_skip = (MAX_X-MIN_X)/STEP_SIZE;
+    // get starting point in the mmapped file - t_num rows in
+    pinfo* local_fptr = FPTR_START+(row_skip*t_num);
+
+    for (FLOAT y=MIN_Y + t_start; y<MAX_Y; y+=(NUM_THREADS*STEP_SIZE))
+    {
+        for (FLOAT x=MIN_X; x<MAX_X; x+=STEP_SIZE)
+        {
+            pinfo* p=local_fptr++;
+            local_inside += mpoint(y, x, p, MAX_ITER);
+        }
+        // skip ahead the relevant number of lines
+        local_fptr += (row_skip * (NUM_THREADS-1));
+    }
+    pthread_mutex_lock(&acc_lock);
+    inside_points += local_inside;
+    pthread_mutex_unlock(&acc_lock);
+}
+
 int main()
 {
     // for multiples of 2^-n (where n is < the number of mantissa digits)
     // floating point can exactly represent the values. Work out the
     // such a value by repeated dividing by two.
+    //
+    // TODO: this is all constant, so use a bit of TMP on it.
     FLOAT step = 1.0;
     unsigned long long invstep = 1;
     for (int i=0; i<BINARY_DIGITS; i++)
         invstep *= 2;
     }
     const FLOAT plot_area = (MAX_Y-MIN_Y)*(MAX_X-MIN_X);  // hopefully an integer.
+    // following assumes we're drawing a square.
     const unsigned long long total_points = plot_area * invstep * invstep;
 
-
     // we will write the data to a file mandel_{BD}.dat where BD
     // is the 'n' in the step size 2^-n
     // a 'single precision' full map (23 bits mantissa) would take
-    // up nearly nearly 10 petabytes (2^50) of storage. Much better
-    // use than most of the data out there, me thinks.
+    // up nearly 10 petabytes (2^50) of storage. Much better use
+    // than most of the data out there, me thinks.
     const unsigned long long mapsize = total_points*sizeof(pinfo) + HEADER_LEN;
     char fname[21];
     snprintf(fname, 20, "mandel_%d.dat", BINARY_DIGITS);
 
     printf("mapped file %d at %p %llu bytes\n", fd, fptr, mapsize);
 
-    unsigned long long z = 0;
-    for (FLOAT y=MIN_Y; y<MAX_Y; y+=step)
+    pthread_t t_handle[NUM_THREADS];
+    int t_num[NUM_THREADS];
+
+    // setup global read-only vars //BADBADBAD (but could be worse)
+    INIT_MAX_ITER = fheader->max_iter;
+    FPTR_START = fptr;
+    STEP_SIZE = step;
+
+    printf("Using %d threads\n", NUM_THREADS);
+
+    pthread_mutex_init(&acc_lock, NULL);
+    for (int nt=0; nt < NUM_THREADS; nt++)
     {
-        for (FLOAT x=MIN_X; x<MAX_X; x+=step)
-        {
-            pinfo* p=fptr++;
-            mpoint(y, x, p, fheader->max_iter);
-        }
+        t_num[nt] = nt;
+        pthread_create(&t_handle[nt],
+                NULL, worker_start, reinterpret_cast<void*>(&t_num[nt]));
+    }
+    for (int nt=0; nt < NUM_THREADS; nt++)
+    {
+        pthread_join(t_handle[nt], NULL);
     }
 
     fheader->max_iter += MAX_ITER;
     // these figures.
     // Should be approx 1.50648
     // https://www.fractalus.com/kerry/articles/area/mandelbrot-area.html
-    printf("%llu %d (%.6f)\n", total_points, inside_points, ((plot_area*inside_points)/float(total_points)));
+    printf("%llu %lu (%.6f)\n", total_points, inside_points, ((plot_area*inside_points)/float(total_points)));
 }
-// MandelMap - copyright Ben Bass 2010
+// MandelMap - copyright Ben Bass 2010-2011
 
 #ifndef MANDEL_MAP_H
 #define MANDEL_MAP_H
 
 typedef double FLOAT;
 
-const int BINARY_DIGITS = 10;
-
-const int MAX_ITER=256;
+const int BINARY_DIGITS = 11;
+const int NUM_THREADS = 2;
+const int MAX_ITER=1024;
 
 struct pinfo
 {