Added files

julia/mandelbrot/cpu_mandelbrot.c

@@ -17,11 +17,10 @@ void cpu_mandelbrot(void *descr[], void *cl_arg)
 
         float centerr = params[0];
         float centeri = params[1];
-
         float offset = params[2];
         float dim = params[3];
         float zoom = width * 0.25296875;
-        float conv_limit = 2.0;
+        float diverge = 4.0;
         int max_iter = (width/2) * 0.049715909 * log10(zoom);
 
         int x,y,n;
@@ -32,27 +31,23 @@ void cpu_mandelbrot(void *descr[], void *cl_arg)
                         float ci = centeri + (y+offset - (dim/2))/zoom;
                         float zr = cr;
                         float zi = ci;
-                        float m = zr * zr + zi * zi;
                         
-                        for (n = 0; n <= max_iter && m < conv_limit * conv_limit; n++) {
-
+                        for (n = 0; n <= max_iter; n++) {
+				if (zr*zr + zi*zi>diverge) break;
                                 float tmp = zr*zr - zi*zi + cr;
                                 zi = 2*zr*zi + ci;
                                 zr = tmp;
-                                m = zr*zr + zi*zi;
                         }
-
+			int color;
+			if (n<max_iter)
+				color = round(15.*n/max_iter);
+			else
+				color = 0;
+			pixels[x*ldP + y] = color;
 		}
-		int color;
-		if (n==max_iter) fprintf(stderr,".");
-		else fprintf(stderr,"%d",n);
-		if (n<max_iter)
-			color = round(15.*n/max_iter);
-		else
-			color = 0;
-		pixels[x*ldP + y] = color;
 	}
 }
+
 char* CPU = "cpu_mandelbrot";
 char* GPU = "gpu_mandelbrot";
 extern char *starpu_find_function(char *name, char *device) {

julia/mandelbrot/mandelbrot.jl

@@ -1,30 +1,114 @@
-function mandelbrotjl(pixels ::Matrix{Int64}, centerr ::Float64, centeri ::Float64)
-    height,width = size(pixels)
-    zoom = width * 0.25296875
-    val_diverge = 2.0
-    max_iterations = (width/2) * 0.049715909 * log10(zoom);
+import Libdl
+using StarPU
+using LinearAlgebra
 
-
-    for y = 1:height
+@target STARPU_CPU+STARPU_CUDA
+@codelet function mandelbrot(pixels ::Matrix{Int64}, params ::Matrix{Float32} ) :: Float32
+    height :: Int64 = height(pixels)
+    width :: Int64 = width(pixels)
+    zoom :: Float64 = width * 0.25296875
+    iz :: Float64 = 1. / zoom
+    diverge :: Float32 = 4.0
+    max_iterations :: Float32 = ((width/2) * 0.049715909 * log10(zoom));
+    imi :: Float32 = 1. / max_iterations
+    centerr :: Float32 = params[1,1]
+    centeri :: Float32 = params[2,1]
+    offset :: Float32 = params[3,1]
+    dim :: Float32 = params[4,1]
+    cr :: Float64 = 0.
+    zr :: Float64 = 0.
+    ci :: Float64 = 0.
+    zi :: Float64 = 0.
+    n :: Int64 = 0
+    tmp :: Float64 = 0.
+    @parallel for y = 1:height
         for x = 1:width
-            cr = centerr + (x - (width / 2))/zoom
+            cr = centerr + (x-1 - (dim / 2)) * iz
             zr = cr
-            ci = centeri + (y - (height / 2))/zoom
+            ci = centeri + (y-1+offset - (dim / 2)) * iz
             zi = ci
-
-            n = 0
-            while ((n < max_iterations) && (zr*zr + zi*zi < val_diverge*val_diverge))
+            for n = 0:max_iterations
+                if (zr*zr + zi*zi > diverge)
+                    break
+                end
                 tmp = zr*zr - zi*zi + cr
                 zi = 2*zr*zi + ci
                 zr = tmp
-                n = n+1
             end
             
             if (n < max_iterations)
-                pixels[y,x] = round(255 * n / max_iterations)
+                pixels[y,x] = round(15 * n * imi)
             else
                 pixels[y,x] = 0
             end
         end
     end
-end
+    return 0. :: Float32
+end
+
+@debugprint "starpu_init"
+starpu_init()
+
+function mandelbrot_with_starpu(A ::Matrix{Int64}, params ::Matrix{Float32}, nslicesx ::Int64)
+    horiz = StarpuDataFilter(STARPU_MATRIX_FILTER_BLOCK, nslicesx)
+    @starpu_block let
+	hA, hP = starpu_data_register(A,params)
+	starpu_data_partition(hA,horiz)
+        starpu_data_partition(hP,horiz)
+        
+	@starpu_sync_tasks for taskx in (1 : nslicesx)
+                @starpu_async_cl mandelbrot(hA[taskx], hP[taskx]) [STARPU_W, STARPU_R]
+	end
+    end
+end
+
+function pixels2img(pixels ::Matrix{Int64}, width ::Int64, height ::Int64, filename ::String)
+    MAPPING = [[66,30,15],[25,7,26],[9,1,47],[4,4,73],[0,7,100],[12,44,138],[24,82,177],[57,125,209],[134,181,229],[211,236,248],[241,233,191],[248,201,95],[255,170,0],[204,128,0],[153,87,0],[106,52,3]]
+    open(filename, "w") do f
+        write(f, "P3\n$width $height\n255\n")
+        for i = 1:height
+            for j = 1:width
+                write(f,"$(MAPPING[1+pixels[i,j]][1]) $(MAPPING[1+pixels[i,j]][2]) $(MAPPING[1+pixels[i,j]][3]) ")
+            end
+            write(f, "\n")
+        end
+    end
+end
+
+function min_times(cr ::Float64, ci ::Float64, dim ::Int64, nslices ::Int64)
+    tmin=0;
+    
+    pixels ::Matrix{Int64} = zeros(dim, dim)
+    params :: Matrix{Float32} = zeros(4*nslices,1)
+    for i=0:(nslices-1)
+        params[4*i+1,1] = cr
+        params[4*i+2,1] = ci
+        params[4*i+3,1] = i*dim/nslices
+        params[4*i+4,1] = dim
+    end
+    for i = 1:10
+        t = time_ns();
+        mandelbrot_with_starpu(pixels, params, nslices)
+        t = time_ns()-t
+        if (tmin==0 || tmin>t)
+            tmin=t
+        end
+    end
+    pixels2img(pixels,dim,dim,"out$(dim).ppm")
+    return tmin
+end
+
+function display_time(cr ::Float64, ci ::Float64, start_dim ::Int64, step_dim ::Int64, stop_dim ::Int64, nslices ::Int64)
+    for dim in (start_dim : step_dim : stop_dim)
+        res = min_times(cr, ci, dim, nslices)
+        res=res/dim/dim; # time per pixel
+        println("$(dim) $(res)")
+    end
+end
+
+
+display_time(-0.800671,-0.158392,32,32,4096,4)
+
+@debugprint "starpu_shutdown"
+starpu_shutdown()
+