From 2c238039c2efda1788ea72c9fb41ff354acc8e97 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 16:04:46 +0100 Subject: Move the "indexed reflection array" thing to where it can't do any harm --- Makefile.am | 8 +- Makefile.in | 28 +- data/defs.h | 23 -- data/diffraction.cl | 6 +- libcrystfel/Makefile.am | 18 +- libcrystfel/Makefile.in | 31 +-- libcrystfel/src/diffraction-gpu.c | 529 ------------------------------------- libcrystfel/src/diffraction-gpu.h | 57 ---- libcrystfel/src/diffraction.c | 463 --------------------------------- libcrystfel/src/diffraction.h | 34 --- libcrystfel/src/list_tmp.h | 106 -------- libcrystfel/src/peaks.c | 1 - libcrystfel/src/reflist-utils.c | 70 ----- libcrystfel/src/reflist-utils.h | 4 - libcrystfel/src/utils.h | 26 -- src/diffraction-gpu.c | 530 ++++++++++++++++++++++++++++++++++++++ src/diffraction-gpu.h | 57 ++++ src/diffraction.c | 464 +++++++++++++++++++++++++++++++++ src/diffraction.h | 34 +++ src/list_tmp.h | 106 ++++++++ src/pattern_sim.c | 71 +++++ src/pattern_sim.h | 50 ++++ 22 files changed, 1363 insertions(+), 1353 deletions(-) delete mode 100644 data/defs.h delete mode 100644 libcrystfel/src/diffraction-gpu.c delete mode 100644 libcrystfel/src/diffraction-gpu.h delete mode 100644 libcrystfel/src/diffraction.c delete mode 100644 libcrystfel/src/diffraction.h delete mode 100644 libcrystfel/src/list_tmp.h create mode 100644 src/diffraction-gpu.c create mode 100644 src/diffraction-gpu.h create mode 100644 src/diffraction.c create mode 100644 src/diffraction.h create mode 100644 src/list_tmp.h create mode 100644 src/pattern_sim.h diff --git a/Makefile.am b/Makefile.am index 5e7d8d19..f0e5053d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -35,12 +35,14 @@ LDADD += $(top_builddir)/libcrystfel/libcrystfel.la src_partial_sim_SOURCES = src/partial_sim.c -src_pattern_sim_SOURCES = src/pattern_sim.c +src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \ + src/diffraction-gpu.c if HAVE_OPENCL TESTS += tests/gpu_sim_check noinst_PROGRAMS += tests/gpu_sim_check -tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c +tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/diffraction.c \ + src/diffraction-gpu.c endif src_process_hkl_SOURCES = src/process_hkl.c @@ -88,7 +90,7 @@ EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \ src/post-refinement.h src/hrs-scaling.h src/scaling-report.h crystfeldir = $(datadir)/crystfel -crystfel_DATA = data/diffraction.cl data/defs.h data/hdfsee.ui +crystfel_DATA = data/diffraction.cl data/hdfsee.ui EXTRA_DIST += $(crystfel_DATA) diff --git a/Makefile.in b/Makefile.in index a1d869fc..87142b72 100644 --- a/Makefile.in +++ b/Makefile.in @@ -142,7 +142,8 @@ src_partialator_OBJECTS = $(am_src_partialator_OBJECTS) src_partialator_LDADD = $(LDADD) src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la -am_src_pattern_sim_OBJECTS = src/pattern_sim.$(OBJEXT) +am_src_pattern_sim_OBJECTS = src/pattern_sim.$(OBJEXT) \ + src/diffraction.$(OBJEXT) src/diffraction-gpu.$(OBJEXT) src_pattern_sim_OBJECTS = $(am_src_pattern_sim_OBJECTS) src_pattern_sim_LDADD = $(LDADD) src_pattern_sim_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ @@ -163,9 +164,12 @@ src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS) src_render_hkl_LDADD = $(LDADD) src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la -am__tests_gpu_sim_check_SOURCES_DIST = tests/gpu_sim_check.c +am__tests_gpu_sim_check_SOURCES_DIST = tests/gpu_sim_check.c \ + src/diffraction.c src/diffraction-gpu.c @HAVE_OPENCL_TRUE@am_tests_gpu_sim_check_OBJECTS = \ -@HAVE_OPENCL_TRUE@ tests/gpu_sim_check.$(OBJEXT) +@HAVE_OPENCL_TRUE@ tests/gpu_sim_check.$(OBJEXT) \ +@HAVE_OPENCL_TRUE@ src/diffraction.$(OBJEXT) \ +@HAVE_OPENCL_TRUE@ src/diffraction-gpu.$(OBJEXT) tests_gpu_sim_check_OBJECTS = $(am_tests_gpu_sim_check_OBJECTS) tests_gpu_sim_check_LDADD = $(LDADD) tests_gpu_sim_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ @@ -567,8 +571,12 @@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib \ LDADD = $(top_builddir)/lib/libgnu.a @IGNORE_UNUSED_LIBRARIES_CFLAGS@ \ $(top_builddir)/libcrystfel/libcrystfel.la src_partial_sim_SOURCES = src/partial_sim.c -src_pattern_sim_SOURCES = src/pattern_sim.c -@HAVE_OPENCL_TRUE@tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c +src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \ + src/diffraction-gpu.c + +@HAVE_OPENCL_TRUE@tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/diffraction.c \ +@HAVE_OPENCL_TRUE@ src/diffraction-gpu.c + src_process_hkl_SOURCES = src/process_hkl.c src_indexamajig_SOURCES = src/indexamajig.c @BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c @@ -588,7 +596,7 @@ tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \ INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data crystfeldir = $(datadir)/crystfel -crystfel_DATA = data/diffraction.cl data/defs.h data/hdfsee.ui +crystfel_DATA = data/diffraction.cl data/hdfsee.ui man_MANS = doc/man/indexamajig.1 doc/man/process_hkl.1 doc/man/pattern_sim.1 \ doc/man/crystfel_geometry.1 doc/man/powder_plot.1 @@ -777,6 +785,10 @@ src/partialator$(EXEEXT): $(src_partialator_OBJECTS) $(src_partialator_DEPENDENC $(AM_V_CCLD)$(LINK) $(src_partialator_OBJECTS) $(src_partialator_LDADD) $(LIBS) src/pattern_sim.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/diffraction.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/diffraction-gpu.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) src/pattern_sim$(EXEEXT): $(src_pattern_sim_OBJECTS) $(src_pattern_sim_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/pattern_sim$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_pattern_sim_OBJECTS) $(src_pattern_sim_LDADD) $(LIBS) @@ -832,6 +844,8 @@ mostlyclean-compile: -rm -f src/calibrate_detector.$(OBJEXT) -rm -f src/check_hkl.$(OBJEXT) -rm -f src/compare_hkl.$(OBJEXT) + -rm -f src/diffraction-gpu.$(OBJEXT) + -rm -f src/diffraction.$(OBJEXT) -rm -f src/dw-hdfsee.$(OBJEXT) -rm -f src/get_hkl.$(OBJEXT) -rm -f src/hdfsee.$(OBJEXT) @@ -857,6 +871,8 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/calibrate_detector.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/check_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/compare_hkl.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction-gpu.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dw-hdfsee.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/get_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/hdfsee.Po@am__quote@ diff --git a/data/defs.h b/data/defs.h deleted file mode 100644 index 088c533a..00000000 --- a/data/defs.h +++ /dev/null @@ -1,23 +0,0 @@ -/* - * defs.h - * - * Constant definitions which must be consistent throughout - * - * (c) 2006-2010 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifndef DEFS_H -#define DEFS_H - - -/* Maxmimum index to hold values up to (can be increased if necessary) */ -#define INDMAX 140 - -/* Array size */ -#define IDIM (INDMAX*2 +1) - - -#endif /* DEFS_H */ diff --git a/data/diffraction.cl b/data/diffraction.cl index 9f939413..0a463007 100644 --- a/data/diffraction.cl +++ b/data/diffraction.cl @@ -10,7 +10,11 @@ */ -#include +/* Maxmimum index to hold values up to (can be increased if necessary) + * WARNING: Altering this value constitutes an ABI change, and means you must + * update src/pattern_sim.h then recompile and reinstall everything. */ +#define INDMAX 140 + #ifndef M_PI #define M_PI ((float)(3.14159265)) #endif diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 576044a6..141842ff 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -1,22 +1,22 @@ lib_LTLIBRARIES = libcrystfel.la + libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/thread-pool.c src/image.c src/hdf5-file.c \ src/beam-parameters.c src/geometry.c src/statistics.c \ src/symmetry.c src/stream.c src/peaks.c \ - src/reflist-utils.c src/filters.c src/diffraction.c \ - src/diffraction-gpu.c src/cl-utils.c src/render.c \ - src/index.c src/dirax.c src/mosflm.c src/reax.c + src/reflist-utils.c src/filters.c src/cl-utils.c \ + src/render.c src/index.c src/dirax.c src/mosflm.c \ + src/reax.c libcrystfel_la_includedir=$(includedir)/crystfel/ -libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ - src/hdf5-file.h src/reflist.h src/symmetry.h \ - src/cell.h src/diffraction.h src/image.h \ + +libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ + src/reflist.h src/symmetry.h src/cell.h \ src/reflist-utils.h src/thread-pool.h \ - src/cl-utils.h src/filters.h src/list_tmp.h \ + src/cl-utils.h src/filters.h \ src/statistics.h src/utils.h src/detector.h \ src/geometry.h src/peaks.h src/stream.h \ - src/render.h src/index.h \ - ../data/defs.h + src/render.h src/index.h src/image.h INCLUDES = "-I$(top_srcdir)/data" AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib diff --git a/libcrystfel/Makefile.in b/libcrystfel/Makefile.in index a0b527e4..8992f9ed 100644 --- a/libcrystfel/Makefile.in +++ b/libcrystfel/Makefile.in @@ -88,9 +88,9 @@ am_libcrystfel_la_OBJECTS = src/reflist.lo src/utils.lo src/cell.lo \ src/detector.lo src/thread-pool.lo src/image.lo \ src/hdf5-file.lo src/beam-parameters.lo src/geometry.lo \ src/statistics.lo src/symmetry.lo src/stream.lo src/peaks.lo \ - src/reflist-utils.lo src/filters.lo src/diffraction.lo \ - src/diffraction-gpu.lo src/cl-utils.lo src/render.lo \ - src/index.lo src/dirax.lo src/mosflm.lo src/reax.lo + src/reflist-utils.lo src/filters.lo src/cl-utils.lo \ + src/render.lo src/index.lo src/dirax.lo src/mosflm.lo \ + src/reax.lo libcrystfel_la_OBJECTS = $(am_libcrystfel_la_OBJECTS) AM_V_lt = $(am__v_lt_$(V)) am__v_lt_ = $(am__v_lt_$(AM_DEFAULT_VERBOSITY)) @@ -351,20 +351,18 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/thread-pool.c src/image.c src/hdf5-file.c \ src/beam-parameters.c src/geometry.c src/statistics.c \ src/symmetry.c src/stream.c src/peaks.c \ - src/reflist-utils.c src/filters.c src/diffraction.c \ - src/diffraction-gpu.c src/cl-utils.c src/render.c \ - src/index.c src/dirax.c src/mosflm.c src/reax.c + src/reflist-utils.c src/filters.c src/cl-utils.c \ + src/render.c src/index.c src/dirax.c src/mosflm.c \ + src/reax.c libcrystfel_la_includedir = $(includedir)/crystfel/ -libcrystfel_la_include_HEADERS = src/beam-parameters.h src/diffraction-gpu.h \ - src/hdf5-file.h src/reflist.h src/symmetry.h \ - src/cell.h src/diffraction.h src/image.h \ +libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ + src/reflist.h src/symmetry.h src/cell.h \ src/reflist-utils.h src/thread-pool.h \ - src/cl-utils.h src/filters.h src/list_tmp.h \ + src/cl-utils.h src/filters.h \ src/statistics.h src/utils.h src/detector.h \ src/geometry.h src/peaks.h src/stream.h \ - src/render.h src/index.h \ - ../data/defs.h + src/render.h src/index.h src/image.h INCLUDES = "-I$(top_srcdir)/data" AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib \ @@ -457,9 +455,6 @@ src/peaks.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/reflist-utils.lo: src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/filters.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/diffraction.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/diffraction-gpu.lo: src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) src/cl-utils.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/render.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/index.lo: src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) @@ -479,10 +474,6 @@ mostlyclean-compile: -rm -f src/cl-utils.lo -rm -f src/detector.$(OBJEXT) -rm -f src/detector.lo - -rm -f src/diffraction-gpu.$(OBJEXT) - -rm -f src/diffraction-gpu.lo - -rm -f src/diffraction.$(OBJEXT) - -rm -f src/diffraction.lo -rm -f src/dirax.$(OBJEXT) -rm -f src/dirax.lo -rm -f src/filters.$(OBJEXT) @@ -525,8 +516,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cell.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/cl-utils.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/detector.Plo@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction-gpu.Plo@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/diffraction.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dirax.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/filters.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/geometry.Plo@am__quote@ diff --git a/libcrystfel/src/diffraction-gpu.c b/libcrystfel/src/diffraction-gpu.c deleted file mode 100644 index 605b1514..00000000 --- a/libcrystfel/src/diffraction-gpu.c +++ /dev/null @@ -1,529 +0,0 @@ -/* - * diffraction-gpu.c - * - * Calculate diffraction patterns by Fourier methods (GPU version) - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - -#ifdef HAVE_CL_CL_H -#include -#else -#include -#endif - -#include "image.h" -#include "utils.h" -#include "cell.h" -#include "diffraction.h" -#include "cl-utils.h" -#include "beam-parameters.h" - - -#define SAMPLING (4) -#define BWSAMPLING (10) -#define DIVSAMPLING (1) -#define SINC_LUT_ELEMENTS (4096) - - -struct gpu_context -{ - cl_context ctx; - cl_command_queue cq; - cl_program prog; - cl_kernel kern; - cl_mem intensities; - cl_mem flags; - - /* Array of sinc LUTs */ - cl_mem *sinc_luts; - cl_float **sinc_lut_ptrs; - int max_sinc_lut; /* Number of LUTs, i.e. one greater than the maximum - * index. This equals the highest allowable "n". */ -}; - - -static void check_sinc_lut(struct gpu_context *gctx, int n) -{ - cl_int err; - cl_image_format fmt; - int i; - - if ( n > gctx->max_sinc_lut ) { - - gctx->sinc_luts = realloc(gctx->sinc_luts, - n*sizeof(*gctx->sinc_luts)); - gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs, - n*sizeof(*gctx->sinc_lut_ptrs)); - - for ( i=gctx->max_sinc_lut; isinc_lut_ptrs[i] = NULL; - } - - gctx->max_sinc_lut = n; - } - - if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return; - - /* Create a new sinc LUT */ - gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float)); - gctx->sinc_lut_ptrs[n-1][0] = n; - if ( n == 1 ) { - for ( i=1; isinc_lut_ptrs[n-1][i] = 1.0; - } - } else { - for ( i=1; isinc_lut_ptrs[n-1][i] = val; - } - } - - fmt.image_channel_order = CL_INTENSITY; - fmt.image_channel_data_type = CL_FLOAT; - - gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - &fmt, SINC_LUT_ELEMENTS, 1, 0, - gctx->sinc_lut_ptrs[n-1], &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create LUT for %i\n", n); - return; - } -} - - -static int set_arg_float(struct gpu_context *gctx, int idx, float val) -{ - cl_int err; - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -static int set_arg_int(struct gpu_context *gctx, int idx, int val) -{ - cl_int err; - - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val) -{ - cl_int err; - - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell) -{ - cl_int err; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - float klow, khigh; - int i; - cl_float16 cell; - cl_int4 ncells; - const int sampling = SAMPLING; - cl_float bwstep; - int n_inf = 0; - int n_neg = 0; - cl_float divxlow, divxstep; - cl_float divylow, divystep; - int n_nan = 0; - int sprod; - - if ( gctx == NULL ) { - ERROR("GPU setup failed.\n"); - return; - } - - cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az; - cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz; - cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz; - - /* Calculate wavelength */ - klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0)); - khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0)); - bwstep = (khigh-klow) / BWSAMPLING; - - /* Calculate divergence stuff */ - divxlow = -image->beam->divergence/2.0; - divylow = -image->beam->divergence/2.0; - divxstep = image->beam->divergence / DIVSAMPLING; - divystep = image->beam->divergence / DIVSAMPLING; - - ncells.s[0] = na; - ncells.s[1] = nb; - ncells.s[2] = nc; - ncells.s[3] = 0; /* unused */ - - /* Ensure all required LUTs are available */ - check_sinc_lut(gctx, na); - check_sinc_lut(gctx, nb); - check_sinc_lut(gctx, nc); - - if ( set_arg_float(gctx, 2, klow) ) return; - if ( set_arg_mem(gctx, 9, gctx->intensities) ) return; - if ( set_arg_int(gctx, 12, sampling) ) return; - if ( set_arg_float(gctx, 14, bwstep) ) return; - if ( set_arg_mem(gctx, 15, gctx->sinc_luts[na-1]) ) return; - if ( set_arg_mem(gctx, 16, gctx->sinc_luts[nb-1]) ) return; - if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return; - if ( set_arg_mem(gctx, 18, gctx->flags) ) return; - if ( set_arg_float(gctx, 23, divxlow) ) return; - if ( set_arg_float(gctx, 24, divxstep) ) return; - if ( set_arg_int(gctx, 25, DIVSAMPLING) ) return; - if ( set_arg_float(gctx, 26, divylow) ) return; - if ( set_arg_float(gctx, 27, divystep) ) return; - if ( set_arg_int(gctx, 28, DIVSAMPLING) ) return; - - /* Unit cell */ - err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set unit cell: %s\n", clError(err)); - return; - } - - /* Local memory for reduction */ - sprod = BWSAMPLING*SAMPLING*SAMPLING*DIVSAMPLING*DIVSAMPLING; - err = clSetKernelArg(gctx->kern, 13, sprod*sizeof(cl_float), NULL); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set local memory: %s\n", clError(err)); - return; - } - - /* Allocate memory for the result */ - image->data = calloc(image->width * image->height, sizeof(float)); - image->twotheta = calloc(image->width * image->height, sizeof(double)); - - /* Iterate over panels */ - for ( i=0; idet->n_panels; i++ ) { - - size_t dims[3]; - size_t ldims[3] = {SAMPLING, SAMPLING, - BWSAMPLING * DIVSAMPLING * DIVSAMPLING}; - struct panel *p; - cl_mem tt; - size_t tt_size; - cl_mem diff; - size_t diff_size; - float *diff_ptr; - float *tt_ptr; - int pan_width, pan_height; - int fs, ss; - - p = &image->det->panels[i]; - - pan_width = 1 + p->max_fs - p->min_fs; - pan_height = 1 + p->max_ss - p->min_ss; - - /* Buffer for the results of this panel */ - diff_size = pan_width * pan_height * sizeof(cl_float); - diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, - diff_size, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate diffraction memory\n"); - return; - } - tt_size = pan_width * pan_height * sizeof(cl_float); - tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size, - NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate twotheta memory\n"); - return; - } - - if ( set_arg_mem(gctx, 0, diff) ) return; - if ( set_arg_mem(gctx, 1, tt) ) return; - if ( set_arg_int(gctx, 3, pan_width) ) return; - if ( set_arg_float(gctx, 4, p->cnx) ) return; - if ( set_arg_float(gctx, 5, p->cny) ) return; - if ( set_arg_float(gctx, 6, p->res) ) return; - if ( set_arg_float(gctx, 7, p->clen) ) return; - if ( set_arg_int(gctx, 10, p->min_fs) ) return; - if ( set_arg_int(gctx, 11, p->min_ss) ) return; - if ( set_arg_float(gctx, 19, p->fsx) ) return; - if ( set_arg_float(gctx, 20, p->fsy) ) return; - if ( set_arg_float(gctx, 21, p->ssx) ) return; - if ( set_arg_float(gctx, 22, p->ssy) ) return; - - dims[0] = pan_width * SAMPLING; - dims[1] = pan_height * SAMPLING; - dims[2] = BWSAMPLING * DIVSAMPLING * DIVSAMPLING; - - err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 3, NULL, - dims, ldims, 0, NULL, NULL); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't enqueue diffraction kernel: %s\n", - clError(err)); - return; - } - - clFinish(gctx->cq); - - diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE, - CL_MAP_READ, 0, diff_size, - 0, NULL, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't map diffraction buffer: %s\n", - clError(err)); - return; - } - tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ, - 0, tt_size, 0, NULL, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't map tt buffer\n"); - return; - } - - for ( fs=0; fsmin_fs + fs; - tss = p->min_ss + ss; - image->data[tfs + image->width*tss] = val; - image->twotheta[tfs + image->width*tss] = tt; - - } - } - - clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr, - 0, NULL, NULL); - clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr, - 0, NULL, NULL); - - clReleaseMemObject(diff); - clReleaseMemObject(tt); - - } - - - if ( n_neg + n_inf + n_nan ) { - ERROR("WARNING: The GPU calculation produced %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg, n_inf, n_nan); - } - -} - - -/* Setup the OpenCL stuff, create buffers, load the structure factor table */ -struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, unsigned char *flags, - const char *sym, int dev_num) -{ - struct gpu_context *gctx; - cl_uint nplat; - cl_platform_id platforms[8]; - cl_context_properties prop[3]; - cl_int err; - cl_device_id dev; - size_t intensities_size; - float *intensities_ptr; - size_t flags_size; - float *flags_ptr; - size_t maxwgsize; - int i; - char cflags[512] = ""; - - STATUS("Setting up GPU...\n"); - - err = clGetPlatformIDs(8, platforms, &nplat); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't get platform IDs: %i\n", err); - return NULL; - } - if ( nplat == 0 ) { - ERROR("Couldn't find at least one platform!\n"); - return NULL; - } - prop[0] = CL_CONTEXT_PLATFORM; - prop[1] = (cl_context_properties)platforms[0]; - prop[2] = 0; - - gctx = malloc(sizeof(*gctx)); - gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU, - NULL, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create OpenCL context: %i\n", err); - free(gctx); - return NULL; - } - - dev = get_cl_dev(gctx->ctx, dev_num); - - gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create OpenCL command queue\n"); - free(gctx); - return NULL; - } - - /* Create a single-precision version of the scattering factors */ - intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float); - intensities_ptr = malloc(intensities_size); - if ( intensities != NULL ) { - for ( i=0; iintensities = clCreateBuffer(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - intensities_size, intensities_ptr, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate intensities memory\n"); - free(gctx); - return NULL; - } - free(intensities_ptr); - - if ( sym != NULL ) { - if ( strcmp(sym, "1") == 0 ) { - strncat(cflags, "-DPG1 ", 511-strlen(cflags)); - } else if ( strcmp(sym, "-1") == 0 ) { - strncat(cflags, "-DPG1BAR ", 511-strlen(cflags)); - } else if ( strcmp(sym, "6/mmm") == 0 ) { - strncat(cflags, "-DPG6MMM ", 511-strlen(cflags)); - } 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 { - 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)); - } - } else { - if ( intensities != NULL ) { - ERROR("You gave me an intensities file but no point" - " group. I'm assuming '1'.\n"); - strncat(cflags, "-DPG1 ", 511-strlen(cflags)); - } - } - - /* Create a flag array */ - flags_size = IDIM*IDIM*IDIM*sizeof(cl_float); - flags_ptr = malloc(flags_size); - if ( flags != NULL ) { - for ( i=0; iflags = clCreateBuffer(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - flags_size, flags_ptr, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate flag buffer\n"); - free(gctx); - return NULL; - } - free(flags_ptr); - - gctx->prog = load_program(DATADIR"/crystfel/diffraction.cl", gctx->ctx, - dev, &err, cflags); - if ( err != CL_SUCCESS ) { - free(gctx); - return NULL; - } - - gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create kernel\n"); - free(gctx); - return NULL; - } - - gctx->max_sinc_lut = 0; - gctx->sinc_lut_ptrs = NULL; - gctx->sinc_luts = NULL; - - clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE, - sizeof(size_t), &maxwgsize, NULL); - STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize); - - return gctx; -} - - -void cleanup_gpu(struct gpu_context *gctx) -{ - int i; - - clReleaseProgram(gctx->prog); - clReleaseMemObject(gctx->intensities); - - /* Release LUTs */ - for ( i=1; i<=gctx->max_sinc_lut; i++ ) { - if ( gctx->sinc_lut_ptrs[i-1] != NULL ) { - clReleaseMemObject(gctx->sinc_luts[i-1]); - free(gctx->sinc_lut_ptrs[i-1]); - } - } - - free(gctx->sinc_luts); - free(gctx->sinc_lut_ptrs); - - clReleaseCommandQueue(gctx->cq); - clReleaseContext(gctx->ctx); - free(gctx); -} diff --git a/libcrystfel/src/diffraction-gpu.h b/libcrystfel/src/diffraction-gpu.h deleted file mode 100644 index a3bde4e1..00000000 --- a/libcrystfel/src/diffraction-gpu.h +++ /dev/null @@ -1,57 +0,0 @@ -/* - * diffraction-gpu.h - * - * Calculate diffraction patterns by Fourier methods (GPU version) - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#ifndef DIFFRACTION_GPU_H -#define DIFFRACTION_GPU_H - -#include "image.h" -#include "cell.h" - -struct gpu_context; - -#if HAVE_OPENCL - -extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell); -extern struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, - const unsigned char *flags, - const char *sym, int dev_num); -extern void cleanup_gpu(struct gpu_context *gctx); - -#else - -static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell) -{ - /* Do nothing */ - ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n"); -} - -static struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, - const unsigned char *flags, - const char *sym, int dev_num) -{ - return NULL; -} - -static void cleanup_gpu(struct gpu_context *gctx) -{ -} - -#endif - -#endif /* DIFFRACTION_GPU_H */ diff --git a/libcrystfel/src/diffraction.c b/libcrystfel/src/diffraction.c deleted file mode 100644 index 9532a6ce..00000000 --- a/libcrystfel/src/diffraction.c +++ /dev/null @@ -1,463 +0,0 @@ -/* - * diffraction.c - * - * Calculate diffraction patterns by Fourier methods - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#include -#include -#include -#include -#include -#include -#include - -#include "image.h" -#include "utils.h" -#include "cell.h" -#include "diffraction.h" -#include "beam-parameters.h" -#include "symmetry.h" - - -#define SAMPLING (4) -#define BWSAMPLING (10) -#define DIVSAMPLING (1) -#define SINC_LUT_ELEMENTS (4096) - - -static double *get_sinc_lut(int n) -{ - int i; - double *lut; - - lut = malloc(SINC_LUT_ELEMENTS*sizeof(double)); - lut[0] = n; - if ( n == 1 ) { - for ( i=1; i= 0.0); - - val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); - val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); - - val1 = val1; - val2 = val2; - - return (1.0-f)*val1 + f*val2; -} - - -static double interpolate_bilinear(const double *ref, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, signed int l) -{ - signed int k; - double val1, val2; - float f; - - k = (signed int)kd; - if ( kd < 0.0 ) k -= 1; - f = kd - (float)k; - assert(f >= 0.0); - - val1 = interpolate_linear(ref, flags, sym, hd, k, l); - val2 = interpolate_linear(ref, flags, sym, hd, k+1, l); - - return (1.0-f)*val1 + f*val2; -} - - -static double interpolate_intensity(const double *ref, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, float ld) -{ - signed int l; - double val1, val2; - float f; - - l = (signed int)ld; - if ( ld < 0.0 ) l -= 1; - f = ld - (float)l; - assert(f >= 0.0); - - val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l); - val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1); - - return (1.0-f)*val1 + f*val2; -} - - -static double complex interpolate_phased_linear(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, - signed int k, signed int l) -{ - signed int h; - double val1, val2; - float f; - double ph1, ph2; - double re1, re2, im1, im2; - double re, im; - - h = (signed int)hd; - if ( hd < 0.0 ) h -= 1; - f = hd - (float)h; - assert(f >= 0.0); - - val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); - val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); - ph1 = sym_lookup_phase(phases, flags, sym, h, k, l); - ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l); - - val1 = val1; - val2 = val2; - - /* Calculate real and imaginary parts */ - re1 = val1 * cos(ph1); - im1 = val1 * sin(ph1); - re2 = val2 * cos(ph2); - im2 = val2 * sin(ph2); - - re = (1.0-f)*re1 + f*re2; - im = (1.0-f)*im1 + f*im2; - - return re + im*I; -} - - -static double complex interpolate_phased_bilinear(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, - signed int l) -{ - signed int k; - double complex val1, val2; - float f; - - k = (signed int)kd; - if ( kd < 0.0 ) k -= 1; - f = kd - (float)k; - assert(f >= 0.0); - - val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l); - val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l); - - return (1.0-f)*val1 + f*val2; -} - - -static double interpolate_phased_intensity(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, float ld) -{ - signed int l; - double complex val1, val2; - float f; - - l = (signed int)ld; - if ( ld < 0.0 ) l -= 1; - f = ld - (float)l; - assert(f >= 0.0); - - val1 = interpolate_phased_bilinear(ref, phases, flags, sym, - hd, kd, l); - val2 = interpolate_phased_bilinear(ref, phases, flags, sym, - hd, kd, l+1); - - return cabs((1.0-f)*val1 + f*val2); -} - - -/* Look up the structure factor for the nearest Bragg condition */ -static double molecule_factor(const double *intensities, const double *phases, - const unsigned char *flags, struct rvec q, - double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - GradientMethod m, const SymOpList *sym) -{ - float hd, kd, ld; - signed int h, k, l; - double r; - - hd = q.u * ax + q.v * ay + q.w * az; - kd = q.u * bx + q.v * by + q.w * bz; - ld = q.u * cx + q.v * cy + q.w * cz; - - /* No flags -> flat intensity distribution */ - if ( flags == NULL ) return 1.0e5; - - switch ( m ) { - case GRADIENT_MOSAIC : - fesetround(1); /* Round to nearest */ - h = (signed int)rint(hd); - k = (signed int)rint(kd); - l = (signed int)rint(ld); - if ( abs(h) > INDMAX ) r = 0.0; - else if ( abs(k) > INDMAX ) r = 0.0; - else if ( abs(l) > INDMAX ) r = 0.0; - else r = sym_lookup_intensity(intensities, flags, sym, h, k, l); - break; - case GRADIENT_INTERPOLATE : - r = interpolate_intensity(intensities, flags, sym, hd, kd, ld); - break; - case GRADIENT_PHASED : - r = interpolate_phased_intensity(intensities, phases, flags, - sym, hd, kd, ld); - break; - default: - ERROR("This gradient method not implemented yet.\n"); - exit(1); - } - - return r; -} - - -void get_diffraction(struct image *image, int na, int nb, int nc, - const double *intensities, const double *phases, - const unsigned char *flags, UnitCell *cell, - GradientMethod m, const SymOpList *sym) -{ - unsigned int fs, ss; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - float klow, khigh, bwstep; - double *lut_a; - double *lut_b; - double *lut_c; - double divxlow, divylow, divxstep, divystep; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - /* Allocate (and zero) the "diffraction array" */ - image->data = calloc(image->width * image->height, sizeof(float)); - - /* Needed later for Lorentz calculation */ - image->twotheta = malloc(image->width * image->height * sizeof(double)); - - klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0)); - khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0)); - bwstep = (khigh-klow) / BWSAMPLING; - - divxlow = -image->beam->divergence/2.0; - divylow = -image->beam->divergence/2.0; - divxstep = image->beam->divergence / DIVSAMPLING; - divystep = image->beam->divergence / DIVSAMPLING; - - lut_a = get_sinc_lut(na); - lut_b = get_sinc_lut(nb); - lut_c = get_sinc_lut(nc); - - for ( fs=0; fswidth; fs++ ) { - for ( ss=0; ssheight; ss++ ) { - - int fs_step, ss_step, kstep; - int divxval, divyval; - int idx = fs + image->width*ss; - - for ( fs_step=0; fs_stepdata[idx] += intensity; - - if ( fs_step + ss_step + kstep == 0 ) { - image->twotheta[idx] = twotheta; - } - - } - } - } - } - } - - image->data[idx] /= (SAMPLING*SAMPLING*BWSAMPLING - *DIVSAMPLING*DIVSAMPLING); - - - } - progress_bar(fs, image->width-1, "Calculating diffraction"); - } - - free(lut_a); - free(lut_b); - free(lut_c); -} diff --git a/libcrystfel/src/diffraction.h b/libcrystfel/src/diffraction.h deleted file mode 100644 index f71d3cce..00000000 --- a/libcrystfel/src/diffraction.h +++ /dev/null @@ -1,34 +0,0 @@ -/* - * diffraction.h - * - * Calculate diffraction patterns by Fourier methods - * - * (c) 2006-2011 Thomas White - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#ifndef DIFFRACTION_H -#define DIFFRACTION_H - -#include "image.h" -#include "cell.h" -#include "symmetry.h" - -typedef enum { - GRADIENT_MOSAIC, - GRADIENT_INTERPOLATE, - GRADIENT_PHASED -} GradientMethod; - -extern void get_diffraction(struct image *image, int na, int nb, int nc, - const double *intensities, const double *phases, - const unsigned char *flags, UnitCell *cell, - GradientMethod m, const SymOpList *sym); - -#endif /* DIFFRACTION_H */ diff --git a/libcrystfel/src/list_tmp.h b/libcrystfel/src/list_tmp.h deleted file mode 100644 index a524b2f9..00000000 --- a/libcrystfel/src/list_tmp.h +++ /dev/null @@ -1,106 +0,0 @@ -/* - * Template for creating indexed 3D lists of a given type, usually indexed - * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX. - * - * These are used, for example, for: - * - a list of 'double complex' values for storing structure factors, - * - a list of 'double' values for storing reflection intensities, - * - a list of 'unsigned int' values for counts of some sort. - * - * When LABEL and TYPE are #defined appropriately, including this header - * defines functions such as: - * - new_list_