| [ix,n,iperm,kflag,ier]=ipsort(ix,n,iperm,kflag,ier); |
function [ix,n,iperm,kflag,ier]=ipsort(ix,n,iperm,kflag,ier);
%***BEGIN PROLOGUE IPSORT
%***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 N6A1A, N6A2A
%***TYPE INTEGER (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H)
%***KEYWORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT
%***AUTHOR Jones, R. E., (SNLA)
% Kahaner, D. K., (NBS)
% Rhoads, G. S., (NBS)
% Wisniewski, J. A., (SNLA)
%***DESCRIPTION
%
% IPSORT returns the permutation vector IPERM generated by sorting
% the array IX and, optionally, rearranges the values in IX. IX may
% be sorted in increasing or decreasing order. A slightly modified
% quicksort algorithm is used.
%
% IPERM is such that IX(IPERM(I)) is the Ith value in the
% rearrangement of IX. IPERM may be applied to another array by
% calling IPPERM, SPPERM, DPPERM or HPPERM.
%
% The main difference between IPSORT and its active sorting equivalent
% ISORT is that the data are referenced indirectly rather than
% directly. Therefore, IPSORT should require approximately twice as
% long to execute as ISORT. However, IPSORT is more general.
%
% Description of Parameters
% IX - input/output -- integer array of values to be sorted.
% If ABS(KFLAG) = 2, then the values in IX will be
% rearranged on output; otherwise, they are unchanged.
% N - input -- number of values in array IX to be sorted.
% IPERM - output -- permutation array such that IPERM(I) is the
% index of the value in the original order of the
% IX array that is in the Ith location in the sorted
% order.
% KFLAG - input -- control parameter:
% = 2 means return the permutation vector resulting from
% sorting IX in increasing order and sort IX also.
% = 1 means return the permutation vector resulting from
% sorting IX in increasing order and do not sort IX.
% = -1 means return the permutation vector resulting from
% sorting IX in decreasing order and do not sort IX.
% = -2 means return the permutation vector resulting from
% sorting IX in decreasing order and sort IX 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.
% 810801 Further modified by David K. Kahaner.
% 870423 Modified by Gregory S. Rhoads for passive sorting with the
% option for the rearrangement of the original data.
% 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 IPSORT
% .. Scalar Arguments ..
persistent gt i ij il indx indx0 istrt itemp iu j k kk l lm lmt m nn r ;
if isempty(gt), gt=zeros(1,3); end;
% .. Array Arguments ..
iperm_shape=size(iperm);iperm=reshape(iperm,1,[]);
ix_shape=size(ix);ix=reshape(ix,1,[]);
% .. Local Scalars ..
if isempty(r), r=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(itemp), itemp=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 IPSORT
ier = 0;
nn = fix(n);
if( nn<1 )
ier = 1;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','IPSORT','The number of values to be sorted, N, is not positive.',ier,1);
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
return;
else;
kk = fix(abs(kflag));
if( kk~=1 && kk~=2 )
ier = 2;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','IPSORT','The sort control parameter, KFLAG, is not 2, 1, -1, or -2.',ier,1);
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
return;
else;
%
% 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 )
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
return;
end;
%
% Alter array IX to get decreasing order if needed
%
if( kflag<=-1 )
for i = 1 : nn;
ix(i) = fix(-ix(i));
end; i = fix(nn+1);
end;
%
% Sort IX only
%
m = 1;
i = 1;
j = fix(nn);
r = .375e0;
end;
end;
%
gt(:)=0;
while( true );
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( i==j )
gt(3)=1;
continue;
end;
if( r<=0.5898437e0 )
r = r + 3.90625e-2;
else;
r = r - 0.21875e0;
end;
%
end;
gt(1)=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( ix(iperm(i))>ix(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( ix(iperm(j))<ix(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( ix(iperm(i))>ix(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( ix(iperm(l))<=ix(lm) )
%
% Find an element in the first half of the array which is greater
% than LM
%
while( true );
k = fix(k + 1);
if( ix(iperm(k))>=ix(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;
gt(2)=0;
if( j-i>=1 )
gt(1)=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( ix(iperm(i))>ix(lm) )
k = fix(i);
%
while( true );
iperm(k+1) = fix(iperm(k));
k = fix(k - 1);
%
if( ix(lm)>=ix(iperm(k)) )
break;
end;
end;
iperm(k+1) = fix(lm);
end;
end;
%
% Begin again on another portion of the unsorted array
%
end;
gt(3)=0;
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
gt(2)=1;
end;
%
% Clean up
%
if( kflag<=-1 )
for i = 1 : nn;
ix(i) = fix(-ix(i));
end; i = fix(nn+1);
end;
%
% Rearrange the values of IX 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);
itemp = fix(ix(istrt));
while( iperm(indx)>0 );
ix(indx) = fix(ix(iperm(indx)));
indx0 = fix(indx);
iperm(indx) = fix(-iperm(indx));
indx = fix(abs(iperm(indx)));
end;
ix(indx0) = fix(itemp);
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;
%
iperm_shape=zeros(iperm_shape);iperm_shape(:)=iperm(1:numel(iperm_shape));iperm=iperm_shape;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
end %subroutine ipsort
%DECK ISAMAX
|
|
