From be2bc936957d4fcbc2001ca48909bb2ce8d3f5c7 Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Fri, 19 Apr 2019 23:32:06 +0200
Subject: Working Futhark! 🎉
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

That was actually easier than expected, and (surprise!) it runs ORDERS
OF MAGNITUDE faster than my C++ implementation...
C++ takes ~1min on my i5-4278U CPU @ 2.60GHz using all 4 virtual cores;
the Futhark takes ~11sec on 1 core. Not sure whether my own code is so
terribly bad or whether Futhark is cool, but this is cool nevertheless.
---
 aberth/aberth.cpp | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 58 insertions(+), 4 deletions(-)

(limited to 'aberth/aberth.cpp')

diff --git a/aberth/aberth.cpp b/aberth/aberth.cpp
index 2adf207..d002ba4 100644
--- a/aberth/aberth.cpp
+++ b/aberth/aberth.cpp
@@ -16,6 +16,10 @@
 #include <cassert>
 #include "../lodepng.h"
 
+extern "C" {
+#include "aberth_kernel.h"
+}
+
 using namespace std;
 
 
@@ -132,7 +136,7 @@ struct AberthState {
 		}
 	}
 
-	// Gauss-Seidel-style step where the updated values are already used in the current iteration
+    // Lagrange-style step where the new elements are computed in parallel from the previous values
 	bool step() {
 		array<Com, N * N> pairs;
 		for (int i = 0; i < N - 1; i++) {
@@ -143,6 +147,7 @@ struct AberthState {
 
 		bool allConverged = true;
 
+		AApprox newapprox;
 		AApprox offsets;
 		for (int i = 0; i < N; i++) {
 			Com pval = eval(poly, approx[i]);
@@ -153,12 +158,15 @@ struct AberthState {
 			for (int j = i + 1; j < N; j++) sum += pairs[N * i + j];
 			offsets[i] = quo / (1.0 - quo * sum);
 
-			approx[i] -= offsets[i];
+			// approx[i] -= offsets[i];
+			newapprox[i] = approx[i] - offsets[i];
 
-			double sval = eval(boundPoly, abs(approx[i]));
+			double sval = eval(boundPoly, abs(newapprox[i]));
 			if (abs(pval) > 1e-9 * sval) allConverged = false;
 		}
 
+		approx = newapprox;
+
 		return allConverged;
 	}
 
@@ -332,6 +340,48 @@ static vector<uint8_t> drawImage(int W, int H, const vector<int> &counts, int ma
 	return image;
 }
 
+static vector<int32_t> invoke_kernel(
+		int32_t width, int32_t height,
+		double left, double top, double right, double bottom,
+		int32_t seed) {
+
+	futhark_context_config *config = futhark_context_config_new();
+	// futhark_context_config_select_device_interactively(config);
+
+	futhark_context *ctx = futhark_context_new(config);
+
+	auto check_ret = [ctx](int ret) {
+		if (ret != 0) {
+			char *str = futhark_context_get_error(ctx);
+			cerr << str << endl;
+			free(str);
+			exit(1);
+		}
+	};
+
+	futhark_i32_1d *dest_arr;
+
+	check_ret(futhark_entry_main_all(
+			ctx, &dest_arr,
+			width, height,
+			left, top, right, bottom,
+			seed));
+
+	check_ret(futhark_context_sync(ctx));
+
+	// Shouldn't free _this_ pointer, apparently
+	int64_t *shape = futhark_shape_i32_1d(ctx, dest_arr);
+	assert(shape[0] == width * height);
+
+	vector<int32_t> buffer(width * height);
+	check_ret(futhark_values_i32_1d(ctx, dest_arr, buffer.data()));
+	check_ret(futhark_free_i32_1d(ctx, dest_arr));
+
+	futhark_context_config_free(config);
+
+	return buffer;
+}
+
 int main(int argc, char **argv) {
 	srandomdev();
 
@@ -339,7 +389,11 @@ int main(int argc, char **argv) {
 	vector<int> counts;
 
 	if (argc <= 1) {
-		tie(W, H, counts) = computeCounts();
+		// tie(W, H, counts) = computeCounts();
+		W = H = 900;
+		Com bottomLeft = Com(-1.5, -1.5);
+		Com topRight = Com(1.5, 1.5);
+		counts = invoke_kernel(W, H, bottomLeft.real(), topRight.imag(), topRight.real(), bottomLeft.imag(), 42);
 		writeCounts(W, H, counts, "out.txt");
 	} else if (argc == 2) {
 		tie(W, H, counts) = readCounts(argv[1]);
-- 
cgit v1.2.3-70-g09d2