//# GenSort.h: General sort functions //# Copyright (C) 1993,1994,1995,1996,1997,1999 //# Associated Universities, Inc. Washington DC, USA. //# //# This library is free software; you can redistribute it and/or modify it //# under the terms of the GNU Library General Public License as published by //# the Free Software Foundation; either version 2 of the License, or (at your //# option) any later version. //# //# This library 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 Library General Public //# License for more details. //# //# You should have received a copy of the GNU Library General Public License //# along with this library; if not, write to the Free Software Foundation, //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. //# //# Correspondence concerning AIPS++ should be addressed as follows: //# Internet email: aips2-request@nrao.edu. //# Postal address: AIPS++ Project Office //# National Radio Astronomy Observatory //# 520 Edgemont Road //# Charlottesville, VA 22903-2475 USA //# //# $Id$ #ifndef CASA_GENSORT_H #define CASA_GENSORT_H #include #include #include namespace casacore { //# NAMESPACE CASACORE - BEGIN //# Forward declarations. template class Block; // General in-place sort functions // // // // // The static member functions of this templated class are highly optimized // sort functions. They do an in-place sort of an array of values. The // functions are templated, so they can in principle be used with any // data type. However, if used with non-builtin data types, their // class must provide certain functions (see Template Type Argument // Requirements). // // If it is impossible or too expensive to define these functions, the // Sort class can be used instead. This sorts // indirectly using an index array. Instead of the functions mentioned // above it requires a comparison routine. // // The GenSort functions can sort: //
    //
  • C-arrays of values; //
  • Arrays of values -- the array can have any shape // and the increment can be >1; //
  • Blocks of values -- there is a special function to // sort less elements than the size of the Block. //
// // The sort order can be specified in the order field: //
//
Sort::Ascending (default), //
Sort::Descending. //
// // Previously the sort algorithm to use could be given in the options field. //
//
Sort::QuickSort (default) //
is the fastest. It is about 4-6 times faster // than the qsort function on the SUN. No worst case has been // found, even not for cases where qsort is terribly slow. //
Sort::HeapSort //
is about twice as slow as quicksort. // It has the advantage that the worst case is always o(n*log(n)), // while quicksort can have hypothetical inputs with o(n*n). //
Sort::InsSort //
is o(n*n) for random inputs. It is, however, the // only stable sort (i.e. equal values remain in the original order). //
// However, these options are not used anymore because the sort now always // uses a merge sort that is equally fast for random input and much faster for // degenerated cases like an already ordered or reversely ordered array. // Furthermore, merge sort is always stable and will be parallelized if OpenMP // support is enabled giving a 6-fold speedup on 8 cores. //
Sort::NoDuplicates in the options field indicates that // duplicate values will be removed (only the first occurrance is kept). //
The previous sort functionality is still available through the functions // quickSort, heapSort, and insSort. //

// All the sort functions return the number of values sorted as their // function value; when duplicate values have been removed, the number of // unique valuess will be returned. //

// The class also provides a function to find the k-th largest value // in an array of values. This uses a stripped-down version of quicksort // and is at least 6 times faster than a full quicksort. // // //

  • operator= to assign when swapping elements //
  • operator<, operator> and // operator== to compare elements //
  • default constructor to allocate a temporary // template class GenSort { public: // Sort a C-array containing nr T-type objects. // The sort is done in place and is always stable (thus equal keys keep // their original order). It returns the number of values, which // can be different if a NoDuplicates sort is done. //
    Insertion sort is used for short arrays (<50 elements). Otherwise, // a merge sort is used which will be parallelized if casacore is built // with OpenMP support. // static uInt sort (T*, uInt nr, Sort::Order = Sort::Ascending, int options = 0); static uInt sort (Array&, Sort::Order = Sort::Ascending, int options = 0); static uInt sort (Block&, uInt nr, Sort::Order = Sort::Ascending, int options = 0); // // Find the k-th largest value. //
    Note: it does a partial quicksort, thus the data array gets changed. static T kthLargest (T* data, uInt nr, uInt k); // Sort C-array using quicksort. static uInt quickSort (T*, uInt nr, Sort::Order = Sort::Ascending, int options = 0); // Sort C-array using heapsort. static uInt heapSort (T*, uInt nr, Sort::Order = Sort::Ascending, int options = 0); // Sort C-array using insertion sort. static uInt insSort (T*, uInt nr, Sort::Order = Sort::Ascending, int options = 0); // Sort C-array using parallel merge sort (using OpenMP). // By default OpenMP determines the number of threads that can be used. static uInt parSort (T*, uInt nr, Sort::Order = Sort::Ascending, int options = 0, int nthread = 0); // Swap 2 elements in array. static inline void swap (T&, T&); // Reverse the elements in res and put them into data. // Care is taken if both pointers reference the same data. static void reverse (T* data, const T* res, uInt nrrec); private: // Thedata buffer is divided in nparts parts. // In each part the values are in ascending order. // The index tells the nr of elements in each part. // Recursively each two subsequent parts are merged until only part is left // (giving the sorted array). Alternately data and tmp // are used for the merge result. The pointer containing the final result // is returned. //
    If possible, merging the parts is done in parallel (using OpenMP). static T* merge (T* data, T* tmp, uInt nrrec, uInt* index, uInt nparts); // Quicksort in ascending order. static void quickSortAsc (T*, Int, Bool multiThread=False, Int rec_lim=128); // Heapsort in ascending order. static void heapSortAsc (T*, Int); // Helper function for ascending heapsort. static void heapAscSiftDown (Int, Int, T*); // Insertion sort in ascending order. static uInt insSortAsc (T*, Int, int option); // Insertion sort in ascending order allowing duplicates. // This is also used by quicksort for its last steps. static uInt insSortAscDup (T*, Int); // Insertion sort in ascending order allowing no duplicates. // This is also used by the other sort algorithms to skip duplicates. static uInt insSortAscNoDup (T*, Int); }; // General indirect sort functions // // // // This class is similar to GenSort. // The only difference is that the functions in this class sort // indirectly instead of in-place. // They return the result of the sort as an sorted vector of indices // This is slower, because an extra indirection is involved in each // comparison. However, this sort allows to sort const data. // Another advantage is that this sort is always stable (i.e. equal // values are kept in their original order). template class GenSortIndirect { public: // Sort a C-array containing nr T-type objects. // The resulting index vector gives the sorted indices. static INX sort (Vector& indexVector, const T* data, INX nr, Sort::Order = Sort::Ascending, int options = Sort::QuickSort); // Sort a C-array containing nr T-type objects. // The resulting index vector gives the sorted indices. static INX sort (Vector& indexVector, const Array& data, Sort::Order = Sort::Ascending, int options = Sort::QuickSort); // Sort a C-array containing nr T-type objects. // The resulting index vector gives the sorted indices. static INX sort (Vector& indexVector, const Block& data, INX nr, Sort::Order = Sort::Ascending, int options = Sort::QuickSort); // Find the index of the k-th largest value. static INX kthLargest (T* data, INX nr, INX k); // Sort container using quicksort. // The argument inx gives the index defining the order of the // values in the data array. Its length must be at least nr // and it must be filled with the index values of the data. // Usually this is 0..nr, but it could contain a selection of the data. static INX quickSort (INX* inx, const T* data, INX nr, Sort::Order, int options); // Sort container using heapsort. static INX heapSort (INX* inx, const T* data, INX nr, Sort::Order, int options); // Sort container using insertion sort. static INX insSort (INX* inx, const T* data, INX nr, Sort::Order, int options); // Sort container using parallel merge sort (using OpenMP). // By default the maximum number of threads is used. static INX parSort (INX* inx, const T* data, INX nr, Sort::Order, int options, int nthreads=0); private: // Swap 2 indices. static inline void swapInx (INX& index1, INX& index2); // Thedata buffer is divided in nparts parts. // In each part the values are in ascending order. // The index tells the nr of elements in each part. // Recursively each two subsequent parts are merged until only part is left // (giving the sorted array). Alternately data and tmp // are used for the merge result. The pointer containing the final result // is returned. //
    If possible, merging the parts is done in parallel (using OpenMP). static INX* merge (const T* data, INX* inx, INX* tmp, INX nrrec, INX* index, INX nparts); // Check if 2 values are in ascending order. // When equal, the order is correct if index1 Global in-place sort functions // The following global functions are easier to use than the static // GenSort member functions. // They do an in-place sort of data, thus the data themselves are moved // ending up in the requested order. //

    // The default sorting method is QuickSort, which is the fastest. // However, there are pathological cases where it can be slow. // HeapSort is about twice a slow, but its speed is guaranteed. // InsSort (insertion sort) can be very, very slow, but it is the only // stable sort method (i.e. equal values are kept in their original order). // However, indirect sorting methods // are available to make QuickSort and HeapSort stable. //

    // All sort methods have an option to skip duplicate values. This is the // only case where the returned number of values can be less than the // original number of values. // template inline uInt genSort (T* data, uInt nr, Sort::Order order = Sort::Ascending, int options=0) { return GenSort::sort (data, nr, order, options); } template inline uInt genSort (Array& data, Sort::Order order = Sort::Ascending, int options=0) { return GenSort::sort (data, order, options); } template inline uInt genSort (Block& data, Sort::Order order = Sort::Ascending, int options=0) { return GenSort::sort (data, data.nelements(), order, options); } template inline uInt genSort (Block& data, uInt nr, Sort::Order order = Sort::Ascending, int options=0) { return GenSort::sort (data, nr, order, options); } // //

    Global indirect sort functions // The following global functions easier to use than the static // GenSortIndirect member functions. // They do an indirect sort of data, thus the data themselves are not moved. // Rather an index vector is returned giving the sorted data indices. //

    // The sorting method used is merge sort, which is always stable. // It is the fastest, especially if it can use multiple threads. //

    // Unlike the in-place sorting methods // , all indirect sorting methods are stable (i.e. equal // values are left in their original order). //

    // All sort methods have an option to skip duplicate values. This is the // only case where the returned number of values can be less than the // original number of values. // template inline uInt genSort (Vector& indexVector, const T* data, INX nr, Sort::Order order = Sort::Ascending, int options=0) { return GenSortIndirect::sort (indexVector, data, nr, order, options); } template inline uInt genSort (Vector& indexVector, const Array& data, Sort::Order order = Sort::Ascending, int options=0) { return GenSortIndirect::sort (indexVector, data, order, options); } template inline uInt genSort (Vector& indexVector, const Block& data, Sort::Order order = Sort::Ascending, int options=0) { return GenSortIndirect::sort (indexVector, data, data.nelements(), order, options); } template inline uInt genSort (Vector& indexVector, const Block& data, INX nr, Sort::Order order = Sort::Ascending, int options=0) { return GenSortIndirect::sort (indexVector, data, nr, order, options); } // // Implement inline member functions. template inline void GenSort::swap (T& l, T& r) { T t = l; l = r; r = t; } template inline void GenSortIndirect::swapInx (INX& i, INX& j) { INX t = i; i = j; j = t; } template inline int GenSortIndirect::isAscending (const T* data, INX i, INX j) { return (data[i] > data[j] || (data[i] == data[j] && i > j)); } } //# NAMESPACE CASACORE - END #ifndef CASACORE_NO_AUTO_TEMPLATES #include #endif //# CASACORE_NO_AUTO_TEMPLATES #endif