Index: /trunk/psLib/src/mathtypes/psVector.c =================================================================== --- /trunk/psLib/src/mathtypes/psVector.c (revision 15710) +++ /trunk/psLib/src/mathtypes/psVector.c (revision 15711) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt, University of Hawaii * -* @version $Revision: 1.99 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-10-29 20:35:31 $ +* @version $Revision: 1.100 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-11-29 02:56:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -344,35 +344,71 @@ // XXX since the key size is fixed (same number of bits) a Radix sort could be // a big win here -- someone should benchmark this -JH -#define PSVECTOR_SORT_CASE(NAME) \ -case PS_TYPE_##NAME: { \ - ps##NAME temp; \ - ps##NAME *value = vector->data.NAME; \ - for (;;) { \ - if (l > 0) { \ - temp = value[--l]; \ + +// The following sort code is based on gsl_heapsort from GSL-1.8 (www.gnu.org/software/gsl), file +// gsl-1.8/sort/sort.c, which is distributed under the GNU General Public License, version 2. +// +// Implement Heap sort -- direct and indirect sorting +// Based on descriptions in Sedgewick "Algorithms in C" +// +// Copyright (C) 1999 Thomas Walter +// +// 18 February 2000: Modified for GSL by Brian Gough + +// Comparison and swap functions for sorting values directly +#define PSVECTOR_SORT_COMPARE_DIRECT(A,B) (value[A] < value[B]) +#define PSVECTOR_SORT_SWAP_DIRECT(TYPE,A,B) { \ + if (A != B) { \ + ps##TYPE temp = value[A]; \ + value[A] = value[B]; \ + value[B] = temp; \ + } \ +} + +// Comparison and swap functions for sorting vector indices +#define PSVECTOR_SORT_COMPARE_INDEX(A,B) (value[index[A]] < value[index[B]]) +#define PSVECTOR_SORT_SWAP_INDEX(TYPE,A,B) { \ + if (A != B) { \ + ps##TYPE temp = index[A]; \ + index[A] = index[B]; \ + index[B] = temp; \ + } \ +} + +// Sort the heap +#define PSVECTOR_SORT_DOWNHEAP(SWAPTYPE, COMPAREFUNC, SWAPFUNC, K) { \ + unsigned long k = K; /* Local version of k --- so the higher-level version isn't modified */ \ + while (k <= N / 2) { \ + unsigned long j = 2 * k; \ + if (j < N && COMPAREFUNC(j, j + 1)) { \ + j++; \ + } \ + if (COMPAREFUNC(k, j)) { \ + SWAPFUNC(SWAPTYPE, j, k); \ } else { \ - temp = value[ir]; \ - value[ir] = value[0]; \ - if (--ir == 0) { \ - value[0] = temp; \ - return true; \ - } \ + break; \ } \ - int i = l; \ - int j = (l << 1) + 1; \ - while (j <= ir) { \ - if (j < ir && value[j] < value[j+1]) \ - ++j; \ - if (temp < value[j]) { \ - value[i]=value[j]; \ - j += (i=j) + 1; \ - } else { \ - j = ir + 1; \ - } \ - } \ - value[i] = temp; \ + k = j; \ } \ -} \ -break; +} + +// Driver for heap sort +#define PSVECTOR_SORT_CASE(ELEMTYPE, SWAPTYPE, COMPAREFUNC, SWAPFUNC) \ +case PS_TYPE_##ELEMTYPE: { \ + ps##ELEMTYPE *value = vector->data.ELEMTYPE; \ + unsigned long N = vector->n - 1; /* Index of last element */ \ + unsigned long i = N / 2 + 1; /* Adding one to compensate for i-- below */ \ + do { \ + i--; \ + PSVECTOR_SORT_DOWNHEAP(SWAPTYPE, COMPAREFUNC, SWAPFUNC, i); \ + } while (i > 0); \ + while (N > 0) { \ + SWAPFUNC(SWAPTYPE, 0, N); /* Swap elements */ \ + /* Process the heap */ \ + N--; \ + PSVECTOR_SORT_DOWNHEAP(SWAPTYPE, COMPAREFUNC, SWAPFUNC, 0); \ + } \ + break; \ +} +// END of heap sort code from GSL bool psVectorSortInPlace(const psVector *vector) @@ -380,21 +416,20 @@ PS_ASSERT_VECTOR_NON_NULL(vector, false); - long N = vector->n; // Number of elements - if (N < 2) { + if (vector->n < 2) { + // Already sorted! return true; } - long l = N >> 1; - long ir = N - 1; + switch (vector->type.type) { - PSVECTOR_SORT_CASE(U8); - PSVECTOR_SORT_CASE(U16); - PSVECTOR_SORT_CASE(U32); - PSVECTOR_SORT_CASE(U64); - PSVECTOR_SORT_CASE(S8); - PSVECTOR_SORT_CASE(S16); - PSVECTOR_SORT_CASE(S32); - PSVECTOR_SORT_CASE(S64); - PSVECTOR_SORT_CASE(F32); - PSVECTOR_SORT_CASE(F64); + PSVECTOR_SORT_CASE(U8 , U8 , PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(U16, U16, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(U32, U32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(U64, U64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(S8 , S8 , PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(S16, S16, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(S32, S32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(S64, S64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(F32, F32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); + PSVECTOR_SORT_CASE(F64, F64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT); default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, @@ -436,76 +471,36 @@ } -// Heap sort of the index array -#define PSVECTOR_SORT_INDEX_CASE(NAME) \ -case PS_TYPE_##NAME: { \ - psS32 temp; \ - psS32 *index = out->data.S32; \ - ps##NAME *value = inVector->data.NAME; \ - for (;;) { \ - if (l > 0) { \ - temp = index[--l]; \ - } else { \ - temp = index[ir]; \ - index[ir] = index[0]; \ - if (--ir == 0) { \ - index[0] = temp; \ - return out; \ - } \ - } \ - long i = l; \ - long j = (l << 1) + 1; \ - while (j <= ir) { \ - if (j < ir && value[index[j]] < value[index[j+1]]) \ - ++j; \ - if (value[temp ] < value[index[j]]) { \ - index[i]=index[j]; \ - j += (i=j) + 1; \ - } else { \ - j = ir + 1; \ - } \ - } \ - index[i] = temp; \ - } \ -} \ -break; - -psVector* psVectorSortIndex(psVector* outVector, - const psVector* inVector) -{ - if (inVector == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("psVectorSort can not sort a NULL psVector.")); - psFree(outVector); - return NULL; - } - - // XXX isn't this line a waste of space? - psVector *out = psVectorRecycle(outVector, inVector->n, PS_TYPE_S32); // Vector for output - out = psVectorCreate(out, 0, inVector->n, 1, PS_TYPE_S32); - long N = out->n; // Number of elements - if (N < 2) { - return out; - } - long l = N >> 1; - long ir = N - 1; - switch (inVector->type.type) { - PSVECTOR_SORT_INDEX_CASE(U8); - PSVECTOR_SORT_INDEX_CASE(U16); - PSVECTOR_SORT_INDEX_CASE(U32); - PSVECTOR_SORT_INDEX_CASE(U64); - PSVECTOR_SORT_INDEX_CASE(S8); - PSVECTOR_SORT_INDEX_CASE(S16); - PSVECTOR_SORT_INDEX_CASE(S32); - PSVECTOR_SORT_INDEX_CASE(S64); - PSVECTOR_SORT_INDEX_CASE(F32); - PSVECTOR_SORT_INDEX_CASE(F64); +psVector* psVectorSortIndex(psVector *out, + const psVector *vector) +{ + PS_ASSERT_VECTOR_NON_NULL(vector, NULL); + + psVector *indexVector = psVectorCreate(out, 0, vector->n, 1, PS_TYPE_S32); // Array of indices + psS32 *index = indexVector->data.S32; // Dereference for convenience + + if (vector->n < 2) { + // Already sorted + return indexVector; + } + switch (vector->type.type) { + PSVECTOR_SORT_CASE(U8 , S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(U16, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(U32, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(U64, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(S8 , S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(S16, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(S32, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(S64, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(F32, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); + PSVECTOR_SORT_CASE(F64, S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX); default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, _("Input psVector is an unsupported type (0x%x)."), - inVector->type.type); - psFree(outVector); - return NULL; - } - return out; + vector->type.type); + psFree(indexVector); + return NULL; + } + + return indexVector; }