[limit,last,maxerr,ermax,elist,iord,nrmax]=dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax);

function [limit,last,maxerr,ermax,elist,iord,nrmax]=dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax);
%***BEGIN PROLOGUE  DQPSRT
%***SUBSIDIARY
%***PURPOSE  This routine maintains the descending ordering in the
%            list of the local error estimated resulting from the
%            interval subdivision process. At each call two error
%            estimates are inserted using the sequential search
%            method, top-down for the largest error estimate and
%            bottom-up for the smallest error estimate.
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (QPSRT-S, DQPSRT-D)
%***KEYWORDS  SEQUENTIAL SORTING
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%           Ordering routine
%           Standard fortran subroutine
%           doubleprecision version
%
%           PARAMETERS (MEANING AT OUTPUT)
%              LIMIT  - Integer
%                       Maximum number of error estimates the list
%                       can contain
%
%              LAST   - Integer
%                       Number of error estimates currently in the list
%
%              MAXERR - Integer
%                       MAXERR points to the NRMAX-th largest error
%                       estimate currently in the list
%
%              ERMAX  - doubleprecision
%                       NRMAX-th largest error estimate
%                       ERMAX = ELIST(MAXERR)
%
%              ELIST  - doubleprecision
%                       Vector of dimension LAST containing
%                       the error estimates
%
%              IORD   - Integer
%                       Vector of dimension LAST, the first K elements
%                       of which contain pointers to the error
%                       estimates, such that
%                       ELIST(IORD(1)),...,  ELIST(IORD(K))
%                       form a decreasing sequence, with
%                       K = LAST if LAST.LE.(LIMIT/2+2), and
%                       K = LIMIT+1-LAST otherwise
%
%              NRMAX  - Integer
%                       MAXERR = IORD(NRMAX)
%
%***SEE ALSO  DQAGE, DQAGIE, DQAGPE, DQAWSE
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800101  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DQPSRT
%
persistent errmax errmin gt i ibeg ido isucc j jbnd jupbn k ; 

if isempty(errmax), errmax=0; end;
if isempty(errmin), errmin=0; end;
if isempty(i), i=0; end;
if isempty(ibeg), ibeg=0; end;
if isempty(ido), ido=0; end;
if isempty(isucc), isucc=0; end;
if isempty(j), j=0; end;
if isempty(jbnd), jbnd=0; end;
if isempty(jupbn), jupbn=0; end;
if isempty(k), k=0; end;
if isempty(gt), gt=0; end;
elist_shape=size(elist);elist=reshape(elist,1,[]);
iord_shape=size(iord);iord=reshape(iord,1,[]);
%
%           CHECK WHETHER THE LIST CONTAINS MORE THAN
%           TWO ERROR ESTIMATES.
%
%***FIRST EXECUTABLE STATEMENT  DQPSRT
gt=0;
if( last>2 )
%
%           THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A
%           DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR
%           ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE SHOULD
%           START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE.
%
errmax = elist(maxerr);
if( nrmax~=1 )
ido = fix(nrmax - 1);
for i = 1 : ido;
isucc = fix(iord(nrmax-1));
% ***JUMP OUT OF DO-LOOP
if( errmax<=elist(isucc) )
break;
end;
iord(nrmax) = fix(isucc);
nrmax = fix(nrmax - 1);
end;
end;
%
%           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO BE MAINTAINED
%           IN DESCENDING ORDER. THIS NUMBER DEPENDS ON THE NUMBER OF
%           SUBDIVISIONS STILL ALLOWED.
%
jupbn = fix(last);
if( last>(fix(limit./2)+2) )
jupbn = fix(limit + 3 - last);
end;
errmin = elist(last);
%
%           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
%           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
%
jbnd = fix(jupbn - 1);
ibeg = fix(nrmax + 1);
if( ibeg<=jbnd )
for i = ibeg : jbnd;
isucc = fix(iord(i));
% ***JUMP OUT OF DO-LOOP
if( errmax>=elist(isucc) )
%
%           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
%
iord(i-1) = fix(maxerr);
k = fix(jbnd);
for j = i : jbnd;
isucc = fix(iord(k));
% ***JUMP OUT OF DO-LOOP
if( errmin<elist(isucc) )
iord(k+1) = fix(last);
gt=1;
break;
else;
iord(k+1) = fix(isucc);
k = fix(k - 1);
end;
end;
if(gt==0)
iord(i) = fix(last);
end;
gt=1;
break;
else;
iord(i-1) = fix(isucc);
end;
end;
end;
if(gt==0)
iord(jbnd) = fix(maxerr);
iord(jupbn) = fix(last);
end;
else;
iord(1) = 1;
iord(2) = 2;
end;
%
%           SET MAXERR AND ERMAX.
%
maxerr = fix(iord(nrmax));
ermax = elist(maxerr);
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
end %subroutine dqpsrt
%DECK DQRDC