Skip to content
GitLab
Commit 03ccff08 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added option KDTEST to test original kdtree sort function.

parent 00ad30aa
Hide whitespace changes
Inline Side-by-side
src/Makefile.am
View file @ 03ccff08
......@@ -150,6 +150,7 @@ libcdo_la_SOURCES += \
kdtreelib/kdtree_cartesian.c \
kdtreelib/kdtree_common.c \
kdtreelib/kdtree_spherical.c \
kdtreelib/qsort.c \
kdtreelib/pmergesort.c \
kdtreelib/pqueue.c \
kdtreelib/pqueue.h
......
src/Makefile.in
View file @ 03ccff08
......@@ -144,7 +144,7 @@ am_libcdo_la_OBJECTS = libcdo_la-cdo_pthread.lo libcdo_la-cdo_vlist.lo \
kdtreelib/libcdo_la-kdtree_cartesian.lo \
kdtreelib/libcdo_la-kdtree_common.lo \
kdtreelib/libcdo_la-kdtree_spherical.lo \
kdtreelib/libcdo_la-pmergesort.lo \
kdtreelib/libcdo_la-qsort.lo kdtreelib/libcdo_la-pmergesort.lo \
kdtreelib/libcdo_la-pqueue.lo clipping/libcdo_la-clipping.lo \
clipping/libcdo_la-area.lo \
clipping/libcdo_la-ensure_array_size.lo \
......@@ -528,12 +528,13 @@ libcdo_la_SOURCES = cdo_int.h compare.h cdo_pthread.c cdo_vlist.c \
table.c text.c text.h timebase.h timer.c userlog.c uthash.h \
util.c util.h zaxis.c kdtreelib/kdtree.h \
kdtreelib/kdtree_cartesian.c kdtreelib/kdtree_common.c \
kdtreelib/kdtree_spherical.c kdtreelib/pmergesort.c \
kdtreelib/pqueue.c kdtreelib/pqueue.h clipping/clipping.c \
clipping/clipping.h clipping/area.c clipping/area.h \
clipping/ensure_array_size.c clipping/ensure_array_size.h \
clipping/geometry.h clipping/grid.h clipping/points.h \
clipping/dep_list.h clipping/grid_cell.c clipping/grid_cell.h \
kdtreelib/kdtree_spherical.c kdtreelib/qsort.c \
kdtreelib/pmergesort.c kdtreelib/pqueue.c kdtreelib/pqueue.h \
clipping/clipping.c clipping/clipping.h clipping/area.c \
clipping/area.h clipping/ensure_array_size.c \
clipping/ensure_array_size.h clipping/geometry.h \
clipping/grid.h clipping/points.h clipping/dep_list.h \
clipping/grid_cell.c clipping/grid_cell.h \
clipping/intersection.c clipping/utils.c clipping/utils.h
#
#if ENABLE_NEARPT3
......@@ -675,6 +676,8 @@ kdtreelib/libcdo_la-kdtree_common.lo: kdtreelib/$(am__dirstamp) \
kdtreelib/$(DEPDIR)/$(am__dirstamp)
kdtreelib/libcdo_la-kdtree_spherical.lo: kdtreelib/$(am__dirstamp) \
kdtreelib/$(DEPDIR)/$(am__dirstamp)
kdtreelib/libcdo_la-qsort.lo: kdtreelib/$(am__dirstamp) \
kdtreelib/$(DEPDIR)/$(am__dirstamp)
kdtreelib/libcdo_la-pmergesort.lo: kdtreelib/$(am__dirstamp) \
kdtreelib/$(DEPDIR)/$(am__dirstamp)
kdtreelib/libcdo_la-pqueue.lo: kdtreelib/$(am__dirstamp) \
......@@ -1067,6 +1070,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@kdtreelib/$(DEPDIR)/libcdo_la-kdtree_spherical.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@kdtreelib/$(DEPDIR)/libcdo_la-pmergesort.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@kdtreelib/$(DEPDIR)/libcdo_la-pqueue.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@kdtreelib/$(DEPDIR)/libcdo_la-qsort.Plo@am__quote@
.c.o:
@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
......@@ -1694,6 +1698,13 @@ kdtreelib/libcdo_la-kdtree_spherical.lo: kdtreelib/kdtree_spherical.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o kdtreelib/libcdo_la-kdtree_spherical.lo `test -f 'kdtreelib/kdtree_spherical.c' || echo '$(srcdir)/'`kdtreelib/kdtree_spherical.c
kdtreelib/libcdo_la-qsort.lo: kdtreelib/qsort.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT kdtreelib/libcdo_la-qsort.lo -MD -MP -MF kdtreelib/$(DEPDIR)/libcdo_la-qsort.Tpo -c -o kdtreelib/libcdo_la-qsort.lo `test -f 'kdtreelib/qsort.c' || echo '$(srcdir)/'`kdtreelib/qsort.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) kdtreelib/$(DEPDIR)/libcdo_la-qsort.Tpo kdtreelib/$(DEPDIR)/libcdo_la-qsort.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='kdtreelib/qsort.c' object='kdtreelib/libcdo_la-qsort.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o kdtreelib/libcdo_la-qsort.lo `test -f 'kdtreelib/qsort.c' || echo '$(srcdir)/'`kdtreelib/qsort.c
kdtreelib/libcdo_la-pmergesort.lo: kdtreelib/pmergesort.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT kdtreelib/libcdo_la-pmergesort.lo -MD -MP -MF kdtreelib/$(DEPDIR)/libcdo_la-pmergesort.Tpo -c -o kdtreelib/libcdo_la-pmergesort.lo `test -f 'kdtreelib/pmergesort.c' || echo '$(srcdir)/'`kdtreelib/pmergesort.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) kdtreelib/$(DEPDIR)/libcdo_la-pmergesort.Tpo kdtreelib/$(DEPDIR)/libcdo_la-pmergesort.Plo
......
src/kdtreelib/kdtree_common.c
View file @ 03ccff08
......@@ -8,7 +8,11 @@
extern int pmergesort(void *base, size_t nmemb, size_t size,
#if defined(KDTEST)
int (*compar) (const void *, const void *, int), int axis,
#else
int (*compar) (const void *, const void *),
#endif
int max_threads);
......@@ -87,6 +91,18 @@ kd_printTree(struct kdNode *node)
******************************************************************** */
static int
_compPoints(const void *p1, const void *p2, int axis)
{
struct kd_point *a = (struct kd_point *) p1;
struct kd_point *b = (struct kd_point *) p2;
int ret = (a->point[axis] > b->point[axis]) ? 1 : (a->point[axis] < b->point[axis]) ? -1 : 0;
if ( ret == 0 ) ret = (a->index > b->index) ? 1 : (a->index < b->index) ? -1 : 0;
return ret;
}
static int
_compPoints0(const void *p1, const void *p2)
{
......@@ -189,7 +205,11 @@ kd_doBuildTree(void *threadarg)
else if ( sortaxis == 1 ) qcomp = _compPoints1;
else qcomp = _compPoints2;
#if defined(KDTEST)
pmergesort(points, nPoints, sizeof(struct kd_point), _compPoints, sortaxis, max_threads);
#else
pmergesort(points, nPoints, sizeof(struct kd_point), qcomp, max_threads);
#endif
pivot = nPoints / 2;
if ((node = kd_allocNode(points, pivot, min, max, sortaxis, dim)) == NULL)
......
src/kdtreelib/pmergesort.c
View file @ 03ccff08
......@@ -10,24 +10,44 @@ typedef struct param_t {
size_t first;
size_t nmemb;
size_t size;
#if defined(KDTEST)
int axis;
int (*cmp) (const void *, const void *, int axis);
#else
int (*cmp) (const void *, const void *);
#endif
int max_threads;
} param_t;
extern void qsortR(const void *base0, size_t n, size_t size,
int (*compar) (const void *, const void *, int axis),
int axis);
void pm_buildparams(struct param_t *p, void *a, void *b, size_t first,
size_t nmemb, size_t size,
#if defined(KDTEST)
int (*cmp) (const void *, const void *, int), int axis,
#else
int (*cmp) (const void *, const void *),
#endif
int max_threads);
int pmergesort(void *base, size_t nmemb, size_t size,
#if defined(KDTEST)
int (*compar) (const void *, const void *, int), int axis,
#else
int (*compar) (const void *, const void *),
#endif
int max_threads);
void *mergesort_t(void *args);
int
pmergesort(void *base, size_t nmemb, size_t size,
#if defined(KDTEST)
int (*compar) (const void *, const void *, int), int axis,
#else
int (*compar) (const void *, const void *),
#endif
int max_threads)
{
void *tmp;
......@@ -43,6 +63,9 @@ pmergesort(void *base, size_t nmemb, size_t size,
args.nmemb = nmemb;
args.size = size;
args.cmp = compar;
#if defined(KDTEST)
args.axis = axis;
#endif
args.max_threads = max_threads;
mergesort_t(&args);
......@@ -54,7 +77,11 @@ pmergesort(void *base, size_t nmemb, size_t size,
void
pm_buildparams(struct param_t *p, void *a, void *b, size_t first,
size_t nmemb, size_t size,
#if defined(KDTEST)
int (*cmp) (const void *, const void *, int), int axis,
#else
int (*cmp) (const void *, const void *),
#endif
int max_threads)
{
......@@ -64,6 +91,9 @@ pm_buildparams(struct param_t *p, void *a, void *b, size_t first,
p->nmemb = nmemb;
p->size = size;
p->cmp = cmp;
#if defined(KDTEST)
p->axis = axis;
#endif
p->max_threads = max_threads;
}
......@@ -81,13 +111,21 @@ mergesort_t(void *args)
* Reached maximum number of threads allocated to this
* branch. Proceed with sequential sort of this chunk.
*/
#if defined(KDTEST)
qsortR(mya->a + mya->first * mya->size, mya->nmemb, mya->size, mya->cmp, mya->axis);
#else
qsort(mya->a + mya->first * mya->size, mya->nmemb, mya->size, mya->cmp);
#endif
} else {
/*
* Start two new threads, each sorting half of array a
*/
pm_buildparams(&larg, mya->a, mya->b, mya->first, mya->nmemb / 2,
mya->size, mya->cmp, mya->max_threads / 2);
mya->size, mya->cmp,
#if defined(KDTEST)
mya->axis,
#endif
mya->max_threads / 2);
/*
* Recursively sort the left half
*/
......@@ -97,8 +135,11 @@ mergesort_t(void *args)
}
pm_buildparams(&rarg, mya->a, mya->b, mya->first + mya->nmemb / 2,
mya->nmemb - mya->nmemb / 2, mya->size,
mya->cmp, mya->max_threads / 2);
mya->nmemb - mya->nmemb / 2, mya->size, mya->cmp,
#if defined(KDTEST)
mya->axis,
#endif
mya->max_threads / 2);
/*
* Recursively sort the right half
*/
......@@ -138,7 +179,11 @@ mergesort_t(void *args)
* element
*/
else if (mya->cmp(mya->a + li * mya->size,
mya->a + ri * mya->size) < 1) {
mya->a + ri * mya->size
#if defined(KDTEST)
, mya->axis
#endif
) < 1) {
memcpy(mya->b + i * mya->size, mya->a + li * mya->size,
mya->size);
li++;
......
src/kdtreelib/qsort.c 0 → 100644
View file @ 03ccff08
/* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details (qsort.c
is used to build libc.a and as such is covered by the LGPL -- Joerg
Dietrich 2008)*/
#include <stdlib.h>
/*-
* Copyright (c) 1980, 1983 The Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted
* provided that: (1) source distributions retain this entire copyright
* notice and comment, and (2) distributions including binaries display
* the following acknowledgement: ``This product includes software
* developed by the University of California, Berkeley and its contributors''
* in the documentation or other materials provided with the distribution
* and in all advertising materials mentioning features or use of this
* software. Neither the name of the University nor the names of its
* contributors may be used to endorse or promote products derived
* from this software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
/*
* qsort.c:
* Our own version of the system qsort routine which is faster by an average
* of 25%, with lows and highs of 10% and 50%.
* The THRESHold below is the insertion sort threshold, and has been adjusted
* for records of size 48 bytes.
* The MTHREShold is where we stop finding a better median.
*/
#define THRESH 4 /* threshold for insertion */
#define MTHRESH 6 /* threshold for median */
static int (*qcmp) (const void *, const void *, int axis); /* the comparison routine */
static int qsz; /* size of each record */
static int thresh; /* THRESHold in chars */
static int mthresh; /* MTHRESHold in chars */
void qsortR(const void *base0, size_t n, size_t size,
int (*compar) (const void *, const void *, int), int axis);
/*
* qst:
* Do a quicksort
* First, find the median element, and put that one in the first place as the
* discriminator. (This "median" is just the median of the first, last and
* middle elements). (Using this median instead of the first element is a big
* win). Then, the usual partitioning/swapping, followed by moving the
* discriminator into the right place. Then, figure out the sizes of the two
* partions, do the smaller one recursively and the larger one via a repeat of
* this code. Stopping when there are less than THRESH elements in a partition
* and cleaning up with an insertion sort (in our caller) is a huge win.
* All data swaps are done in-line, which is space-losing but time-saving.
* (And there are only three places where this is done).
*/
static void
qst(char *base, char *max, int axis)
{
char c, *i, *j, *jj;
int ii;
char *mid, *tmp;
int lo, hi;
/*
* At the top here, lo is the number of characters of elements in the
* current partition. (Which should be max - base).
* Find the median of the first, last, and middle element and make
* that the middle element. Set j to largest of first and middle.
* If max is larger than that guy, then it's that guy, else compare
* max with loser of first and take larger. Things are set up to
* prefer the middle, then the first in case of ties.
*/
lo = max - base; /* number of elements as chars */
do {
mid = i = base + qsz * ((lo / qsz) >> 1);
if (lo >= mthresh) {
j = (qcmp((jj = base), i, axis) > 0 ? jj : i);
if (qcmp(j, (tmp = max - qsz), axis) > 0) {
/*
* switch to first loser
*/
j = (j == jj ? i : jj);
if (qcmp(j, tmp, axis) < 0)
j = tmp;
}
if (j != i) {
ii = qsz;
do {
c = *i;
*i++ = *j;
*j++ = c;
} while (--ii);
}
}
/*
* Semi-standard quicksort partitioning/swapping
*/
for(i = base, j = max - qsz;;) {
while (i < mid && qcmp(i, mid, axis) <= 0)
i += qsz;
while (j > mid) {
if (qcmp(mid, j, axis) <= 0) {
j -= qsz;
continue;
}
tmp = i + qsz; /* value of i after swap */
if (i == mid) {
/*
* j <-> mid, new mid is j
*/
mid = jj = j;
} else {
/*
* i <-> j
*/
jj = j;
j -= qsz;
}
goto swap;
}
if (i == mid) {
break;
} else {
/*
* i <-> mid, new mid is i
*/
jj = mid;
tmp = mid = i; /* value of i after swap */
j -= qsz;
}
swap:
ii = qsz;
do {
c = *i;
*i++ = *jj;
*jj++ = c;
} while (--ii);
i = tmp;
}
/*
* Look at sizes of the two partitions, do the smaller
* one first by recursion, then do the larger one by
* making sure lo is its size, base and max are update
* correctly, and branching back. But only repeat
* (recursively or by branching) if the partition is
* of at least size THRESH.
*/
i = (j = mid) + qsz;
if ((lo = j - base) <= (hi = max - i)) {
if (lo >= thresh)
qst(base, j, axis);
base = i;
lo = hi;
} else {
if (hi >= thresh)
qst(i, max, axis);
max = j;
}
} while (lo >= thresh);
}
/*
* qsort:
* First, set up some global parameters for qst to share. Then, quicksort
* with qst(), and then a cleanup insertion sort ourselves. Sound simple?
* It's not...
*/
void
qsortR(const void *base0, size_t n, size_t size,
int (*compar) (const void *, const void *, int axis), int axis)
{
char *base = (char *) base0;
char c, *i, *j, *lo, *hi;
char *min, *max;
if (n <= 1)
return;
qsz = size;
qcmp = compar;
thresh = qsz * THRESH;
mthresh = qsz * MTHRESH;
max = base + n * qsz;
if (n >= THRESH) {
qst(base, max, axis);
hi = base + thresh;
} else {
hi = max;
}
/*
* First put smallest element, which must be in the first THRESH, in
* the first position as a sentinel. This is done just by searching
* the first THRESH elements (or the first n if n < THRESH), finding
* the min, and swapping it into the first position.
*/
for(j = lo = base; (lo += qsz) < hi;)
if (qcmp(j, lo, axis) > 0)
j = lo;
if (j != base) {
/*
* swap j into place
*/
for(i = base, hi = base + qsz; i < hi;) {
c = *j;
*j++ = *i;
*i++ = c;
}
}
/*
* With our sentinel in place, we now run the following hyper-fast
* insertion sort. For each remaining element, min, from [1] to [n-1],
* set hi to the index of the element AFTER which this one goes.
* Then, do the standard insertion sort shift on a character at a time
* basis for each element in the frob.
*/
for(min = base; (hi = min += qsz) < max;) {
while (qcmp(hi -= qsz, min, axis) > 0)
/*
* void
*/ ;
if ((hi += qsz) != min) {
for(lo = min + qsz; --lo >= min;) {
c = *lo;
for(i = j = lo; (j -= qsz) >= hi; i = j)
*i = *j;
*i = c;
}
}
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Imprint | Privacy Policy