added qsort_r replacement

darcs-hash:20081109061514-c8de0-209e83f287c63329b8e0a609c8eb237cbc0ada29.gz

diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c
index 20b7650fd68237050c3cc2d2250fb812c005d663..a09b3f47824c24215701aa0b669fb2a7e933d18a 100644 (file)
--- 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 0f50a2d0fd0bee9d42eedd12c9083e5f40e5c0f3..9ec9e71de9fef233d560f5db38379df5ccaddc8e 100644 (file)
--- 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 <math.h>

diff --git a/util/Makefile.am b/util/Makefile.am
index fe6e1591579755a39fe03e55eb01f3297b63a50e..41b501fc5897f16bbbb7c527f9ec08c153156f5d 100644 (file)
--- 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 6cb45f8f99909ef7114735fa9ce6f7dc0fcd6a6c..46e2ca6a16bab123074182b8eb30ec25fc842568 100644 (file)
--- a/util/nlopt-util.h
+++ b/util/nlopt-util.h
@@ -23,6 +23,8 @@
 #ifndef NLOPT_UTIL_H
 #define NLOPT_UTIL_H
 
+#include <stdlib.h>
+
 #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 (file)
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 <nlopt-util.h>
+#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 */
+}