From bca7a90bc0bb8e5aecbbdd946e950bcaa24f1c77 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 29 Jul 2011 16:07:16 +0200 Subject: More ReAx stuff --- config.h.in | 3 ++ configure | 70 +++++++++++++++++++++++++++++++++- configure.ac | 15 +++++++- src/reax.c | 122 +++++++++++++++++++++++++++++++++++++++++++++++++++++++---- 4 files changed, 200 insertions(+), 10 deletions(-) diff --git a/config.h.in b/config.h.in index ea2b3689..b09f25d8 100644 --- a/config.h.in +++ b/config.h.in @@ -38,6 +38,9 @@ /* Define to 1 if you have the header file. */ #undef HAVE_FCNTL_H +/* Define to 1 if FFTW is available. */ +#undef HAVE_FFTW + /* Define to 1 if you have the `forkpty' function. */ #undef HAVE_FORKPTY diff --git a/configure b/configure index e47de008..569d2a6b 100755 --- a/configure +++ b/configure @@ -616,6 +616,8 @@ HTML_DIR GTKDOC_MKPDF GTKDOC_REBASE GTKDOC_CHECK +HAVE_FFTW_FALSE +HAVE_FFTW_TRUE HAVE_CAIRO_FALSE HAVE_CAIRO_TRUE BUILD_CUBEIT_FALSE @@ -7537,6 +7539,60 @@ fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for fftw_plan_dft_r2c_1d in -lfftw3" >&5 +$as_echo_n "checking for fftw_plan_dft_r2c_1d in -lfftw3... " >&6; } +if test "${ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d+set}" = set; then : + $as_echo_n "(cached) " >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lfftw3 $LIBS" +cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + +/* Override any GCC internal prototype to avoid an error. + Use char because int might match the return type of a GCC + builtin and then its argument prototype would still apply. */ +#ifdef __cplusplus +extern "C" +#endif +char fftw_plan_dft_r2c_1d (); +int +main () +{ +return fftw_plan_dft_r2c_1d (); + ; + return 0; +} +_ACEOF +if ac_fn_c_try_link "$LINENO"; then : + ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d=yes +else + ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d=no +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" >&5 +$as_echo "$ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" >&6; } +if test "x$ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" = x""yes; then : + + +$as_echo "#define HAVE_FFTW 1" >>confdefs.h + + FFTW_LIBS="-lfftw3" + have_fftw=true + +else + + { $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: ReAx indexing will not be available." >&5 +$as_echo "$as_me: WARNING: ReAx indexing will not be available." >&2;} + have_fftw=false + +fi + + + if test x$have_opencl = xtrue; then HAVE_OPENCL_TRUE= HAVE_OPENCL_FALSE='#' @@ -7576,6 +7632,14 @@ else fi + if test x$have_fftw = xtrue; then + HAVE_FFTW_TRUE= + HAVE_FFTW_FALSE='#' +else + HAVE_FFTW_TRUE='#' + HAVE_FFTW_FALSE= +fi + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for C compiler flag to ignore unused libraries" >&5 @@ -7642,7 +7706,7 @@ CFLAGS="$CFLAGS $GDK_pixbuf_CFLAGS $GDK_pixbuf_2_CFLAGS" LIBS="$LIBS $HDF5_LIBS -lm -lz $GSL_LIBS $GTK_LIBS $OPENCL_LIBS -pthread" LIBS="$LIBS $LIBTIFF_LIBS $libPNG_LIBS $Cairo_LIBS $GDK_pixbuf_LIBS" -LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $LDFLAGS" +LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $FFTW_LIBS $LDFLAGS" @@ -8043,6 +8107,10 @@ if test -z "${HAVE_CAIRO_TRUE}" && test -z "${HAVE_CAIRO_FALSE}"; then as_fn_error $? "conditional \"HAVE_CAIRO\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 fi +if test -z "${HAVE_FFTW_TRUE}" && test -z "${HAVE_FFTW_FALSE}"; then + as_fn_error $? "conditional \"HAVE_FFTW\" was never defined. +Usually this means the macro was only invoked conditionally." "$LINENO" 5 +fi if test -z "${ENABLE_GTK_DOC_TRUE}" && test -z "${ENABLE_GTK_DOC_FALSE}"; then as_fn_error $? "conditional \"ENABLE_GTK_DOC\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 diff --git a/configure.ac b/configure.ac index 1b7ccb70..01a89553 100644 --- a/configure.ac +++ b/configure.ac @@ -230,6 +230,18 @@ AC_CHECK_LIB([rt], [clock_gettime], ]) +AC_CHECK_LIB([fftw3], [fftw_plan_dft_r2c_1d], +[ + AC_DEFINE([HAVE_FFTW], [1], + [Define to 1 if FFTW is available.]) + FFTW_LIBS="-lfftw3" + have_fftw=true +], [ + AC_MSG_WARN([ReAx indexing will not be available.]) + have_fftw=false +]) + + dnl Conditionals... AM_CONDITIONAL([HAVE_OPENCL], test x$have_opencl = xtrue) @@ -242,6 +254,7 @@ AM_CONDITIONAL([BUILD_CUBEIT], test x$have_cairo = xtrue \ AM_CONDITIONAL([HAVE_CAIRO], test x$have_cairo = xtrue) +AM_CONDITIONAL([HAVE_FFTW], test x$have_fftw = xtrue) gl_IGNORE_UNUSED_LIBRARIES @@ -252,7 +265,7 @@ CFLAGS="$CFLAGS $GDK_pixbuf_CFLAGS $GDK_pixbuf_2_CFLAGS" LIBS="$LIBS $HDF5_LIBS -lm -lz $GSL_LIBS $GTK_LIBS $OPENCL_LIBS -pthread" LIBS="$LIBS $LIBTIFF_LIBS $libPNG_LIBS $Cairo_LIBS $GDK_pixbuf_LIBS" -LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $LDFLAGS" +LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $FFTW_LIBS $LDFLAGS" GTK_DOC_CHECK([1.11],[--flavour no-tmpl]) diff --git a/src/reax.c b/src/reax.c index 77ca124f..20e8c94c 100644 --- a/src/reax.c +++ b/src/reax.c @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include "image.h" #include "utils.h" @@ -39,9 +41,116 @@ struct reax_private { IndexingPrivate base; struct dvec *directions; + int n_dir; }; +static double check_dir(struct dvec *dir, ImageFeatureList *flist, double modv, + double pmax, double *fft_in, fftw_complex *fft_out, + fftw_plan plan) +{ + int n, i, v; + double *vals; + + n = image_feature_count(flist); + + vals = malloc(n*sizeof(double)); + v = 0; + for ( i=0; irx*dir->x + f->ry*dir->y + f->rz*dir->z; + + v++; /* Yuk, yuk, yuk */ + + } + + + + + free(vals); + + return 0.0; +} + + +void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +{ + int i; + struct reax_private *p; + double fom; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double mod_as; + double dx, dy, dz; + int nel, n; + double pmax; + double *fft_in; + fftw_complex *fft_out; + fftw_plan plan; + const double cellmax = 50.0e9; /* 50 nm max cell size */ + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + mod_as = modulus(asx, asy, asz); + + n = image_feature_count(image->features); + pmax = 0.0; + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + val = modulus(f->rx, f->ry, f->rz); + + if ( val > pmax ) pmax = val; + + } + nel = 2.0*pmax*5.0*cellmax; + + fft_in = fftw_malloc(nel*sizeof(double)); + fft_out = fftw_malloc((nel/2 + 1)*sizeof(fftw_complex)); + + plan = fftw_plan_dft_r2c_1d(nel, fft_in, fft_out, FFTW_ESTIMATE); + + /* Search for a* */ + fom = 0.0; dx = 0.0; dy = 0.0; dz = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom; + + new_fom = check_dir(&p->directions[i], image->features, mod_as, + pmax, fft_in, fft_out, plan); + if ( new_fom > fom ) { + fom = new_fom; + dx = p->directions[i].x; + dy = p->directions[i].x; + dz = p->directions[i].x; + } + + } + + fftw_destroy_plan(plan); + fftw_free(fft_in); + fftw_free(fft_out); + + /* No improvement from zero? */ + if ( fom == 0.0 ) return; +} + + IndexingPrivate *reax_prepare() { struct reax_private *priv; @@ -58,12 +167,14 @@ IndexingPrivate *reax_prepare() angular_inc = 0.03; /* From Steller (1997) */ samp = (2.0 * M_PI) / angular_inc; - priv->directions = malloc(samp*samp*sizeof(struct dvec)); + priv->n_dir = samp*samp; + priv->directions = malloc(priv->n_dir*sizeof(struct dvec)); if ( priv == NULL) { free(priv); return NULL; } + /* Generate vectors for 1D Fourier transforms */ for ( ui=0; uidirections[ui + vi*samp]; @@ -88,9 +200,3 @@ IndexingPrivate *reax_prepare() return (IndexingPrivate *)priv; } - - -void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell) -{ - -} -- cgit v1.2.3