| [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
|
|
