dqpsrt.f
4.24 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
!
! fvn comment :
! Unmodified quadpack routine from http://www.netlib.org/quadpack
!
subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
!***begin prologue dqpsrt
!***refer to dqage,dqagie,dqagpe,dqawse
!***routines called (none)
!***revision date 810101 (yymmdd)
!***keywords sequential sorting
!***author piessens,robert,appl. math. & progr. div. - k.u.leuven
! de doncker,elise,appl. math. & progr. div. - k.u.leuven
!***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.
!***description
!
! ordering routine
! standard fortran subroutine
! double precision 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 - double precision
! nrmax-th largest error estimate
! ermax = elist(maxerr)
!
! elist - double precision
! 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)
!
!***end prologue dqpsrt
!
double precision elist,ermax,errmax,errmin
integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr, &
nrmax
dimension elist(last),iord(last)
!
! check whether the list contains more than
! two error estimates.
!
!***first executable statement dqpsrt
if(last.gt.2) go to 10
iord(1) = 1
iord(2) = 2
go to 90
!
! 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.
!
10 errmax = elist(maxerr)
if(nrmax.eq.1) go to 30
ido = nrmax-1
do 20 i = 1,ido
isucc = iord(nrmax-1)
! ***jump out of do-loop
if(errmax.le.elist(isucc)) go to 30
iord(nrmax) = isucc
nrmax = nrmax-1
20 continue
!
! compute the number of elements in the list to be maintained
! in descending order. this number depends on the number of
! subdivisions still allowed.
!
30 jupbn = last
if(last.gt.(limit/2+2)) jupbn = limit+3-last
errmin = elist(last)
!
! insert errmax by traversing the list top-down,
! starting comparison from the element elist(iord(nrmax+1)).
!
jbnd = jupbn-1
ibeg = nrmax+1
if(ibeg.gt.jbnd) go to 50
do 40 i=ibeg,jbnd
isucc = iord(i)
! ***jump out of do-loop
if(errmax.ge.elist(isucc)) go to 60
iord(i-1) = isucc
40 continue
50 iord(jbnd) = maxerr
iord(jupbn) = last
go to 90
!
! insert errmin by traversing the list bottom-up.
!
60 iord(i-1) = maxerr
k = jbnd
do 70 j=i,jbnd
isucc = iord(k)
! ***jump out of do-loop
if(errmin.lt.elist(isucc)) go to 80
iord(k+1) = isucc
k = k-1
70 continue
iord(i) = last
go to 90
80 iord(k+1) = last
!
! set maxerr and ermax.
!
90 maxerr = iord(nrmax)
ermax = elist(maxerr)
return
end subroutine