A better way of doing GPU symmetry

4 files changed, 71 insertions, 204 deletions

diff --git a/data/diffraction.cl b/data/diffraction.cl
index 76c6b9d9..994e27db 100644
--- a/data/diffraction.cl
+++ b/data/diffraction.cl
@@ -136,7 +136,6 @@ float molecule_factor(global float *intensities, global float *flags,
 	float hf, kf, lf;
 	int h, k, l;
 	float val = 0.0;
-	signed int i;
 
 	#ifdef FLAT_INTENSITIES
 	return 100.0;
@@ -149,150 +148,9 @@ float molecule_factor(global float *intensities, global float *flags,
 	h = round(hf);
 	k = round(kf);
 	l = round(lf);
-	i = -h-k;
-
-	#ifdef PG1
-	val += lookup_flagged_intensity(intensities, flags, h, k, l);
-	#endif /* PG1 */
-
-	#ifdef PG1BAR
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, -l);
-	#endif /* PG1BAR */
-
-	#ifdef PGMMM
-	val += lookup_flagged_intensity(intensities, flags, h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, k, -l);
-	val += lookup_flagged_intensity(intensities, flags, h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags, h, -k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags, h, k, -l);
-	#endif /* PGMMM */
-
-	#ifdef PG4
-	val += lookup_flagged_intensity(intensities, flags, -k, h, l);
-	val += lookup_flagged_intensity(intensities, flags, h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, l);
-	val += lookup_flagged_intensity(intensities, flags, k, -h, l);
-	#endif /* PG4 */
-
-	#ifdef PG422
-	val += lookup_flagged_intensity(intensities, flags, -k, h, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -h, -l);
-	val += lookup_flagged_intensity(intensities, flags, k, h, -l);
-	val += lookup_flagged_intensity(intensities, flags, k, -h, -l);
-	val += lookup_flagged_intensity(intensities, flags, h, -k, -l);
-	#endif /* PG422 */
-
-	#ifdef PG321H
-	val += lookup_flagged_intensity(intensities, flags, h, k, l);
-	val += lookup_flagged_intensity(intensities, flags, i, h, l);
-	val += lookup_flagged_intensity(intensities, flags, k, h, -l);
-	val += lookup_flagged_intensity(intensities, flags, k, i, l);
-	val += lookup_flagged_intensity(intensities, flags, i, k, -l);
-	val += lookup_flagged_intensity(intensities, flags, h, i, -l);
-	#endif /* PG321H */
-
-	#ifdef PG6
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  h,  l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  i,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -h,  l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -i,  l);
-	#endif /* PG6 */
-
-	#ifdef PG6M
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  h,  l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  i,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -h,  l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -i,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -h, -l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -i, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  h, -l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  i, -l);
-	#endif /* PG6M */
-
-	#ifdef PG6MMM
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  h,  l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  i,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -h,  l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -i,  l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  h, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h,  i, -l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -h, -l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -i, -l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -h, -l);
-	val += lookup_flagged_intensity(intensities, flags, -k,  i, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  h, -l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  i, -l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -h,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -i,  l);
-	val += lookup_flagged_intensity(intensities, flags, -i, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  h,  l);
-	val += lookup_flagged_intensity(intensities, flags,  h,  i,  l);
-	val += lookup_flagged_intensity(intensities, flags,  i,  k,  l);
-	#endif /* PG6MMM */
-
-	#ifdef PG23
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  l,  h);
-	val += lookup_flagged_intensity(intensities, flags,  l,  h,  k);
-	val += lookup_flagged_intensity(intensities, flags, -k, -l,  h);
-	val += lookup_flagged_intensity(intensities, flags, -l, -h,  k);
-	val += lookup_flagged_intensity(intensities, flags, -k,  l, -h);
-	val += lookup_flagged_intensity(intensities, flags, -l,  h, -k);
-	val += lookup_flagged_intensity(intensities, flags,  k, -l, -h);
-	val += lookup_flagged_intensity(intensities, flags,  l, -h, -k);
-	#endif /* PG23 */
-
-	#ifdef PGM3
-	val += lookup_flagged_intensity(intensities, flags,  h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  k,  l,  h);
-	val += lookup_flagged_intensity(intensities, flags,  l,  h,  k);
-	val += lookup_flagged_intensity(intensities, flags, -k, -l,  h);
-	val += lookup_flagged_intensity(intensities, flags, -l, -h,  k);
-	val += lookup_flagged_intensity(intensities, flags, -k,  l, -h);
-	val += lookup_flagged_intensity(intensities, flags, -l,  h, -k);
-	val += lookup_flagged_intensity(intensities, flags,  k, -l, -h);
-	val += lookup_flagged_intensity(intensities, flags,  l, -h, -k);
-
-	val += lookup_flagged_intensity(intensities, flags, -h, -k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h,  k, -l);
-	val += lookup_flagged_intensity(intensities, flags,  h, -k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -h,  k,  l);
-	val += lookup_flagged_intensity(intensities, flags, -k, -l, -h);
-	val += lookup_flagged_intensity(intensities, flags, -l, -h, -k);
-	val += lookup_flagged_intensity(intensities, flags,  k,  l, -h);
-	val += lookup_flagged_intensity(intensities, flags,  l,  h, -k);
-	val += lookup_flagged_intensity(intensities, flags,  k, -l,  h);
-	val += lookup_flagged_intensity(intensities, flags,  l, -h,  k);
-	val += lookup_flagged_intensity(intensities, flags, -k,  l,  h);
-	val += lookup_flagged_intensity(intensities, flags, -l,  h,  k);
-	#endif /* PGM3 */
-
-	/* FIXME: Add the remaining point groups */
+
+	/* Symmetry stuff goes here */
+	INSERT_HERE
 
 	return val;
 	#endif /* FLAT_INTENSITIIES */
diff --git a/src/cl-utils.c b/src/cl-utils.c
index 97efe2f8..9f984169 100644
--- a/src/cl-utils.c
+++ b/src/cl-utils.c
@@ -3,11 +3,11 @@
  *
  * OpenCL utility functions
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -215,7 +215,8 @@ static void show_build_log(cl_program prog, cl_device_id dev)
 
 
 cl_program load_program(const char *filename, cl_context ctx,
-                        cl_device_id dev, cl_int *err, const char *extra_cflags)
+                        cl_device_id dev, cl_int *err, const char *extra_cflags,
+                        const char *insert_stuff)
 {
 	FILE *fh;
 	cl_program prog;
@@ -223,6 +224,8 @@ cl_program load_program(const char *filename, cl_context ctx,
 	size_t len;
 	cl_int r;
 	char cflags[1024] = "";
+	char *insert_pos;
+	size_t il;
 
 	fh = fopen(filename, "r");
 	if ( fh == NULL ) {
@@ -231,10 +234,31 @@ cl_program load_program(const char *filename, cl_context ctx,
 		return 0;
 	}
 	source = malloc(16384);
+	if ( source == NULL ) return 0;
 	len = fread(source, 1, 16383, fh);
 	fclose(fh);
 	source[len] = '\0';
 
+	if ( insert_stuff != NULL ) {
+		insert_pos = strstr(source, "INSERT_HERE");
+
+		if ( insert_pos != NULL ) {
+
+			char *source2;
+			source2 = malloc(strlen(source)+strlen(insert_stuff)+1);
+			if ( source2 == NULL ) return 0;
+
+			il = insert_pos - source;
+			memcpy(source2, source, il);
+			memcpy(source2+il, insert_stuff, strlen(insert_stuff)+1);
+			memcpy(source2+il+strlen(insert_stuff),
+			       source+il+11, strlen(source+il+11)+1);
+			source = source2;
+			free(source);
+
+		}
+	}
+
 	prog = clCreateProgramWithSource(ctx, 1, (const char **)&source,
 	                                 NULL, err);
 	if ( *err != CL_SUCCESS ) {
diff --git a/src/cl-utils.h b/src/cl-utils.h
index e5546b01..6e76748a 100644
--- a/src/cl-utils.h
+++ b/src/cl-utils.h
@@ -3,11 +3,11 @@
  *
  * OpenCL utility functions
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -38,7 +38,7 @@ extern const char *clError(cl_int err);
 extern cl_device_id get_cl_dev(cl_context ctx, int n);
 extern cl_program load_program(const char *filename, cl_context ctx,
                         cl_device_id dev, cl_int *err,
-                        const char *extra_cflags);
+                        const char *extra_cflags, const char *insert_stuff);
 
 
 #endif	/* CLUTILS_H */
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 8b165856..f031e07f 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -381,6 +381,7 @@ struct gpu_context *setup_gpu(int no_sfac,
 	size_t maxwgsize;
 	int i;
 	char cflags[512] = "";
+	char *insert_stuff = NULL;
 
 	STATUS("Setting up GPU...\n");
 
@@ -440,57 +441,41 @@ struct gpu_context *setup_gpu(int no_sfac,
 
 	if ( sym != NULL ) {
 
-		/* Triclinic */
-		if ( strcmp(sym, "1") == 0 ) {
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "-1") == 0 ) {
-			strncat(cflags, "-DPG1BAR ", 511-strlen(cflags));
-
-		/* Monoclinic */
-		/* FIXME: 2, 2/m, (m)*/
-
-		/* Orthorhombic */
-		} else if ( strcmp(sym, "mmm") == 0 ) {
-			strncat(cflags, "-DPGMMM ", 511-strlen(cflags));
-		/* FIXME: 222, mmm, (mm2) */
-
-		/* Tetragonal */
-		} else if ( strcmp(sym, "4") == 0 ) {
-			strncat(cflags, "-DPG4 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "422") == 0 ) {
-			strncat(cflags, "-DPG422 ", 511-strlen(cflags));
-		/* FIXME: 4/m, 4/mmm, (-42m, -4m2, -4, 4mm) */
-
-		/* Trigonal (rhombohedral) */
-		/* FIXME: 3, 32, -3m, (3m) */
-
-		/* Trigonal (hexagonal) */
-		} else if ( strcmp(sym, "321_H") == 0 ) {
-			strncat(cflags, "-DPG321H ", 511-strlen(cflags));
-		/* FIXME: 3, -3, 312, -31m -3m1, (3m1, 31m) */
-
-		/* Hexagonal */
-		} else if ( strcmp(sym, "6") == 0 ) {
-			strncat(cflags, "-DPG6 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "6/m") == 0 ) {
-			strncat(cflags, "-DPG6M ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "6/mmm") == 0 ) {
-			strncat(cflags, "-DPG6MMM ", 511-strlen(cflags));
-		/* FIXME: 622, (-6, -6m2, -62m) */
-
-		/* Cubic */
-		} else if ( strcmp(sym, "23") == 0 ) {
-			strncat(cflags, "-DPG23 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "m-3") == 0 ) {
-			strncat(cflags, "-DPGM3 ", 511-strlen(cflags));
-		/* FIXME: 432, m-3m, (-43m) */
-
-		} else {
-			ERROR("Sorry!  Point group '%s' is not currently"
-			      " supported on the GPU."
-			      " I'm using '1' instead.\n", sym);
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
+		int i, n;
+		SymOpList *pg;
+		size_t islen = 0;
+
+		insert_stuff = malloc(16384);
+		if ( insert_stuff == NULL ) return NULL;
+		insert_stuff[0] = '\0';
+
+		pg = get_pointgroup(sym);
+		n = num_equivs(pg, NULL);
+		for ( i=0; i<n; i++ ) {
+
+			IntegerMatrix *op = get_symop(pg, NULL, i);
+			char line[1024];
+
+			snprintf(line, 1023,
+			         "val += lookup_flagged_intensity(intensities, "
+			         "flags, %s, %s, %s);\n\t",
+			         get_matrix_name(op, 0),
+				 get_matrix_name(op, 1),
+				 get_matrix_name(op, 2));
+
+			islen += strlen(line);
+			if ( islen > 16383 ) {
+				ERROR("Too many symmetry operators.\n");
+				return NULL;
+			}
+			strcat(insert_stuff, line);
+
 		}
+
+		free_symoplist(pg);
+
+		printf("Inserting --->%s<---\n", insert_stuff);
+
 	} else {
 		if ( intensities != NULL ) {
 			ERROR("You gave me an intensities file but no point"
@@ -522,7 +507,7 @@ struct gpu_context *setup_gpu(int no_sfac,
 	free(flags_ptr);
 
 	gctx->prog = load_program(DATADIR"/crystfel/diffraction.cl", gctx->ctx,
-	                          dev, &err, cflags);
+	                          dev, &err, cflags, insert_stuff);
 	if ( err != CL_SUCCESS ) {
 		free(gctx);
 		return NULL;