By Bill McKeeman, MathWorks and Loren Shure, MathWorks

While choosing how to upgrade sort for MATLAB 7, we used MATLAB to prototype `sort`

algorithms and to add the newsort direction qualifiers `'ascend'`

and `'descend'`

. We translated the prototype from MATLAB into C++ to add this latest version of sort to Release 14—and found a few surprises in the process.

Algorithms in MathWorks products and user applications use the MATLAB `sort`

There are 12 types of sortable matrices in MATLAB: the eight integer types `([u]int8/16/32/64);`

plus single, double, char, and logical. Single and double may be real or complex; double and logical can be sparse. Full n-dimensional (ND) arrays can be sorted along any selected dimension, in either ascending or descending order. `sort`

can also return the permutation vector relating the input values to the outputs, with call syntax ranging from the most simple

r = sort(v);

to the most complete

[r,i] = sort(v, dim, dirflag);

The most readily available upgrade would use qsort from the C library, but qsort was too specific for the sort tasks that MATLAB requires. The C library `qsort` has a hard-coded stride of 1 (data is assumed to be arranged in contiguous memory). MATLAB also needs to sort data containing NaN values, which `qsort` does not accommodate.

The available C `qsort` implementations share a common characteristic: to achieve type independence, data access is indirect (via void* pointers) and the user supplies a comparison function. As a result, general-purpose sort algorithms written in C are inelegant and somewhat slow.

Sorting experts recommend a three-level algorithm,
*insert/quick/bin* , where each level calls the one to its left if the size of the array being sorted passes below some experimentally discovered cutoff. When we updated
*insert/quick/bin/quick* ). To achieve speed proportional to the size of the data (order
*N* )
*bin sort* requires as much auxiliary storage for the bins as it does for the input data. Because
*quicksort* manipulates the data in place, it can sort larger arrays, albeit somewhat more slowly. (Hoare’s 1962 quick sort algorithm has undergone a generation of changes, particularly to improve performance when the data is not random.)

*Resources*

Bentley, J. L. and M. D. McIlroy, 1993. Engineering a sort function. Software—Practice and Experience 23, 11, 1249--1265.

Knuth, D. E. 1998. The Art of Computer Programming, Volume 3: Sorting and Searching, 2nd ed. Boston: Addison-Wesley.

Cormen et al. 2001. Introduction to Algorithms, 2nd ed. Columbus: McGraw-Hill.

`sort`

, we expected to require four levels (
Bentley, J. L. and M. D. McIlroy, 1993. Engineering a sort function. Software—Practice and Experience 23, 11, 1249--1265.

Knuth, D. E. 1998. The Art of Computer Programming, Volume 3: Sorting and Searching, 2nd ed. Boston: Addison-Wesley.

Cormen et al. 2001. Introduction to Algorithms, 2nd ed. Columbus: McGraw-Hill.

We explored speed and memory tradeoffs by prototyping the basic algorithms in M. The prototypes work only on vectors. (Full ND arrays are discussed in Extending Vector Sorts to N-Dimensions later in this article.)

```
function v = viSort(v) % vector insert
for i=1:numel(v)-1
[m, k] = min(v(i:end));
v(i+1:i+k-1) = v(i:i+k-2) % make room
v(i) = m; % place min
end
```

```
function v = vqSort(v) % vector quick sort
if numel(v) > 1
piv = median(v); % divide
lss = v(v < piv); % and
eql = v(v == piv); % conquer
gtr = v(v > piv);
v = [vqSort(lss), eql, vqSort(gtr)];
end
```

```
function v = vbSort(v) % vector bin sort
n = numel(v);
if n > 1
mn = min(v);
mx = max(v);
binAvg = 200; % experimental
binCt = ceil(n/binAvg);
slope = binCt/(mx-mn);
binNo = floor(slope.*(v-mn))+1; % assign bin
bins = cell(1, binCt+1); % allocate bins
for i=1:binCt+1
bins{i} = vqSort(v(binNo==i));
end
v = [bins{:}]; % reassemble
end
```

As we planned to have `vqSort`

call `viSort`

when there were few enough data points, we modified the *quick sort.* design.

```
function v = vqSort(v) % quick sort
if numel(v) < 32
v = viSort(v); % experimental
else
piv = median(v); % need approximation
lss = v(v < piv);
eql = v(v == piv);
gtr = v(v > piv);
v = [vqSort(lss), eql, vqSort(gtr)];
end
```

We translated the M-code vector sorts into C in-place sorts so that we could test our four-level design hypothesis. The results surprised us. For vectors of random data on Intel x86 hardware, the most modern *quick sorts* were slower than the already-implemented *quick sort* in MATLAB! Algorithmic changes made to the modern sorts to ensure good behavior on less-than-random input were costly, while the existing MATLAB sort cleverly uses its pivot-picking approximation to do some sorting work.

We finally decided that best performance on random data was more useful than adequate performance on all data, and kept the existing MATLAB *quick sort*. (Performance results depend somewhat on the underlying hardware.)

We took into account that the C implementation of *bin sort* needs three passes: one pass to decide on the number of bins, one to bin the data, and one to call qsort on each bin. The results are left in place, conveniently ordered for output. We can justify this considerable overhead because it results in an order *N*, instead of *N*log*N*, run-time algorithm.

Another surprise! We experimented to find a cutover from *quick sort* to *bin sort*, but never found it. Our best effort at a C-implemented *bin sort* was slower than our C-implemented *quick sort*, even on datasets containing up to 10 million elements. As a result, we abbreviated our development plan to a two-stage sort: *insert/quick*. Both these algorithms already have satisfactory implementations in MATLAB.

`NaN`

and `Inf`

Next, we extended the double sorts to the sortable MATLAB types. We recoded the `insert/quick`

sort pair as C++ template functions parameterized by the underlying MATLAB types.

MATLAB sorts `NaN`

values high, beyond all other single/double values. Because `NaN`

and `Inf`

need to be dealt with only for the real types, we implemented a function (called a jacket) to handle the cases for single and double. The jacket calls the usual sort routine on the finite values and handles the nonfinite values correctly for MATLAB. This jacket’s sorting prepass counts the `NaN`

s and replaces them with `+Inf`

. After the sort, the jacket overwrites an appropriate number of `+Inf`

values with `NaN`

s at the top end of the sorted vector. (Had we still been doing *bin sort*, this jacket would also have had to count and remove `Inf`

values, as the bin selection algorithm depends on finite values.)

```
function v = vSort(v) % sort NaN jacket
if isfloat(v)
nans = isnan(v);
v(nans) = inf;
v = vqSort(v);
v(end-sum(nans)+1:end) = nan;
else
v = vqSort(v);
end
```

We extended the vector sorts to sort along a user-selected dimension of an ND array. Array data is laid out in memory in column-major order: elements in a column are adjacent in memory. Elements along another dimension are separated by a constant distance (the so-called stride, usually larger than 1). We prototyped in M first

```
function v = strideSort(v, d) % stride jacket
s = size(v);
stride = prod(s(1:d-1));
for i = 1 : stride*s(d) : numel(v) % #of elements
for j = i : i+stride-1
idx = j : stride : j+stride*s(d)-1;
v(idx) = vSort(v(idx));
end
end
```

We recoded the C++ in-place *insert/quick* sort pair to step by a stride, and then wrote a second level of C++ jacket to call them.

We extended the array sorts to provide the permutation vectors, which required that the sort be stable, (i.e., sorting may not reorder equal values). The algorithms `viSort`

, `vqSort`

, and `vbSort`

are stable, but our C++ implementations were not. Without requesting the permutation vector, you cannot tell whether the sort is unstable.

Rather than recode the *insert/quick* sorts, we chose to make a post-sort pass over the sorted data, detect runs of equal values, and then sort the corresponding subvector in the parallel index vector. Any `NaN`

s presented a special problem because they had been interspersed with the `Inf`

values. We found the original indices of the `NaN`

s and `Inf`

s by making a pass over the original unsorted data.

The MATLAB algorithm for complex values first sorts by absolute value, then sorts by phase angle, and finally sorts to ensure stability for duplicate values. We decided to leave the MATLAB C code for sorting complex double values in place and call it from the jackets described above.

We passed the direction flags '`ascend`

' and '`descend`

' to the jackets. If the flag is '`descend`

', the output is reversed just before being returned. After the reversal, we ran the index vector sorts to restore stability.

We tested for the following criteria:

- The sorts should accept all valid inputs (12 types of data, shapes from empty arrays to large ND ensembles, and valid dimension and direction flags) and otherwise fail with an error message.
- If requested, the sorts should produce the permutation vector.
- The sorts should “run as fast as C.” Suppose T means “type” and S means “shape.” Then we guarantee the following.

Statement | Guarantees |
---|---|

r = sort(v, dim, flags) | T(r)=T(v) _{^} S(r)=S(v) |

[r,i] = sort(v, dim, flags) | T(r)=T(v) _{^} T(i)=double _{^} S(r)=S(i)=S(v) |

We checked that in each case the outputs were sorted and the permutation vectors applied to the inputs, given the outputs. We coded all the test routines in M, using the MATLAB functions `class`

and `size`

.

To test performance, we raced our algorithm against the C standard library `qsort`

and the C++ Standard Template Library `sort`

for simple cases. The MATLAB `sort`

function is about five times faster than C `qsort`

and a little faster than the C++ STL `sort`

.

MATLAB provided a productive environment to explore enhancements to the `sort`

algorithm. Prototyping the possibilities in MATLAB before recoding the algorithms in C++ helped us quickly understand our programming choices. We can still use the prototype M versions of `sort`

as test components when we verify the MATLAB `sort`

function’s behavior. Using MATLAB to upgrade MATLAB enabled us to deliver a higher performance `sort`

, and we can use M-code to explain how we did it.

`Similar to the MATLAB function sort with options 'ascend' and 'descend’, the function find has added functionality in MATLAB 7 that lets you find the 'first' or 'last' elements that meet some criterion. For example, if you want to find the first element in array X less than 42, use`

find`find(X<42,1,'first')`

This operation stops as soon as the first element meeting the condition is encountered, which can result in significant time savings if the array is large and the value you find happens to be far from the end. If, instead, you want to find the final 10 elements > 17, then use`find(X>17,10,'last')`

Again, the operation returns as soon as 10 values are found. If fewer than the requested number is found after the entire array is scanned, the return vector will have fewer than the requested number of elements. Before`sort`

had the option to return items in descending order, you had to perform the`sort`

and then index from the end. With MATLAB 7,`sort`

has the new '`descend'`

option, and`find`

has the new`'last'`

option to combine the find with the indexing into a single function call.

Published 2004