Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[dx,dy,n,kflag]=dsort(dx,dy,n,kflag);
```function [dx,dy,n,kflag]=dsort(dx,dy,n,kflag);
%***BEGIN PROLOGUE  DSORT
%***PURPOSE  Sort an array and optionally make the same interchanges in
%            an auxiliary array.  The array may be sorted in increasing
%            or decreasing order.  A slightly modified QUICKSORT
%            algorithm is used.
%***LIBRARY   SLATEC
%***CATEGORY  N6A2B
%***TYPE      doubleprecision (SSORT-S, DSORT-D, ISORT-I)
%***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
%***AUTHOR  Jones, R. E., (SNLA)
%           Wisniewski, J. A., (SNLA)
%***DESCRIPTION
%
%   DSORT sorts array DX and optionally makes the same interchanges in
%   array DY.  The array DX may be sorted in increasing order or
%   decreasing order.  A slightly modified quicksort algorithm is used.
%
%   Description of Parameters
%      DX - array of values to be sorted   (usually abscissas)
%      DY - array to be (optionally) carried along
%      N  - number of values in array DX to be sorted
%      KFLAG - control parameter
%            =  2  means sort DX in increasing order and carry DY along.
%            =  1  means sort DX in increasing order (ignoring DY)
%            = -1  means sort DX in decreasing order (ignoring DY)
%            = -2  means sort DX in decreasing order and carry DY along.
%
%***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
%                 for sorting with minimal storage, Communications of
%                 the ACM, 12, 3 (1969), pp. 185-187.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   761101  DATE WRITTEN
%   761118  Modified to use the Singleton quicksort algorithm.  (JAW)
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891009  Removed unreferenced statement labels.  (WRB)
%   891024  Changed category.  (WRB)
%   891024  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   901012  Declared all variables; changed X,Y to DX,DY; changed
%           code to parallel SSORT. (M. McClain)
%   920501  Reformatted the REFERENCES section.  (DWL, WRB)
%   920519  Clarified error messages.  (DWL)
%   920801  Declarations section rebuilt and code restructured to use
%           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
%***end PROLOGUE  DSORT
%     .. Scalar Arguments ..
%     .. Array Arguments ..
persistent gt i ij il iu j k kk l m nn r t tt tty ty ;

dx_shape=size(dx);dx=reshape(dx,1,[]);
dy_shape=size(dy);dy=reshape(dy,1,[]);
%     .. Local Scalars ..
if isempty(r), r=0; end;
if isempty(t), t=0; end;
if isempty(tt), tt=0; end;
if isempty(tty), tty=0; end;
if isempty(ty), ty=0; end;
if isempty(i), i=0; end;
if isempty(ij), ij=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
if isempty(nn), nn=0; end;
%     .. Local Arrays ..
if isempty(il), il=zeros(1,21); end;
if isempty(iu), iu=zeros(1,21); end;
if isempty(gt), gt=zeros(1,16); end;
%     .. External Subroutines ..
%     .. Intrinsic Functions ..
%***FIRST EXECUTABLE STATEMENT  DSORT
nn = fix(n);
if( nn<1 )
xermsg('SLATEC','DSORT','The number of values to be sorted is not positive.',1,1);
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
dy_shape=zeros(dy_shape);dy_shape(:)=dy(1:numel(dy_shape));dy=dy_shape;
return;
end;
%
kk = fix(abs(kflag));
if( kk~=1 && kk~=2 )
xermsg('SLATEC','DSORT','The sort control parameter, K, is not 2, 1, -1, or -2.',2,1);
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
dy_shape=zeros(dy_shape);dy_shape(:)=dy(1:numel(dy_shape));dy=dy_shape;
return;
end;
%
%     Alter array DX to get decreasing order if needed
%
gt(:)=0;
while (1);
if(gt(16)==0)
if(gt(15)==0)
if(gt(14)==0)
if(gt(13)==0)
if(gt(12)==0)
if(gt(11)==0)
if(gt(10)==0)
if(gt(9)==0)
if(gt(8)==0)
if(gt(7)==0)
if(gt(6)==0)
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( kflag<=-1 )
for i = 1 : nn;
dx(i) = -dx(i);
end; i = fix(nn+1);
end;
%
if( kk==2 )
%
%     Sort DX and carry DY along
%
m = 1;
i = 1;
j = fix(nn);
r = 0.375d0;
gt(9)=1;
continue;
else;
%
%     Sort DX only
%
m = 1;
i = 1;
j = fix(nn);
r = 0.375d0;
end;
%
end;
gt(1)=0;
if( i==j )
gt(5)=1;
continue;
end;
if( r<=0.5898437d0 )
r = r + 3.90625d-2;
else;
r = r - 0.21875d0;
end;
%
end;
gt(2)=0;
k = fix(i);
%
%     Select a central element of the array and savemlv it in location T
%
ij = fix(i + fix((j-i).*r));
t = dx(ij);
%
%     If first element of array is greater than T, interchange with T
%
if( dx(i)>t )
dx(ij) = dx(i);
dx(i) = t;
t = dx(ij);
end;
l = fix(j);
%
%     If last element of array is less than than T, interchange with T
%
if( dx(j)<t )
dx(ij) = dx(j);
dx(j) = t;
t = dx(ij);
%
%        If first element of array is greater than T, interchange with T
%
if( dx(i)>t )
dx(ij) = dx(i);
dx(i) = t;
t = dx(ij);
end;
end;
%
%     Find an element in the second half of the array which is smaller
%     than T
%
end;
gt(3)=0;
l = fix(l - 1);
if( dx(l)>t )
gt(3)=1;
continue;
end;
%
%     Find an element in the first half of the array which is greater
%     than T
%
end;
gt(4)=0;
k = fix(k + 1);
if( dx(k)<t )
gt(4)=1;
continue;
end;
%
%     Interchange these elements
%
if( k<=l )
tt = dx(l);
dx(l) = dx(k);
dx(k) = tt;
gt(3)=1;
continue;
end;
%
%     Save upper and lower subscripts of the array yet to be sorted
%
if( l-i>j-k )
il(m) = fix(i);
iu(m) = fix(l);
i = fix(k);
m = fix(m + 1);
else;
il(m) = fix(k);
iu(m) = fix(j);
j = fix(l);
m = fix(m + 1);
end;
gt(6)=1;
continue;
%
%     Begin again on another portion of the unsorted array
%
end;
gt(5)=0;
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
%
end;
gt(6)=0;
if( j-i>=1 )
gt(2)=1;
continue;
end;
if( i==1 )
gt(1)=1;
continue;
end;
i = fix(i - 1);
%
end;
gt(7)=0;
i = fix(i + 1);
if( i==j )
gt(5)=1;
continue;
end;
t = dx(i+1);
if( dx(i)<=t )
gt(7)=1;
continue;
end;
k = fix(i);
%
end;
gt(8)=0;
dx(k+1) = dx(k);
k = fix(k - 1);
if( t<dx(k) )
gt(8)=1;
continue;
end;
dx(k+1) = t;
gt(7)=1;
continue;
%
end;
gt(9)=0;
if( i==j )
gt(13)=1;
continue;
end;
if( r<=0.5898437d0 )
r = r + 3.90625d-2;
else;
r = r - 0.21875d0;
end;
%
end;
gt(10)=0;
k = fix(i);
%
%     Select a central element of the array and savemlv it in location T
%
ij = fix(i + fix((j-i).*r));
t = dx(ij);
ty = dy(ij);
%
%     If first element of array is greater than T, interchange with T
%
if( dx(i)>t )
dx(ij) = dx(i);
dx(i) = t;
t = dx(ij);
dy(ij) = dy(i);
dy(i) = ty;
ty = dy(ij);
end;
l = fix(j);
%
%     If last element of array is less than T, interchange with T
%
if( dx(j)<t )
dx(ij) = dx(j);
dx(j) = t;
t = dx(ij);
dy(ij) = dy(j);
dy(j) = ty;
ty = dy(ij);
%
%        If first element of array is greater than T, interchange with T
%
if( dx(i)>t )
dx(ij) = dx(i);
dx(i) = t;
t = dx(ij);
dy(ij) = dy(i);
dy(i) = ty;
ty = dy(ij);
end;
end;
%
%     Find an element in the second half of the array which is smaller
%     than T
%
end;
gt(11)=0;
l = fix(l - 1);
if( dx(l)>t )
gt(11)=1;
continue;
end;
%
%     Find an element in the first half of the array which is greater
%     than T
%
end;
gt(12)=0;
k = fix(k + 1);
if( dx(k)<t )
gt(12)=1;
continue;
end;
%
%     Interchange these elements
%
if( k<=l )
tt = dx(l);
dx(l) = dx(k);
dx(k) = tt;
tty = dy(l);
dy(l) = dy(k);
dy(k) = tty;
gt(11)=1;
continue;
end;
%
%     Save upper and lower subscripts of the array yet to be sorted
%
if( l-i>j-k )
il(m) = fix(i);
iu(m) = fix(l);
i = fix(k);
m = fix(m + 1);
else;
il(m) = fix(k);
iu(m) = fix(j);
j = fix(l);
m = fix(m + 1);
end;
gt(14)=1;
continue;
%
%     Begin again on another portion of the unsorted array
%
end;
gt(13)=0;
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
%
end;
gt(14)=0;
if( j-i>=1 )
gt(10)=1;
continue;
end;
if( i==1 )
gt(9)=1;
continue;
end;
i = fix(i - 1);
%
end;
gt(15)=0;
i = fix(i + 1);
if( i==j )
gt(13)=1;
continue;
end;
t = dx(i+1);
ty = dy(i+1);
if( dx(i)<=t )
gt(15)=1;
continue;
end;
k = fix(i);
%
end;
gt(16)=0;
dx(k+1) = dx(k);
dy(k+1) = dy(k);
k = fix(k - 1);
if( t<dx(k) )
gt(16)=1;
continue;
end;
dx(k+1) = t;
dy(k+1) = ty;
gt(15)=1;
continue;
%
%     Clean up
%
break;
end;
if( kflag<=-1 )
for i = 1 : nn;
dx(i) = -dx(i);
end; i = fix(nn+1);
end;
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
dy_shape=zeros(dy_shape);dy_shape(:)=dy(1:numel(dy_shape));dy=dy_shape;
end %subroutine dsort
%DECK DSOSEQ
```