diff options
author | tomsmeding <tom.smeding@gmail.com> | 2019-04-19 23:32:06 +0200 |
---|---|---|
committer | tomsmeding <tom.smeding@gmail.com> | 2019-04-19 23:34:01 +0200 |
commit | be2bc936957d4fcbc2001ca48909bb2ce8d3f5c7 (patch) | |
tree | 9b2402ff15cc074e9c7be6de14affc57c735f1bf /aberth/aberth.cpp | |
parent | bc72c77a8a3d6a82cd234d40308dd8c47aa0e742 (diff) |
Working Futhark! 🎉
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.
Diffstat (limited to 'aberth/aberth.cpp')
-rw-r--r-- | aberth/aberth.cpp | 62 |
1 files changed, 58 insertions, 4 deletions
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]); |