Changeset 15711

Timestamp:

Nov 28, 2007, 4:56:01 PM (18 years ago)

Author:

Paul Price

Message:

Removing heap sort code that was straight out of Numerical Recipes.
Replacing with heap sort code from GNU Scientific Library,
reformatted. Advantage of reformatting is that the direct sort and
index sort are now implemented using the same code, with different
functions for comparison and swapping. Code passes tests in
tap_psVectorSort and tap_psVectorSortIndex.

File:

• trunk/psLib/src/mathtypes/psVector.c (modified) (4 diffs)

trunk/psLib/src/mathtypes/psVector.c

@@ -10 +10 @@
 *  @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 +344 @@
 // 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 +416 @@
     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 +471 @@
 }
 
-// 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;
 }