Index: /trunk/psLib/src/math/Makefile.am =================================================================== --- /trunk/psLib/src/math/Makefile.am (revision 15712) +++ /trunk/psLib/src/math/Makefile.am (revision 15713) @@ -57,5 +57,6 @@ noinst_HEADERS = \ - md5.h + md5.h \ + psSort.h CLEANFILES = *~ *.bb *.bbg *.da Index: /trunk/psLib/src/math/psSort.h =================================================================== --- /trunk/psLib/src/math/psSort.h (revision 15713) +++ /trunk/psLib/src/math/psSort.h (revision 15713) @@ -0,0 +1,76 @@ +// Heap sort of the array +// Based on descriptions in Sedgewick "Algorithms in C" +// +// Copyright (C) 1999 Thomas Walter +// Copyright (C) 2007 Paul Price, Institute for Astronomy, University of Hawaii +// +// 18 February 2000: Modified for GSL by Brian Gough +// 29 November 2007: Modified for psLib by Paul Price +// +// +// This is free software; you can redistribute it and/or modify it +// under the terms of the GNU General Public License as published by the +// Free Software Foundation; either version 2, or (at your option) any +// later version. +// +// This source is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +// for more details. +// + +#ifndef PS_SORT_H +#define PS_SORT_H + +// 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 + +// 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. + + +// Use of the PSSORT() macro requires additional macros defined: +// COMPAREEXPR(A,B): expression to compare element A with element B; true if element A is smaller than B. +// SWAPFUNC(TYPE,A,B): swap element A with element B; the type is provided so that a temporary variable can +// be defined. + +// Sort the heap +#define PSSORT_DOWNHEAP(COMPAREFUNC, SWAPFUNC, SWAPTYPE, 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 { \ + break; \ + } \ + k = j; \ + } \ +} + +// Driver for heap sort +// SIZE: Size of the array +// COMPAREEXPR: Macro with expression for comparison +// SWAPFUNC: Macro to swap elements +// SWAPTYPE: Type for swapping, passed to SWAPFUNC +#define PSSORT(SIZE, COMPAREEXPR, SWAPFUNC, SWAPTYPE) { \ + unsigned long N = (SIZE) - 1; /* Index of last element */ \ + unsigned long i = N / 2 + 1; /* Adding one to compensate for i-- below */ \ + do { \ + i--; \ + PSSORT_DOWNHEAP(COMPAREEXPR, SWAPFUNC, SWAPTYPE, i); \ + } while (i > 0); \ + while (N > 0) { \ + SWAPFUNC(SWAPTYPE, 0, N); /* Swap elements */ \ + /* Process the heap */ \ + N--; \ + PSSORT_DOWNHEAP(COMPAREEXPR, SWAPFUNC, SWAPTYPE, 0); \ + } \ +} +// END of heap sort code from GSL + + +#endif Index: /trunk/psLib/src/mathtypes/psVector.c =================================================================== --- /trunk/psLib/src/mathtypes/psVector.c (revision 15712) +++ /trunk/psLib/src/mathtypes/psVector.c (revision 15713) @@ -10,6 +10,6 @@ * @author Joshua Hoblitt, University of Hawaii * -* @version $Revision: 1.100 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-11-29 02:56:01 $ +* @version $Revision: 1.101 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-11-29 21:46:54 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -33,4 +33,5 @@ #include "psAssert.h" #include "psString.h" +#include "psSort.h" static void vectorFree(psVector* psVec) @@ -341,17 +342,4 @@ } -// Heap sort of the array -// 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 - -// 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 @@ -375,40 +363,10 @@ } -// 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 { \ - break; \ - } \ - k = j; \ - } \ -} - -// Driver for heap sort -#define PSVECTOR_SORT_CASE(ELEMTYPE, SWAPTYPE, COMPAREFUNC, SWAPFUNC) \ +#define PSVECTOR_SORT_CASE(ELEMTYPE, COMPAREEXPR, SWAPFUNC, SWAPTYPE) \ 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); \ - } \ + PSSORT(vector->n, COMPAREEXPR, SWAPFUNC, SWAPTYPE); \ break; \ } -// END of heap sort code from GSL bool psVectorSortInPlace(const psVector *vector) @@ -422,14 +380,14 @@ switch (vector->type.type) { - 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); + PSVECTOR_SORT_CASE(U8 , PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, U8 ); + PSVECTOR_SORT_CASE(U16, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, U16); + PSVECTOR_SORT_CASE(U32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, U32); + PSVECTOR_SORT_CASE(U64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, U64); + PSVECTOR_SORT_CASE(S8 , PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, S8 ); + PSVECTOR_SORT_CASE(S16, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, S16); + PSVECTOR_SORT_CASE(S32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, S32); + PSVECTOR_SORT_CASE(S64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, S64); + PSVECTOR_SORT_CASE(F32, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, F32); + PSVECTOR_SORT_CASE(F64, PSVECTOR_SORT_COMPARE_DIRECT, PSVECTOR_SORT_SWAP_DIRECT, F64); default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, @@ -484,14 +442,14 @@ } 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); + PSVECTOR_SORT_CASE(U8 , PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(U16, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(U32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(U64, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(S8 , PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(S16, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(S32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(S64, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(F32, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); + PSVECTOR_SORT_CASE(F64, PSVECTOR_SORT_COMPARE_INDEX, PSVECTOR_SORT_SWAP_INDEX, S32); default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, Index: /trunk/psLib/src/types/psArray.c =================================================================== --- /trunk/psLib/src/types/psArray.c (revision 15712) +++ /trunk/psLib/src/types/psArray.c (revision 15713) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.65 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-05 23:42:46 $ + * @version $Revision: 1.66 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-29 21:46:55 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -33,4 +33,5 @@ #include "psAbort.h" #include "psAssert.h" +#include "psSort.h" #define DEFAULT_ARRAY_ADD 10 // Default number to add to an array when not specified @@ -243,4 +244,14 @@ } +// Comparison and swap functions for sorting values directly +#define PSARRAY_SORT_COMPARE_INDEX(A,B) (func(array[index[A]], array[index[B]]) < 0) +#define PSARRAY_SORT_SWAP_INDEX(TYPE,A,B) { \ + if (A != B) { \ + ps##TYPE temp = index[A]; \ + index[A] = index[B]; \ + index[B] = temp; \ + } \ +} + // Heap sort of the index array psVector *psArraySortIndex (psVector *out, psArray *in, psCompareFunc func) { @@ -250,52 +261,9 @@ } - // XXX isn't this line a waste of space? - out = psVectorRecycle(out, in->n, PS_TYPE_S32); // Vector for output out = psVectorCreate(out, 0, in->n, 1, PS_TYPE_S32); - long N = out->n; // Number of elements - if (N < 2) { - return out; - } - - long I = N >> 1; - long J = N - 1; - - int result; - - psS32 temp; - psS32 *index = out->data.S32; - while (1) { - if (I > 0) { - // I--; - temp = index[--I]; - } else { - temp = index[J]; - index[J] = index[0]; - // J--; - if (--J == 0) { - index[0] = temp; - return out; - } - } - long i = I; - long j = (I << 1) + 1; - while (j <= J) { - if (j < J) { - result = func (in->data[index[j]], in->data[index[j+1]]); - if (result < 0) { - j++; - } - } - result = func (in->data[temp], in->data[index[j]]); - if (result < 0) { - index[i]=index[j]; - j += (i=j) + 1; - } else { - j = J + 1; - } - } - index[i] = temp; - } - psAbort ("impossible to reach here"); + psS32 *index = out->data.S32; // Dereference index vector + psPtr *array = in->data; // Dereference input array + PSSORT(out->n, PSARRAY_SORT_COMPARE_INDEX, PSARRAY_SORT_SWAP_INDEX, S32); + return out; }