for-test.c

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#if HAVE_CILK
#include <cilk/cilk.h>
#include <cilk/cilk_api.h>
#endif

typedef struct {
	int max_iter, w, h;
	double xmin, xmax, ymin, ymax;
} global_t;

typedef struct {
	int i, j, k;
} local_t;

int n_calls = 0;

static int compute(void *_g, void *_l)
{
	global_t *g = (global_t*)_g;
	local_t *l = (local_t*)_l;
	double x, x0 = g->xmin + (g->xmax - g->xmin) * l->i / g->w;
	double y, y0 = g->ymin + (g->ymax - g->ymin) * l->j / g->h;
	int k;

	x = x0, y = y0;
	assert(l->k < 0);
	for (k = 0; k < g->max_iter; ++k) {
		double z = x * y;
		x *= x; y *= y;
		if (x + y >= 4) break;
		x = x - y + x0;
		y = z + z + y0; 
	}
	l->k = k;
	__sync_fetch_and_add(&n_calls, 1);
	return 0;
}

void kt_for(int n_threads, int (*func)(void*,void*), void *shared, int n_items, int item_size, void *items);

int main(int argc, char *argv[])
{
	int i, j, tot, type = 0, n_threads = 2;
	global_t global = { 10240*100, 800, 600, -2., -1.2, -1.2, 1.2 };
	local_t *local;

	if (argc > 1) {
		type = argv[1][0] == 'o'? 1 : argv[1][0] == 'c'? 2 : 0;
		if (argv[1][0] >= '0' && argv[1][0] <= '9')
			n_threads = atoi(argv[1]);
	} else {
		fprintf(stderr, "Usage: ./a.out [openmp | cilk | #threads]\n");
	}
	tot = global.w * global.h;
	local = (local_t*)calloc(tot, sizeof(local_t));
	for (j = 0; j < global.h; ++j) {
		local_t *lj = &local[j * global.w];
		for (i = 0; i < global.w; ++i)
			lj[i].i = i, lj[i].j = j, lj[i].k = -1;
	}
	if (type == 0) {
		kt_for(n_threads, compute, &global, tot, sizeof(local_t), local);
	} else if (type == 1) {
		#pragma omp parallel for
		for (i = 0; i < tot; ++i)
			compute(&global, &local[i]);
	} else if (type == 2) {
		#if HAVE_CILK
		cilk_for (i = 0; i < tot; ++i)
			compute(&global, &local[i]);
		#endif
	}
	free(local);
	assert(n_calls == tot);
	return 0;
}