add openacc compatible version

Makefile

@@ -2,27 +2,30 @@
 
 OPTIMIZATION=-O3 -flto
 CC=gcc
-CFLAGS=-Isrc/ -lm -lpthread -Wall -Wextra -pedantic-errors $(OPTIMIZATION)
+CFLAGS=-Isrc/ -lm -lpthread -Wall -Wextra -pedantic-errors $(OPTIMIZATION) -DPOINT_DTYPE=float
 
 .PHONY: directories
 
 directories:
 	mkdir -p obj out
 
-obj/point.o: src/point.c src/point.h
+clean:
-	$(CC) $(CFLAGS) -c -o $@ src/point.c
+	rm -rf obj/* out/*
 
-obj/scene.o: src/scene.c src/scene.h
+obj/scene.o: src/scene.c src/scene.h src/point.h
 	$(CC) $(CFLAGS) -c -o $@ src/scene.c
 
-obj/camera.o: src/camera.c src/camera.h
+obj/camera.o: src/camera.c src/camera.h src/point.h
 	$(CC) $(CFLAGS) -c -o $@ src/camera.c
 
-obj/images.o: images/src/images.c images/src/images.h
+obj/images.o: images/src/images.c images/src/images.h src/point.h
 	$(CC) $(CFLAGS) -c -o $@ images/src/images.c
 
-march: obj/camera.o obj/scene.o obj/point.o obj/images.o
+march: obj/camera.o obj/scene.o obj/images.o src/point.h
 	$(CC) $(CFLAGS) -o out/march $^ marcher.c
 
-bench: obj/camera.o obj/scene.o obj/point.o obj/images.o
+bench: obj/camera.o obj/scene.o obj/images.o src/point.h
 	$(CC) $(CFLAGS) -o out/bench $^ bench.c
+
+gpu: obj/camera.o obj/images.o src/point.h
+	$(CC) -fopenacc $(CFLAGS) -o out/gpu $^ gpu.c

gpu.c

@@ -0,0 +1,186 @@
+// use floats instead of doubles
+#define DTYPE float
+// #define POINT_DTYPE DTYPE
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "src/point.h"
+#include "src/camera.h"
+#include "images/src/images.h"
+
+#define ITERS 100
+#define POWER 3
+
+#define SIZE 1000
+#define STEP_WIDTH (2 / ((DTYPE) (SIZE - 1)))
+
+#define CAM_POSITION 1.35
+
+#define STEPS 100
+#define THRESHOLD 0.001
+
+#include <openacc.h>
+
+
+float mandelbulb_dist(struct point pt) 
+{
+    DTYPE power = POWER;
+
+    struct point z = pt;
+	DTYPE dr = 1.0;
+	DTYPE r = 0.0;
+	for (int i = 0; i < ITERS ; i++) {
+		r = pt_length(z);
+
+		if (r>2) {
+            break;
+        }
+		
+		// convert to polar coordinates
+		DTYPE theta = acos(z.z/r);
+		DTYPE phi = atan2(z.y,z.x);
+		dr =  pow(r, power-1.0)*power*dr + 1.0;
+		
+		// scale and rotate the point
+		DTYPE zr = pow(r, power);
+		theta = theta*power;
+		phi = phi*power;
+		
+		// convert back to cartesian coordinates
+        z = (struct point) {
+                .x = sin(theta)*cos(phi) * zr + pt.x,
+                .y = sin(phi)*sin(theta) * zr + pt.y,
+                .z = cos(theta) * zr + pt.z
+        };
+	}
+	
+    return 0.5*log(r)*r/dr;
+}
+
+struct setup {
+    struct point p0;            // origin
+    struct point direction;     // ray direction
+    struct point x;             // ray movement in col
+    struct point y;             // ray movement in row
+};
+
+struct setup make_setup()
+{
+    // set up camera
+    struct camera cam;
+    cam.fov = 90;
+
+    camera_set_looking_at(&cam, (struct point){.x=CAM_POSITION, .y= CAM_POSITION, .z = CAM_POSITION}, PT_ZERO);
+
+    struct point span_z, span_xy;
+
+    // get rotation axis
+    pt_orthogonal_plane(cam.direction, &span_z, &span_xy);
+
+    printf("rendering %ix%ipx\n", SIZE, SIZE);
+
+    // distance each ray has from anothe on the ortogonal plane
+    //DTYPE step_dist = 2 / (DTYPE) (width - 1);
+
+    // vectors to move on the projection plane
+    struct point move_right = pt_scale(span_xy, STEP_WIDTH);
+    struct point move_up =    pt_scale(span_z,  STEP_WIDTH);;
+
+    // set starting point
+    struct point starting_point = pt_normalize(cam.direction);
+
+    // rotate starting point to (0,0)
+    starting_point = pt_add(starting_point, pt_mult(move_right, - SIZE / (DTYPE) 2));
+    starting_point = pt_add(starting_point, pt_mult(move_up, - SIZE / (DTYPE) 2));
+
+    return (struct setup) {
+        cam.location, starting_point, span_xy, span_z
+    };
+}
+
+int main()
+{
+    struct setup setup = make_setup();
+    struct point start = setup.p0;
+    
+    //printf("device num acc_device_current: %d\n", acc_get_num_devices(acc_device_current));
+    //printf("device num acc_device_none: %d\n", acc_get_num_devices(acc_device_none));
+    //printf("device num acc_device_default: %d\n", acc_get_num_devices(acc_device_default));
+    //printf("device num acc_device_host: %d\n", acc_get_num_devices(acc_device_host));
+    //printf("device num acc_device_not_host: %d\n", acc_get_num_devices(acc_device_not_host));
+    //printf("device num acc_device_nvidia: %d\n", acc_get_num_devices(acc_device_nvidia));
+    //printf("device num acc_device_radeon: %d\n", acc_get_num_devices(acc_device_radeon));
+
+    // get backing buff
+    int* buff = calloc(sizeof(int), SIZE * SIZE);
+
+    // indicate malloc failure
+    if (buff == NULL)
+        return -255;
+
+    printf("Before kernel\n");
+
+#pragma acc data copy(buff)
+{
+    // kernel goes brr
+    #pragma acc kernels
+    for (int x = 0; x < SIZE; x++) {
+        for (int y = 0; y < SIZE; y++) {
+            // get direction
+            struct point offset = pt_add(pt_mult(setup.x, STEP_WIDTH * x), pt_mult(setup.y, STEP_WIDTH * y));
+            struct point direction = pt_add(setup.direction, offset);
+            // get start
+            struct point loc = start;
+            // march!
+            DTYPE dist;
+            int res = -1;
+
+            for (int i = 0; i < STEPS; i++) {
+                dist = mandelbulb_dist(loc);
+                if (dist < THRESHOLD) {
+                    res = i;
+                    break;
+                }
+                if (dist > 100) {
+                    break;
+                }
+                loc = pt_add(loc, pt_scale(direction, dist));
+            }
+            buff[y * SIZE + x] = res;
+        }
+    }
+}
+
+    printf("after kernel\n");
+
+    // convert distance field into image
+
+    Image img;
+
+    // initialize shared pixel buffer
+    image_new(SIZE, SIZE, &img);
+    
+    #pragma acc parallel
+    for (unsigned int x = 0; x < SIZE; x++) {
+        #pragma acc loop
+        for (unsigned int y = 0; y < SIZE; y++) {
+            if (buff[y * SIZE + x] < 0) {
+                image_set_px(img, x, y, 0, 0, 0);
+            } else {
+                // float in range of [0,1]
+                DTYPE fac = buff[y * SIZE + x] / (DTYPE) STEPS;
+                // calc shade 
+                int shade = ((1-fac) * 255);
+                image_set_px(img, x, y, shade, shade, shade);
+            }
+        }
+    }
+
+    printf("after encoding\n");
+
+
+    image_save_bmp(img, "gpu-goes-brrrrrrr.bmp");
+
+    return 0;
+}

src/point.c

@@ -1,118 +0,0 @@
-#include <math.h>
-#include <stdio.h>
-#include "point.h"
-
-// get the length of vector
-inline double pt_length_inline (struct point pt) __attribute__((always_inline));
-
-inline struct point pt_mult(struct point pt, double scalar) {
-    return (struct point) {
-        .x = pt.x * scalar,
-        .y = pt.y * scalar,
-        .z = pt.z * scalar
-    };
-}
-
-// return internal angle between a and b
-inline double pt_angle(struct point a, struct point b) {
-    return acos(pt_dot(
-        pt_normalize(a),
-        pt_normalize(b)
-    ));
-}
-
-// get the length of vector
-inline double pt_length_inline (struct point pt) {
-    return sqrt((pt.x * pt.x) + (pt.y * pt.y) + (pt.z * pt.z));
-}
-
-double pt_length(struct point pt) {
-    return pt_length_inline(pt);
-}
-
-// add the vector add to the vector pt
-inline struct point pt_add(struct point pt, struct point add) {
-    return (struct point) {
-        .x = pt.x + add.x,
-        .y = pt.y + add.y,
-        .z = pt.z + add.z,
-    };
-}
-
-// add the vector add to the vector pt
-inline struct point pt_sub(struct point pt, struct point sub) {
-    return (struct point) {
-        .x = pt.x - sub.x,
-        .y = pt.y - sub.y,
-        .z = pt.z - sub.z,
-    };
-}
-
-inline double pt_dist(struct point p1, struct point p2) {
-    return pt_length_inline(pt_sub(p1, p2));
-}
-
-// normalize a vector
-inline struct point pt_normalize(struct point pt) {
-    double length = pt_length(pt);
-
-    return (struct point) {
-        .x = pt.x / length,
-        .y = pt.y / length,
-        .z = pt.z / length
-    };
-}
-
-// scale vector to length
-inline struct point pt_scale(struct point pt, double length) {
-    double f = length / pt_length_inline(pt);
-    return (struct point) {
-            .x = pt.x * f,
-            .y = pt.y * f,
-            .z = pt.z * f
-    };
-}
-
-// dot product of two vectors
-inline double pt_dot(struct point a, struct point b) {
-    return a.x*b.x + a.y*b.y + a.z*b.z;
-}
-
-// cross product of two vectors
-inline struct point pt_cross(struct point a, struct point b) {
-    return (struct point) {
-        .x = a.y*b.z - a.z*b.y,
-        .y = a.z*b.x - a.x*b.z,
-        .z = a.x*b.y - a.y*b.x
-    };
-}
-
-inline void pt_print(struct point pt) {
-    printf("(%f, %f, %f)\n", pt.x, pt.y, pt.z);
-}
-
-inline void pt_print_n(const char* name, struct point pt) {
-    printf("%s: (%f, %f, %f)\n", name, pt.x, pt.y, pt.z);
-}
-
-// find two vectors that span the orthogonal plane, where 
-// span_xy is a vector lying on the xy-plane (and pointing left)
-// and span_z is orthogonal to span_xy pointing "upwards"
-inline void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy) {
-    pt = pt_normalize(pt);
-
-    // get the vector lying on the xy axis
-    // this is done by
-    *span_xy = pt_normalize(pt_cross(PT_NEW(0,0,1), pt)); // points to the "left" (of the viewing direction)
-
-    // now use this, to find the vector 
-    *span_z = pt_normalize(pt_cross(pt, *span_xy));
-}
-
-inline struct point pt_mod(struct point pt, double mod) {
-    return (struct point) {
-        .x = fabs(fmod(pt.x, mod)),
-        .y = fabs(fmod(pt.y, mod)),
-        .z = fabs(fmod(pt.z, mod))
-    };
-}

src/point.h

@@ -1,26 +1,149 @@
 #pragma once
 
+#ifndef POINT_DTYPE
+#define POINT_DTYPE double
+#endif
+
 
 struct point {
-    double x;
+    POINT_DTYPE x;
-    double y;
+    POINT_DTYPE y;
-    double z;
+    POINT_DTYPE z;
 };
 
 #define PT_ZERO ((struct point) {.x=0, .y=0, .z=0})
 #define PT_NEW(x,y,z) ((struct point){(x), (y), (z)})
 
-struct point pt_scale(struct point pt, double length);
+static struct point pt_scale(struct point pt, POINT_DTYPE length);
-struct point pt_normalize(struct point pt);
+static struct point pt_normalize(struct point pt);
-struct point pt_mult(struct point pt, double scalar);
+static struct point pt_mult(struct point pt, POINT_DTYPE scalar);
-double pt_length(struct point pt);
+static POINT_DTYPE pt_length(struct point pt);
-struct point pt_add(struct point pt, struct point add);
+static struct point pt_add(struct point pt, struct point add);
-struct point pt_sub(struct point pt, struct point sub);
+static struct point pt_sub(struct point pt, struct point sub);
-double pt_dist(struct point p1, struct point p2);
+static POINT_DTYPE pt_dist(struct point p1, struct point p2);
-struct point pt_mod(struct point pt, double mod);
+static struct point pt_mod(struct point pt, POINT_DTYPE mod);
-double pt_dot(struct point a, struct point b);
+static POINT_DTYPE pt_dot(struct point a, struct point b);
-struct point pt_cross(struct point a, struct point b);
+static struct point pt_cross(struct point a, struct point b);
-double pt_angle(struct point a, struct point b);
+static POINT_DTYPE pt_angle(struct point a, struct point b);
-void pt_print(struct point pt);
+static void pt_print(struct point pt);
-void pt_print_n(const char* name, struct point pt);
+static void pt_print_n(const char* name, struct point pt);
-void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy);
+static void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy);
+
+#include <math.h>
+#include <stdio.h>
+#include "point.h"
+
+// get the length of vector
+static inline POINT_DTYPE pt_length_inline (struct point pt) __attribute__((always_inline));
+
+static inline struct point pt_mult(struct point pt, POINT_DTYPE scalar) {
+    return (struct point) {
+        .x = pt.x * scalar,
+        .y = pt.y * scalar,
+        .z = pt.z * scalar
+    };
+}
+
+// return internal angle between a and b
+static inline POINT_DTYPE pt_angle(struct point a, struct point b) {
+    return acos(pt_dot(
+        pt_normalize(a),
+        pt_normalize(b)
+    ));
+}
+
+// get the length of vector
+static inline POINT_DTYPE pt_length_inline (struct point pt) {
+    return sqrt((pt.x * pt.x) + (pt.y * pt.y) + (pt.z * pt.z));
+}
+
+static inline POINT_DTYPE pt_length(struct point pt) {
+    return pt_length_inline(pt);
+}
+
+// add the vector add to the vector pt
+static inline struct point pt_add(struct point pt, struct point add) {
+    return (struct point) {
+        .x = pt.x + add.x,
+        .y = pt.y + add.y,
+        .z = pt.z + add.z,
+    };
+}
+
+// add the vector add to the vector pt
+static inline struct point pt_sub(struct point pt, struct point sub) {
+    return (struct point) {
+        .x = pt.x - sub.x,
+        .y = pt.y - sub.y,
+        .z = pt.z - sub.z,
+    };
+}
+
+static inline POINT_DTYPE pt_dist(struct point p1, struct point p2) {
+    return pt_length_inline(pt_sub(p1, p2));
+}
+
+// normalize a vector
+static inline struct point pt_normalize(struct point pt) {
+    POINT_DTYPE length = pt_length(pt);
+
+    return (struct point) {
+        .x = pt.x / length,
+        .y = pt.y / length,
+        .z = pt.z / length
+    };
+}
+
+// scale vector to length
+static inline struct point pt_scale(struct point pt, POINT_DTYPE length) {
+    POINT_DTYPE f = length / pt_length_inline(pt);
+    return (struct point) {
+            .x = pt.x * f,
+            .y = pt.y * f,
+            .z = pt.z * f
+    };
+}
+
+// dot product of two vectors
+static inline POINT_DTYPE pt_dot(struct point a, struct point b) {
+    return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+// cross product of two vectors
+static inline struct point pt_cross(struct point a, struct point b) {
+    return (struct point) {
+        .x = a.y*b.z - a.z*b.y,
+        .y = a.z*b.x - a.x*b.z,
+        .z = a.x*b.y - a.y*b.x
+    };
+}
+
+static inline void pt_print(struct point pt) {
+    printf("(%f, %f, %f)\n", pt.x, pt.y, pt.z);
+}
+
+static inline void pt_print_n(const char* name, struct point pt) {
+    printf("%s: (%f, %f, %f)\n", name, pt.x, pt.y, pt.z);
+}
+
+// find two vectors that span the orthogonal plane, where 
+// span_xy is a vector lying on the xy-plane (and pointing left)
+// and span_z is orthogonal to span_xy pointing "upwards"
+static inline void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy) {
+    pt = pt_normalize(pt);
+
+    // get the vector lying on the xy axis
+    // this is done by
+    *span_xy = pt_normalize(pt_cross(PT_NEW(0,0,1), pt)); // points to the "left" (of the viewing direction)
+
+    // now use this, to find the vector 
+    *span_z = pt_normalize(pt_cross(pt, *span_xy));
+}
+
+static inline struct point pt_mod(struct point pt, POINT_DTYPE mod) {
+    return (struct point) {
+        .x = fabs(fmod(pt.x, mod)),
+        .y = fabs(fmod(pt.y, mod)),
+        .z = fabs(fmod(pt.z, mod))
+    };
+}