[ia,ja,a,n,kflag]=qs2i1r(ia,ja,a,n,kflag);

function [ia,ja,a,n,kflag]=qs2i1r(ia,ja,a,n,kflag);
%***BEGIN PROLOGUE  QS2I1R
%***SUBSIDIARY
%***PURPOSE  Sort an integer array, moving an integer and real array.
%            This routine sorts the integer array IA and makes the same
%            interchanges in the integer array JA and the real array A.
%            The array IA may be sorted in increasing order or decreas-
%            ing order.  A slightly modified QUICKSORT algorithm is
%            used.
%***LIBRARY   SLATEC (SLAP)
%***CATEGORY  N6A2A
%***TYPE      SINGLE PRECISION (QS2I1R-S, QS2I1D-D)
%***KEYWORDS  SINGLETON QUICKSORT, SLAP, SORT, SORTING
%***AUTHOR  Jones, R. E., (SNLA)
%           Kahaner, D. K., (NBS)
%           Seager, M. K., (LLNL) [e-mail]
%           Wisniewski, J. A., (SNLA)
%***DESCRIPTION
%     Written by Rondall E Jones
%     Modified by John A. Wisniewski to use the Singleton QUICKSORT
%     algorithm. date 18 November 1976.
%
%     Further modified by David K. Kahaner
%     National Bureau of Standards
%     August, 1981
%
%     Even further modification made to bring the code up to the
%     Fortran 77 level and make it more readable and to carry
%     along one integer array and one real array during the sort by
%     Mark K. Seager
%     Lawrence Livermore National Laboratory
%     November, 1987
%     This routine was adapted from the ISORT routine.
%
%     ABSTRACT
%         This routine sorts an integer array IA and makes the same
%         interchanges in the integer array JA and the real array A.
%         The array IA may be sorted in increasing order or decreasing
%         order.  A slightly modified quicksort algorithm is used.
%
%     DESCRIPTION OF PARAMETERS
%        IA - Integer array of values to be sorted.
%        JA - Integer array to be carried along.
%         A - Real array to be carried along.
%         N - Number of values in integer array IA to be sorted.
%     KFLAG - Control parameter
%           = 1 means sort IA in INCREASING order.
%           =-1 means sort IA in DECREASING order.
%
%***SEE ALSO  SS2Y
%***REFERENCES  R. C. Singleton, Algorithm 347, An Efficient Algorithm
%                 for Sorting With Minimal Storage, Communications ACM
%                 12:3 (1969), pp.185-7.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   761118  DATE WRITTEN
%   890125  Previous REVISION DATE
%   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
%   890922  Numerous changes to prologue to make closer to SLATEC
%           standard.  (FNF)
%   890929  Numerous changes to reduce SP/DP differences.  (FNF)
%   900805  Changed XERROR calls to calls to XERMSG.  (RWC)
%   910411  Prologue converted to Version 4.0 format.  (BAB)
%   910506  Made subsidiary to SS2Y and corrected reference.  (FNF)
%   920511  Added complete declaration section.  (WRB)
%   920929  Corrected format of reference.  (FNF)
%   921012  Added E0's to f.p. constants.  (FNF)
%***end PROLOGUE  QS2I1R
%VD$R NOVECTOR
%VD$R NOCONCUR
%     .. Scalar Arguments ..
%     .. Array Arguments ..
%     .. Local Scalars ..
persistent gt i iit ij il it iu j jjt jt k kk l m nn r ta tta ; 

if isempty(r), r=0; end;
if isempty(ta), ta=0; end;
if isempty(tta), tta=0; end;
if isempty(i), i=0; end;
if isempty(iit), iit=0; end;
if isempty(ij), ij=0; end;
if isempty(it), it=0; end;
if isempty(j), j=0; end;
if isempty(jjt), jjt=0; end;
if isempty(jt), jt=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;
if isempty(gt), gt=zeros(1,2); 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  QS2I1R
nn = fix(n);
gt(:)=0;
if( nn<1 )
xermsg('SLATEC','QS2I1R','The number of values to be sorted was not positive.',1,1);
elseif( n~=1 ) ;
kk = fix(abs(kflag));
if( kk==1 )
%
%     Alter array IA to get decreasing order if needed.
%
if( kflag<1 )
for i = 1 : nn;
ia(i) = fix(-ia(i));
end; i = fix(nn+1);
end;
%
%     Sort IA and carry JA and A along.
%     And now...Just a little black magic...
m = 1;
i = 1;
j = fix(nn);
r = .375e0;
while( true );
if( r<=0.5898437e0 )
r = r + 3.90625e-2;
else;
r = r - .21875e0;
end;
while( true );
k = fix(i);
%
%     Select a central element of the array and savemlv it in location
%     it, jt, at.
%
ij = fix(i + fix((j-i).*r));
it = fix(ia(ij));
jt = fix(ja(ij));
ta = a(ij);
%
%     If first element of array is greater than it, interchange with it.
%
if( ia(i)>it )
ia(ij) = fix(ia(i));
ia(i) = fix(it);
it = fix(ia(ij));
ja(ij) = fix(ja(i));
ja(i) = fix(jt);
jt = fix(ja(ij));
a(ij) = a(i);
a(i) = ta;
ta = a(ij);
end;
l = fix(j);
%
%     If last element of array is less than it, swap with it.
%
if( ia(j)<it )
ia(ij) = fix(ia(j));
ia(j) = fix(it);
it = fix(ia(ij));
ja(ij) = fix(ja(j));
ja(j) = fix(jt);
jt = fix(ja(ij));
a(ij) = a(j);
a(j) = ta;
ta = a(ij);
%
%     If first element of array is greater than it, swap with it.
%
if( ia(i)>it )
ia(ij) = fix(ia(i));
ia(i) = fix(it);
it = fix(ia(ij));
ja(ij) = fix(ja(i));
ja(i) = fix(jt);
jt = fix(ja(ij));
a(ij) = a(i);
a(i) = ta;
ta = a(ij);
end;
end;
%
%     Find an element in the second half of the array which is
%     smaller than it.
%
while( true );
l = fix(l - 1);
if( ia(l)<=it )
%
%     Find an element in the first half of the array which is
%     greater than it.
%
while( true );
k = fix(k + 1);
if( ia(k)>=it )
break;
end;
end;
%
%     Interchange these elements.
%
if( k>l )
break;
end;
iit = fix(ia(l));
ia(l) = fix(ia(k));
ia(k) = fix(iit);
jjt = fix(ja(l));
ja(l) = fix(ja(k));
ja(k) = fix(jjt);
tta = a(l);
a(l) = a(k);
a(k) = tta;
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;
gt(:)=0;
while( j-i<1 );
if( i~=j )
if( i==1 )
gt(1)=1;
break;
end;
i = fix(i - 1);
while( true );
i = fix(i + 1);
if( i==j )
break;
end;
it = fix(ia(i+1));
jt = fix(ja(i+1));
ta = a(i+1);
if( ia(i)>it )
k = fix(i);
while( true );
ia(k+1) = fix(ia(k));
ja(k+1) = fix(ja(k));
a(k+1) = a(k);
k = fix(k - 1);
if( it>=ia(k) )
break;
end;
end;
ia(k+1) = fix(it);
ja(k+1) = fix(jt);
a(k+1) = ta;
end;
end;
end;
%
%     Begin again on another portion of the unsorted array.
%
m = fix(m - 1);
if( m==0 )
gt(2)=1;
break;
end;
i = fix(il(m));
j = fix(iu(m));
end;
if(any(gt==1))
break;
end;
end;
gt(1)=0;
if(gt(2)==1)
break;
end;
end;
%
%     Clean up, if necessary.
%
if( kflag<1 )
for i = 1 : nn;
ia(i) = fix(-ia(i));
end; i = fix(nn+1);
end;
else;
xermsg('SLATEC','QS2I1R','The sort control parameter, K, was not 1 or -1.',2,1);
end;
end;
%------------- LAST LINE OF QS2I1R FOLLOWS ----------------------------
end %subroutine qs2i1r
%DECK QWGTC