| [dx,n,iperm,kflag,ier]=dpsort(dx,n,iperm,kflag,ier); |
function [dx,n,iperm,kflag,ier]=dpsort(dx,n,iperm,kflag,ier);
%***BEGIN PROLOGUE DPSORT
%***PURPOSE Return the permutation vector generated by sorting a given
% array and, optionally, rearrange the elements of the array.
% The array may be sorted in increasing or decreasing order.
% A slightly modified quicksort algorithm is used.
%***LIBRARY SLATEC
%***CATEGORY N6A1B, N6A2B
%***TYPE doubleprecision (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H)
%***KEYWORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT
%***AUTHOR Jones, R. E., (SNLA)
% Rhoads, G. S., (NBS)
% Wisniewski, J. A., (SNLA)
%***DESCRIPTION
%
% DPSORT returns the permutation vector IPERM generated by sorting
% the array DX and, optionally, rearranges the values in DX. DX may
% be sorted in increasing or decreasing order. A slightly modified
% quicksort algorithm is used.
%
% IPERM is such that DX(IPERM(I)) is the Ith value in the
% rearrangement of DX. IPERM may be applied to another array by
% calling IPPERM, SPPERM, DPPERM or HPPERM.
%
% The main difference between DPSORT and its active sorting equivalent
% DSORT is that the data are referenced indirectly rather than
% directly. Therefore, DPSORT should require approximately twice as
% long to execute as DSORT. However, DPSORT is more general.
%
% Description of Parameters
% DX - input/output -- doubleprecision array of values to be
% sorted. If ABS(KFLAG) = 2, then the values in DX will be
% rearranged on output; otherwise, they are unchanged.
% N - input -- number of values in array DX to be sorted.
% IPERM - output -- permutation array such that IPERM(I) is the
% index of the value in the original order of the
% DX array that is in the Ith location in the sorted
% order.
% KFLAG - input -- control parameter:
% = 2 means return the permutation vector resulting from
% sorting DX in increasing order and sort DX also.
% = 1 means return the permutation vector resulting from
% sorting DX in increasing order and do not sort DX.
% = -1 means return the permutation vector resulting from
% sorting DX in decreasing order and do not sort DX.
% = -2 means return the permutation vector resulting from
% sorting DX in decreasing order and sort DX also.
% IER - output -- error indicator:
% = 0 if no error,
% = 1 if N is zero or negative,
% = 2 if KFLAG is not 2, 1, -1, or -2.
%***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 by John A. Wisniewski to use the Singleton
% quicksort algorithm.
% 870423 Modified by Gregory S. Rhoads for passive sorting with the
% option for the rearrangement of the original data.
% 890619 doubleprecision version of SPSORT created by D. W. Lozier.
% 890620 Algorithm for rearranging the data vector corrected by R.
% Boisvert.
% 890622 Prologue upgraded to Version 4.0 style by D. Lozier.
% 891128 Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert.
% 920507 Modified by M. McClain to revise prologue text.
% 920818 Declarations section rebuilt and code restructured to use
% IF-THEN-ELSE-ENDIF. (SMR, WRB)
%***end PROLOGUE DPSORT
% .. Scalar Arguments ..
persistent gt100 gt150 gt50 i ij il indx indx0 istrt iu j k kk l lm lmt m nn r temp ;
if isempty(gt50), gt50=0; end;
if isempty(gt100), gt100=0; end;
if isempty(gt150), gt150=0; end;
% .. Array Arguments ..
dx_shape=size(dx);dx=reshape(dx,1,[]);
iperm_shape=size(iperm);iperm=reshape(iperm,1,[]);
% .. Local Scalars ..
if isempty(r), r=0; end;
if isempty(temp), temp=0; end;
if isempty(i), i=0; end;
if isempty(ij), ij=0; end;
if isempty(indx), indx=0; end;
if isempty(indx0), indx0=0; end;
if isempty(istrt), istrt=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(lm), lm=0; end;
if isempty(lmt), lmt=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;
% .. External Subroutines ..
% .. Intrinsic Functions ..
%***FIRST EXECUTABLE STATEMENT DPSORT
ier = 0;
nn = fix(n);
if( nn<1 )
ier = 1;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DPSORT','The number of values to be sorted, N, is not positive.',ier,1);
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
return;
end;
%
kk = fix(abs(kflag));
if( kk~=1 && kk~=2 )
ier = 2;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','DPSORT','The sort control parameter, KFLAG, is not 2, 1, -1, or -2.',ier,1);
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
return;
end;
%
% Initialize permutation vector
%
for i = 1 : nn;
iperm(i) = fix(i);
end; i = fix(nn+1);
%
% Return if only one value is to be sorted
%
if( nn==1 )
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
return;
end;
%
% Alter array DX to get decreasing order if needed
%
if( kflag<=-1 )
for i = 1 : nn;
dx(i) = -dx(i);
end; i = fix(nn+1);
end;
%
% Sort DX only
%
m = 1;
i = 1;
j = fix(nn);
r = .375d0;
%
gt50=0;
gt100=0;
gt150=0;
while( true );
if(gt150==0)
if(gt100==0)
if(gt50==0)
if( i==j )
gt150=1;
continue;
end;
if( r<=0.5898437d0 )
r = r + 3.90625e-2;
else;
r = r - 0.21875d0;
end;
%
end;
gt50=0;
k = fix(i);
%
% Select a central element of the array and savemlv it in location L
%
ij = fix(i + fix((j-i).*r));
lm = fix(iperm(ij));
%
% If first element of array is greater than LM, interchange with LM
%
if( dx(iperm(i))>dx(lm) )
iperm(ij) = fix(iperm(i));
iperm(i) = fix(lm);
lm = fix(iperm(ij));
end;
l = fix(j);
%
% If last element of array is less than LM, interchange with LM
%
if( dx(iperm(j))<dx(lm) )
iperm(ij) = fix(iperm(j));
iperm(j) = fix(lm);
lm = fix(iperm(ij));
%
% If first element of array is greater than LM, interchange
% with LM
%
if( dx(iperm(i))>dx(lm) )
iperm(ij) = fix(iperm(i));
iperm(i) = fix(lm);
lm = fix(iperm(ij));
end;
end;
%
% Find an element in the second half of the array which is smaller
% than LM
%
while( true );
l = fix(l - 1);
if( dx(iperm(l))<=dx(lm) )
%
% Find an element in the first half of the array which is greater
% than LM
%
while( true );
k = fix(k + 1);
if( dx(iperm(k))>=dx(lm) )
break;
end;
end;
%
% Interchange these elements
%
if( k>l )
break;
end;
lmt = fix(iperm(l));
iperm(l) = fix(iperm(k));
iperm(k) = fix(lmt);
end;
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;
%
end;
gt100=0;
if( j-i>=1 )
gt50=1;
continue;
end;
if( i==1 )
continue;
end;
i = fix(i - 1);
%
while( true );
i = fix(i + 1);
if( i==j )
break;
end;
lm = fix(iperm(i+1));
if( dx(iperm(i))>dx(lm) )
k = fix(i);
%
while( true );
iperm(k+1) = fix(iperm(k));
k = fix(k - 1);
%
if( dx(lm)>=dx(iperm(k)) )
break;
end;
end;
iperm(k+1) = fix(lm);
end;
end;
end;
gt150=0;
%
% Begin again on another portion of the unsorted array
%
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
gt100=1;
end;
%
% Clean up
%
if( kflag<=-1 )
for i = 1 : nn;
dx(i) = -dx(i);
end; i = fix(nn+1);
end;
%
% Rearrange the values of X if desired
%
if( kk==2 )
%
% use the IPERM vector as a flag.
% If IPERM(I) < 0, then the I-th value is in correct location
%
for istrt = 1 : nn;
if( iperm(istrt)>=0 )
indx = fix(istrt);
indx0 = fix(indx);
temp = dx(istrt);
while( iperm(indx)>0 );
dx(indx) = dx(iperm(indx));
indx0 = fix(indx);
iperm(indx) = fix(-iperm(indx));
indx = fix(abs(iperm(indx)));
end;
dx(indx0) = temp;
end;
end; istrt = fix(nn+1);
%
% Revert the signs of the IPERM values
%
for i = 1 : nn;
%
iperm(i) = fix(-iperm(i));
end; i = fix(nn+1);
end;
%
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
end %subroutine dpsort
%!!
%!!100 IF ( i==j ) GOTO 500
%!!IF ( r<=0.5898437D0 ) THEN
%!! r = r + 3.90625D-2
%!!ELSE
%!! r = r - 0.21875D0
%!!ENDIF
%!!!
%!!200 k = i
%!!!
%!!! Select a central element of the array and savemlv it in location L
%!!!
%!!ij = i + INT((j-i)*r)
%!!lm = Iperm(ij)
%!!!
%!!! If first element of array is greater than LM, interchange with LM
%!!!
%!!IF ( Dx(Iperm(i))>Dx(lm) ) THEN
%!! Iperm(ij) = Iperm(i)
%!! Iperm(i) = lm
%!! lm = Iperm(ij)
%!!ENDIF
%!!l = j
%!!!
%!!! If last element of array is less than LM, interchange with LM
%!!!
%!!IF ( Dx(Iperm(j))<Dx(lm) ) THEN
%!! Iperm(ij) = Iperm(j)
%!! Iperm(j) = lm
%!! lm = Iperm(ij)
%!! !
%!! ! If first element of array is greater than LM, interchange
%!! ! with LM
%!! !
%!! IF ( Dx(Iperm(i))>Dx(lm) ) THEN
%!! Iperm(ij) = Iperm(i)
%!! Iperm(i) = lm
%!! lm = Iperm(ij)
%!! ENDIF
%!!ENDIF
%!!!
%!!! Find an element in the second half of the array which is smaller
%!!! than LM
%!!!
%!!300 l = l - 1
%!!IF ( Dx(Iperm(l))>Dx(lm) ) GOTO 300
%!!!
%!!! Find an element in the first half of the array which is greater
%!!! than LM
%!!!
%!!400 k = k + 1
%!!IF ( Dx(Iperm(k))<Dx(lm) ) GOTO 400
%!!!
%!!! Interchange these elements
%!!!
%!!IF ( k<=l ) THEN
%!! lmt = Iperm(l)
%!! Iperm(l) = Iperm(k)
%!! Iperm(k) = lmt
%!! GOTO 300
%!!ELSE
%!! !
%!! ! Save upper and lower subscripts of the array yet to be sorted
%!! !
%!! IF ( l-i>j-k ) THEN
%!! il(m) = i
%!! iu(m) = l
%!! i = k
%!! m = m + 1
%!! ELSE
%!! il(m) = k
%!! iu(m) = j
%!! j = l
%!! m = m + 1
%!! ENDIF
%!! GOTO 600
%!!ENDIF
%!!!
%!!! Begin again on another portion of the unsorted array
%!!!
%!!500 m = m - 1
%!!IF ( m==0 ) THEN
%!! !
%!! ! Clean up
%!! !
%!! IF ( Kflag<=-1 ) THEN
%!! DO i = 1 , nn
%!! Dx(i) = -Dx(i)
%!! ENDDO
%!! ENDIF
%!! !
%!! ! Rearrange the values of DX if desired
%!! !
%!! IF ( kk==2 ) THEN
%!! !
%!! ! use the IPERM vector as a flag.
%!! ! If IPERM(I) < 0, then the I-th value is in correct location
%!! !
%!! DO istrt = 1 , nn
%!! IF ( Iperm(istrt)>=0 ) THEN
%!! indx = istrt
%!! indx0 = indx
%!! temp = Dx(istrt)
%!!505 IF ( Iperm(indx)>0 ) THEN
%!! Dx(indx) = Dx(Iperm(indx))
%!! indx0 = indx
%!! Iperm(indx) = -Iperm(indx)
%!! indx = ABS(Iperm(indx))
%!! GOTO 505
%!! ENDIF
%!! Dx(indx0) = temp
%!! ENDIF
%!! ENDDO
%!! !
%!! ! Revert the signs of the IPERM values
%!! !
%!! DO i = 1 , nn
%!! Iperm(i) = -Iperm(i)
%!! ENDDO
%!! !
%!! ENDIF
%!! RETURN
%!!ELSE
%!! i = il(m)
%!! j = iu(m)
%!!ENDIF
%!!!
%!!600 IF ( j-i>=1 ) GOTO 200
%!!IF ( i==1 ) GOTO 100
%!!i = i - 1
%!!!
%!!700 i = i + 1
%!!IF ( i==j ) GOTO 500
%!!lm = Iperm(i+1)
%!!IF ( Dx(Iperm(i))<=Dx(lm) ) GOTO 700
%!!k = i
%!!!
%!!800 Iperm(k+1) = Iperm(k)
%!!k = k - 1
%!!IF ( Dx(lm)<Dx(Iperm(k)) ) GOTO 800
%!!Iperm(k+1) = lm
%!!GOTO 700
%!!!
%!!99999 end subroutine DPSORT
%DECK DPTSL
|
|
