Commit 76a6d396 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

kdtree: added new version for qsortR().

parent ef648758
......@@ -50,15 +50,18 @@ typedef struct kd_point {
static inline
int qcmp(const void *p1, const void *p2, int axis)
int qcmp(struct kd_point *a, struct kd_point *b, 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;
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;
}
return ret;
static inline
int qcmp2(struct kd_point *a, struct kd_point *b, int axis)
{
return (a->point[axis] < b->point[axis]) ? -1 : 0;
}
......
......@@ -127,27 +127,6 @@ _compPoints2(const void *p1, const void *p2)
return ret;
}
/*
void kd_quick_sort(struct kd_point *a, int n, int idx)
{
int i, j;
struct kd_point p, t;
if ( n < 2 ) return;
p = a[n / 2];
for ( i = 0, j = n - 1;; i++, j-- )
{
while ( a[i].point[idx] < p.point[idx] ) i++;
while ( p.point[idx] < a[j].point[idx] ) j--;
if ( i >= j ) break;
t = a[i];
a[i] = a[j];
a[j] = t;
}
kd_quick_sort(a, i, idx);
kd_quick_sort(a + i, n - i, idx);
}
*/
void *
kd_doBuildTree(void *threadarg)
{
......
......@@ -20,7 +20,7 @@ typedef struct param_t {
int max_threads;
} param_t;
extern void qsortR(const void *base0, size_t n, size_t size, int axis);
extern void qsortR(void *base0, size_t n, int axis);
void pm_buildparams(struct param_t *p, void *a, void *b, size_t first,
size_t nmemb, size_t size,
......@@ -113,7 +113,7 @@ mergesort_t(void *args)
* branch. Proceed with sequential sort of this chunk.
*/
#if defined(KDTEST)
qsortR(mya->a + mya->first * mya->size, mya->nmemb, mya->size, mya->axis);
qsortR(mya->a + mya->first * mya->size, mya->nmemb, mya->axis);
#else
qsort(mya->a + mya->first * mya->size, mya->nmemb, mya->size, mya->cmp);
#endif
......@@ -185,8 +185,8 @@ mergesort_t(void *args)
*/
else if (
#if defined(KDTEST)
qcmp(mya->a + li * mya->size,
mya->a + ri * mya->size,
qcmp((struct kd_point *)(mya->a + li * mya->size),
(struct kd_point *)(mya->a + ri * mya->size),
mya->axis)
#else
mya->cmp(mya->a + li * mya->size,
......
......@@ -36,11 +36,10 @@
#define THRESH 4 /* threshold for insertion */
#define MTHRESH 6 /* threshold for median */
static int qsz; /* size of each record */
static int thresh; /* THRESHold in chars */
static int mthresh; /* MTHRESHold in chars */
static int thresh; /* THRESHold */
static int mthresh; /* MTHRESHold */
void qsortR(const void *base0, size_t n, size_t size, int axis);
void qsortR(void *base0, size_t n, int axis);
/*
* qst:
......@@ -57,110 +56,122 @@ void qsortR(const void *base0, size_t n, size_t size, int axis);
* (And there are only three places where this is done).
*/
static void
qst(char *base, char *max, int axis)
static
void qst(struct kd_point *base, struct kd_point *max, int axis)
{
char c, *i, *j, *jj;
int ii;
char *mid, *tmp;
int lo, hi;
struct kd_point c, *i, *j, *jj;
int ii;
struct kd_point *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 */
do {
mid = i = base + (lo >> 1);
if ( lo >= mthresh )
{
j = (qcmp((jj = base), i, axis) > 0 ? jj : i);
if ( qcmp(j, (tmp = max - 1), axis) > 0 )
{
/*
* switch to first loser
*/
j = (j == jj ? i : jj);
if ( qcmp(j, tmp, axis) < 0 ) j = tmp;
}
if ( j != i )
{
ii = 1;
do {
c = *i;
*i++ = *j;
*j++ = c;
} while (--ii);
}
}
/*
* 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.
* Semi-standard quicksort partitioning/swapping
*/
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) {
for ( i = base, j = max - 1;; )
{
while ( i < mid && qcmp(i, mid, axis) <= 0 ) ++i;
while ( j > mid )
{
if ( qcmp(mid, j, axis) <= 0 )
{
--j;
continue;
}
tmp = i + 1; /* value of i after swap */
if ( i == mid )
{
/*
* switch to first loser
* j <-> mid, new mid is j
*/
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 {
mid = jj = j;
}
else
{
/*
* i <-> mid, new mid is i
* i <-> j
*/
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);
jj = j;
--j;
}
goto swap;
}
if ( i == mid )
{
break;
}
else
{
/*
* i <-> mid, new mid is i
*/
jj = mid;
tmp = mid = i; /* value of i after swap */
--j;
}
swap:
ii = 1;
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) + 1;
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 );
}
/*
......@@ -170,63 +181,140 @@ qst(char *base, char *max, int axis)
* It's not...
*/
void
qsortR(const void *base0, size_t n, size_t size, int axis)
void XXXqsortR(void *base0, size_t n, int axis)
{
char *base = (char *) base0;
char c, *i, *j, *lo, *hi;
char *min, *max;
if (n <= 1)
return;
qsz = size;
thresh = qsz * THRESH;
mthresh = qsz * MTHRESH;
max = base + n * qsz;
if (n >= THRESH) {
qst(base, max, axis);
hi = base + thresh;
} else {
hi = max;
struct kd_point *base = (struct kd_point *) base0;
struct kd_point c, *i, *j, *lo, *hi;
struct kd_point *min, *max;
if ( n <= 1 ) return;
thresh = THRESH;
mthresh = MTHRESH;
max = base + n;
if ( n >= THRESH )
{
qst(base, max, axis);
hi = base + thresh;
}
/*
* 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) {
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 += 1) < hi; )
if ( qcmp(j, lo, axis) > 0 )
j = lo;
if ( j != base )
{
/*
* swap j into place
*/
for ( i = base, hi = base + 1; 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 += 1) < max; )
{
while ( qcmp(hi -= 1, min, axis) > 0 )
/*
* swap j into place
*/
for(i = base, hi = base + qsz; i < hi;) {
c = *j;
*j++ = *i;
*i++ = c;
* void
*/ ;
if ( (hi += 1) != min )
{
for ( lo = min + 1; --lo >= min; )
{
c = *lo;
for ( i = j = lo; (j -= 1) >= hi; i = j ) *i = *j;
*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;
}
// quickSort
//
// This public-domain C implementation by Darel Rex Finley.
//
// * This function assumes it is called with valid parameters.
//
// * Example calls:
// quickSort(&myArray[0],5); // sorts elements 0, 1, 2, 3, and 4
// quickSort(&myArray[3],5); // sorts elements 3, 4, 5, 6, and 7
void XqsortR(void *base0, size_t n, int axis)
{
struct kd_point *base = (struct kd_point *) base0;
struct kd_point piv;
#define MAX_LEVELS 300
int beg[MAX_LEVELS], end[MAX_LEVELS], L, R, swap, i = 0;
beg[0] = 0; end[0] = (int)n;
while ( i >= 0 )
{
L = beg[i]; R = end[i]-1;
if ( L < R )
{
piv = base[L];
while ( L < R )
{
while ( L < R && qcmp(&base[R], &piv, axis) >= 0 ) R--;
if ( L < R ) base[L++] = base[R];
while ( L < R && qcmp(&base[L], &piv, axis) <= 0 ) L++;
if ( L < R ) base[R--] = base[L];
}
base[L] = piv; beg[i+1] = L+1; end[i+1] = end[i]; end[i++] = L;
if ( end[i]-beg[i] > end[i-1]-beg[i-1] )
{
swap = beg[i]; beg[i] = beg[i-1]; beg[i-1] = swap;
swap = end[i]; end[i] = end[i-1]; end[i-1] = swap;
}
}
else
{
i--;
}
}
}
void qsortR(void *base0, size_t n, int idx)
{
if ( n < 2 ) return;
struct kd_point *a = (struct kd_point *) base0;
struct kd_point t;
kdata_t p = a[n / 2].point[idx];
int i, j;
for ( i = 0, j = n - 1;; i++, j-- )
{
while ( a[i].point[idx] < p ) i++;
while ( p < a[j].point[idx] ) j--;
if ( i >= j ) break;
t = a[i];
a[i] = a[j];
a[j] = t;
}
qsortR(a, i, idx);
qsortR(a + i, n - i, idx);
}
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