From fa1b4e33a15894870229eb594f5e159a840c5fb6 Mon Sep 17 00:00:00 2001 From: stevenj Date: Sun, 9 Nov 2008 01:15:14 -0500 Subject: [PATCH] added qsort_r replacement darcs-hash:20081109061514-c8de0-209e83f287c63329b8e0a609c8eb237cbc0ada29.gz --- cdirect/cdirect.c | 12 +++--- configure.ac | 2 +- util/Makefile.am | 2 +- util/nlopt-util.h | 6 +++ util/qsort_r.c | 103 ++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 116 insertions(+), 9 deletions(-) create mode 100644 util/qsort_r.c diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c index 20b7650..a09b3f4 100644 --- a/cdirect/cdirect.c +++ b/cdirect/cdirect.c @@ -114,12 +114,12 @@ static double rect_diameter(int n, const double *w, const params *p) #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY -static double *fv_qsort = 0; -static int sort_fv_compare(const void *a_, const void *b_) +static int sort_fv_compare(void *fv_, const void *a_, const void *b_) { + const double *fv = (const double *) fv_; int a = *((const int *) a_), b = *((const int *) b_); - double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]); - double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]); + double fa = MIN(fv[2*a], fv[2*a+1]); + double fb = MIN(fv[2*b], fv[2*b+1]); if (fa < fb) return -1; else if (fa > fb) @@ -131,9 +131,7 @@ static void sort_fv(int n, double *fv, int *isort) { int i; for (i = 0; i < n; ++i) isort[i] = i; - fv_qsort = fv; /* not re-entrant, sigh... */ - qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare); - fv_qsort = 0; + nlopt_qsort_r(isort, (unsigned) n, sizeof(int), fv, sort_fv_compare); } static double function_eval(const double *x, params *p) { diff --git a/configure.ac b/configure.ac index 0f50a2d..9ec9e71 100644 --- a/configure.ac +++ b/configure.ac @@ -53,7 +53,7 @@ AC_CHECK_TYPES(uint32_t, [], [], [$ac_includes_default dnl Checks for libraries and functions AC_CHECK_LIB(m, sin) -AC_CHECK_FUNCS([BSDgettimeofday gettimeofday time]) +AC_CHECK_FUNCS([BSDgettimeofday gettimeofday time qsort_r]) AC_MSG_CHECKING([for isnan]) AC_TRY_LINK([#include diff --git a/util/Makefile.am b/util/Makefile.am index fe6e159..41b501f 100644 --- a/util/Makefile.am +++ b/util/Makefile.am @@ -1,5 +1,5 @@ noinst_LTLIBRARIES = libutil.la -libutil_la_SOURCES = mt19937ar.c sobolseq.c soboldata.h timer.c stop.c nlopt-util.h redblack.c redblack.h +libutil_la_SOURCES = mt19937ar.c sobolseq.c soboldata.h timer.c stop.c nlopt-util.h redblack.c redblack.h qsort_r.c noinst_PROGRAMS = redblack_test redblack_test_SOURCES = redblack_test.c diff --git a/util/nlopt-util.h b/util/nlopt-util.h index 6cb45f8..46e2ca6 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -23,6 +23,8 @@ #ifndef NLOPT_UTIL_H #define NLOPT_UTIL_H +#include + #ifdef __cplusplus extern "C" { @@ -30,6 +32,10 @@ extern "C" int nlopt_isinf(double x); +/* re-entrant qsort */ +extern void nlopt_qsort_r(void *base_, size_t nmemb, size_t size, void *thunk, + int (*compar)(void *, const void *, const void *)); + /* seconds timer */ extern double nlopt_seconds(void); extern unsigned long nlopt_time_seed(void); diff --git a/util/qsort_r.c b/util/qsort_r.c new file mode 100644 index 0000000..583cdc5 --- /dev/null +++ b/util/qsort_r.c @@ -0,0 +1,103 @@ +/* Copyright (c) 2007-2008 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include "config.h" + +/* Simple replacement for the BSD qsort_r function (re-entrant sorting), + if it is not available. (It looks like glibc 2.8 will include qsort_r + as well.) + + (Actually, with glibc 2.3.6 on my Intel Core Duo, my implementation + below seems to be significantly faster than qsort. Go figure. + Even so, I'd rather use the libc version of qsort_r if it is available, + on general principles, and to reduce code bloat.) */ + +#ifndef HAVE_QSORT_R +/* swap size bytes between a_ and b_ */ +static void swap(void *a_, void *b_, size_t size) +{ + if (a_ == b_) return; + { + size_t i, nlong = size / sizeof(long); + long *a = (long *) a_, *b = (long *) b_; + for (i = 0; i < nlong; ++i) { + long c = a[i]; + a[i] = b[i]; + b[i] = c; + } + a_ = (void*) (a + nlong); + b_ = (void*) (b + nlong); + } + { + size_t i; + char *a = (char *) a_, *b = (char *) b_; + size = size % sizeof(long); + for (i = 0; i < size; ++i) { + char c = a[i]; + a[i] = b[i]; + b[i] = c; + } + } +} +#endif /* HAVE_QSORT_R */ + +void nlopt_qsort_r(void *base_, size_t nmemb, size_t size, void *thunk, + int (*compar)(void *, const void *, const void *)) +{ +#ifdef HAVE_QSORT_R + qsort_r(base_, nmemb, size, thunk, compar); +#else + char *base = (char *) base_; + if (nmemb < 10) { /* use O(nmemb^2) algorithm for small enough nmemb */ + size_t i, j; + for (i = 0; i+1 < nmemb; ++i) + for (j = i+1; j < nmemb; ++j) + if (compar(thunk, base+i*size, base+j*size) > 0) + swap(base+i*size, base+j*size, size); + } + else { + size_t i, pivot, npart; + /* pick median of first/middle/last elements as pivot */ + { + const char *a = base, *b = base + (nmemb/2)*size, + *c = base + (nmemb-1)*size; + pivot = compar(thunk,a,b) < 0 + ? (compar(thunk,b,c) < 0 ? nmemb/2 : + (compar(thunk,a,c) < 0 ? nmemb-1 : 0)) + : (compar(thunk,a,c) < 0 ? 0 : + (compar(thunk,b,c) < 0 ? nmemb-1 : nmemb/2)); + } + /* partition array */ + swap(base + pivot*size, base + (nmemb-1) * size, size); + pivot = (nmemb - 1) * size; + for (i = npart = 0; i < nmemb-1; ++i) + if (compar(thunk, base+i*size, base+pivot) <= 0) + swap(base+i*size, base+(npart++)*size, size); + swap(base+npart*size, base+pivot, size); + /* recursive sort of two partitions */ + nlopt_qsort_r(base, npart, size, thunk, compar); + npart++; /* don't need to sort pivot */ + nlopt_qsort_r(base+npart*size, nmemb-npart, size, thunk, compar); + } +#endif /* !HAVE_QSORT_R */ +} -- 2.30.2