← Back to branch summary

~ubuntu-branches/ubuntu/karmic/maxima/karmic

« back to all changes in this revision

Viewing changes to src/numerical/slatec/dqagpe.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-11-13 18:39:14 UTC
  • mto: (2.1.2 hoary) (3.2.1 sid) (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 3.
  • Revision ID: james.westby@ubuntu.com-20041113183914-ttig0evwuatnqosl
Tags: upstream-5.9.1
Import upstream version 5.9.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/numerical/slatec/dqagpe.lisp
 
1
;;; Compiled by f2cl version 2.0 beta 2002-05-06
 
2
;;; 
 
3
;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
 
4
;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
 
5
;;;           (:array-slicing t) (:declare-common nil)
 
6
;;;           (:float-format double-float))
 
7
 
 
8
(in-package "SLATEC")
 
9
 
 
10
 
 
11
(defun dqagpe
 
12
       (f a b npts2 points epsabs epsrel limit result abserr neval ier alist
 
13
        blist rlist elist pts iord level ndin last$)
 
14
  (declare (type (array f2cl-lib:integer4 (*)) ndin level iord)
 
15
   (type (array double-float (*)) pts elist rlist blist alist points)
 
16
   (type f2cl-lib:integer4 last$ ier neval limit npts2)
 
17
   (type double-float abserr result epsrel epsabs b a))
 
18
  (prog ((res3la (make-array 3 :element-type 'double-float))
 
19
         (rlist2 (make-array 52 :element-type 'double-float)) (extrap nil)
 
20
         (noext nil) (i 0) (id 0) (ierro 0) (ind1 0) (ind2 0) (ip1 0)
 
21
         (iroff1 0) (iroff2 0) (iroff3 0) (j 0) (jlow 0) (jupbnd 0) (k 0)
 
22
         (ksgn 0) (ktmin 0) (levcur 0) (levmax 0) (maxerr 0) (f2cl-lib:nint 0)
 
23
         (nintp1 0) (npts 0) (nres 0) (nrmax 0) (numrl2 0) (abseps 0.0)
 
24
         (area 0.0) (area1 0.0) (area12 0.0) (area2 0.0) (a1 0.0) (a2 0.0)
 
25
         (b1 0.0) (b2 0.0) (correc 0.0) (defabs 0.0) (defab1 0.0) (defab2 0.0)
 
26
         (dres 0.0) (epmach 0.0) (erlarg 0.0) (erlast 0.0) (errbnd 0.0)
 
27
         (errmax 0.0) (error1 0.0) (erro12 0.0) (error2 0.0) (errsum 0.0)
 
28
         (ertest 0.0) (oflow 0.0) (resa 0.0) (resabs 0.0) (reseps 0.0)
 
29
         (f2cl-lib:sign 0.0) (temp 0.0) (uflow 0.0) (abs$ 0.0f0))
 
30
    (declare (type single-float abs$)
 
31
     (type (simple-array double-float (52)) rlist2)
 
32
     (type (simple-array double-float (3)) res3la)
 
33
     (type double-float uflow temp f2cl-lib:sign reseps resabs resa oflow
 
34
      ertest errsum error2 erro12 error1 errmax errbnd erlast erlarg epmach
 
35
      dres defab2 defab1 defabs correc b2 b1 a2 a1 area2 area12 area1 area
 
36
      abseps)
 
37
     (type f2cl-lib:integer4 numrl2 nrmax nres npts nintp1 f2cl-lib:nint maxerr
 
38
      levmax levcur ktmin ksgn k jupbnd jlow j iroff3 iroff2 iroff1 ip1 ind2
 
39
      ind1 ierro id i)
 
40
     (type f2cl-lib:logical noext extrap))
 
41
    (setf epmach (f2cl-lib:d1mach 4))
 
42
    (setf ier 0)
 
43
    (setf neval 0)
 
44
    (setf last$ 0)
 
45
    (setf result 0.0)
 
46
    (setf abserr 0.0)
 
47
    (f2cl-lib:fset (f2cl-lib:fref alist (1) ((1 *))) a)
 
48
    (f2cl-lib:fset (f2cl-lib:fref blist (1) ((1 *))) b)
 
49
    (f2cl-lib:fset (f2cl-lib:fref rlist (1) ((1 *))) 0.0)
 
50
    (f2cl-lib:fset (f2cl-lib:fref elist (1) ((1 *))) 0.0)
 
51
    (f2cl-lib:fset (f2cl-lib:fref iord (1) ((1 *))) 0)
 
52
    (f2cl-lib:fset (f2cl-lib:fref level (1) ((1 *))) 0)
 
53
    (setf npts (f2cl-lib:int-sub npts2 2))
 
54
    (if
 
55
     (or (< npts2 2)
 
56
         (<= limit npts)
 
57
         (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29))))
 
58
     (setf ier 6))
 
59
    (if (= ier 6) (go label999))
 
60
    (setf f2cl-lib:sign 1.0)
 
61
    (if (> a b) (setf f2cl-lib:sign -1.0))
 
62
    (f2cl-lib:fset (f2cl-lib:fref pts (1) ((1 *))) (min a b))
 
63
    (if (= npts 0) (go label15))
 
64
    (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 
65
                  ((> i npts) nil)
 
66
      (tagbody
 
67
        (f2cl-lib:fset (f2cl-lib:fref pts ((f2cl-lib:int-add i 1)) ((1 *)))
 
68
                       (f2cl-lib:fref points (i) ((1 *))))
 
69
       label10))
 
70
   label15
 
71
    (f2cl-lib:fset (f2cl-lib:fref pts ((f2cl-lib:int-add npts 2)) ((1 *)))
 
72
                   (max a b))
 
73
    (setf f2cl-lib:nint (f2cl-lib:int-add npts 1))
 
74
    (setf a1 (f2cl-lib:fref pts (1) ((1 *))))
 
75
    (if (= npts 0) (go label40))
 
76
    (setf nintp1 (f2cl-lib:int-add f2cl-lib:nint 1))
 
77
    (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 
78
                  ((> i f2cl-lib:nint) nil)
 
79
      (tagbody
 
80
        (setf ip1 (f2cl-lib:int-add i 1))
 
81
        (f2cl-lib:fdo (j ip1 (f2cl-lib:int-add j 1))
 
82
                      ((> j nintp1) nil)
 
83
          (tagbody
 
84
            (if
 
85
             (<= (f2cl-lib:fref pts (i) ((1 *)))
 
86
                 (f2cl-lib:fref pts (j) ((1 *))))
 
87
             (go label20))
 
88
            (setf temp (f2cl-lib:fref pts (i) ((1 *))))
 
89
            (f2cl-lib:fset (f2cl-lib:fref pts (i) ((1 *)))
 
90
                           (f2cl-lib:fref pts (j) ((1 *))))
 
91
            (f2cl-lib:fset (f2cl-lib:fref pts (j) ((1 *))) temp)))))
 
92
   label20
 
93
    (if
 
94
     (or (/= (f2cl-lib:fref pts (1) ((1 *))) (min a b))
 
95
         (/= (f2cl-lib:fref pts (nintp1) ((1 *))) (max a b)))
 
96
     (setf ier 6))
 
97
    (if (= ier 6) (go label999))
 
98
   label40
 
99
    (setf resabs 0.0)
 
100
    (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 
101
                  ((> i f2cl-lib:nint) nil)
 
102
      (tagbody
 
103
        (setf b1 (f2cl-lib:fref pts ((f2cl-lib:int-add i 1)) ((1 *))))
 
104
        (multiple-value-bind
 
105
            (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 
106
            (dqk21 f a1 b1 area1 error1 defabs resa)
 
107
          (declare (ignore var-0 var-1 var-2))
 
108
          (setf area1 var-3)
 
109
          (setf error1 var-4)
 
110
          (setf defabs var-5)
 
111
          (setf resa var-6))
 
112
        (setf abserr (+ abserr error1))
 
113
        (setf result (+ result area1))
 
114
        (f2cl-lib:fset (f2cl-lib:fref ndin (i) ((1 *))) 0)
 
115
        (if (and (= error1 resa) (/= error1 0.0))
 
116
            (f2cl-lib:fset (f2cl-lib:fref ndin (i) ((1 *))) 1))
 
117
        (setf resabs (+ resabs defabs))
 
118
        (f2cl-lib:fset (f2cl-lib:fref level (i) ((1 *))) 0)
 
119
        (f2cl-lib:fset (f2cl-lib:fref elist (i) ((1 *))) error1)
 
120
        (f2cl-lib:fset (f2cl-lib:fref alist (i) ((1 *))) a1)
 
121
        (f2cl-lib:fset (f2cl-lib:fref blist (i) ((1 *))) b1)
 
122
        (f2cl-lib:fset (f2cl-lib:fref rlist (i) ((1 *))) area1)
 
123
        (f2cl-lib:fset (f2cl-lib:fref iord (i) ((1 *))) i)
 
124
        (setf a1 b1)
 
125
       label50))
 
126
    (setf errsum 0.0)
 
127
    (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 
128
                  ((> i f2cl-lib:nint) nil)
 
129
      (tagbody
 
130
        (if (= (f2cl-lib:fref ndin (i) ((1 *))) 1)
 
131
            (f2cl-lib:fset (f2cl-lib:fref elist (i) ((1 *))) abserr))
 
132
        (setf errsum (+ errsum (f2cl-lib:fref elist (i) ((1 *)))))
 
133
       label55))
 
134
    (setf last$ f2cl-lib:nint)
 
135
    (setf neval (f2cl-lib:int-mul 21 f2cl-lib:nint))
 
136
    (setf dres (coerce (abs result) 'double-float))
 
137
    (setf errbnd (max epsabs (* epsrel dres)))
 
138
    (if (and (<= abserr (* 100.0 epmach resabs)) (> abserr errbnd))
 
139
        (setf ier 2))
 
140
    (if (= f2cl-lib:nint 1) (go label80))
 
141
    (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 
142
                  ((> i npts) nil)
 
143
      (tagbody
 
144
        (setf jlow (f2cl-lib:int-add i 1))
 
145
        (setf ind1 (f2cl-lib:fref iord (i) ((1 *))))
 
146
        (f2cl-lib:fdo (j jlow (f2cl-lib:int-add j 1))
 
147
                      ((> j f2cl-lib:nint) nil)
 
148
          (tagbody
 
149
            (setf ind2 (f2cl-lib:fref iord (j) ((1 *))))
 
150
            (if
 
151
             (> (f2cl-lib:fref elist (ind1) ((1 *)))
 
152
                (f2cl-lib:fref elist (ind2) ((1 *))))
 
153
             (go label60))
 
154
            (setf ind1 ind2)
 
155
            (setf k j)
 
156
           label60))
 
157
        (if (= ind1 (f2cl-lib:fref iord (i) ((1 *)))) (go label70))
 
158
        (f2cl-lib:fset (f2cl-lib:fref iord (k) ((1 *)))
 
159
                       (f2cl-lib:fref iord (i) ((1 *))))
 
160
        (f2cl-lib:fset (f2cl-lib:fref iord (i) ((1 *))) ind1)
 
161
       label70))
 
162
    (if (< limit npts2) (setf ier 1))
 
163
   label80
 
164
    (if (or (/= ier 0) (<= abserr errbnd)) (go label999))
 
165
    (f2cl-lib:fset (f2cl-lib:fref rlist2 (1) ((1 52))) result)
 
166
    (setf maxerr (f2cl-lib:fref iord (1) ((1 *))))
 
167
    (setf errmax (f2cl-lib:fref elist (maxerr) ((1 *))))
 
168
    (setf area result)
 
169
    (setf nrmax 1)
 
170
    (setf nres 0)
 
171
    (setf numrl2 1)
 
172
    (setf ktmin 0)
 
173
    (setf extrap f2cl-lib:%false%)
 
174
    (setf noext f2cl-lib:%false%)
 
175
    (setf erlarg errsum)
 
176
    (setf ertest errbnd)
 
177
    (setf levmax 1)
 
178
    (setf iroff1 0)
 
179
    (setf iroff2 0)
 
180
    (setf iroff3 0)
 
181
    (setf ierro 0)
 
182
    (setf uflow (f2cl-lib:d1mach 1))
 
183
    (setf oflow (f2cl-lib:d1mach 2))
 
184
    (setf abserr oflow)
 
185
    (setf ksgn -1)
 
186
    (if (>= dres (* (- 1.0 (* 50.0 epmach)) resabs)) (setf ksgn 1))
 
187
    (f2cl-lib:fdo (last$ npts2 (f2cl-lib:int-add last$ 1))
 
188
                  ((> last$ limit) nil)
 
189
      (tagbody
 
190
        (setf levcur
 
191
                (f2cl-lib:int-add (f2cl-lib:fref level (maxerr) ((1 *))) 1))
 
192
        (setf a1 (f2cl-lib:fref alist (maxerr) ((1 *))))
 
193
        (setf b1
 
194
                (* 0.5
 
195
                   (+ (f2cl-lib:fref alist (maxerr) ((1 *)))
 
196
                      (f2cl-lib:fref blist (maxerr) ((1 *))))))
 
197
        (setf a2 b1)
 
198
        (setf b2 (f2cl-lib:fref blist (maxerr) ((1 *))))
 
199
        (setf erlast errmax)
 
200
        (multiple-value-bind
 
201
            (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 
202
            (dqk21 f a1 b1 area1 error1 resa defab1)
 
203
          (declare (ignore var-0 var-1 var-2))
 
204
          (setf area1 var-3)
 
205
          (setf error1 var-4)
 
206
          (setf resa var-5)
 
207
          (setf defab1 var-6))
 
208
        (multiple-value-bind
 
209
            (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 
210
            (dqk21 f a2 b2 area2 error2 resa defab2)
 
211
          (declare (ignore var-0 var-1 var-2))
 
212
          (setf area2 var-3)
 
213
          (setf error2 var-4)
 
214
          (setf resa var-5)
 
215
          (setf defab2 var-6))
 
216
        (setf neval (f2cl-lib:int-add neval 42))
 
217
        (setf area12 (+ area1 area2))
 
218
        (setf erro12 (+ error1 error2))
 
219
        (setf errsum (- (+ errsum erro12) errmax))
 
220
        (setf area (- (+ area area12) (f2cl-lib:fref rlist (maxerr) ((1 *)))))
 
221
        (if (or (= defab1 error1) (= defab2 error2)) (go label95))
 
222
        (if
 
223
         (or
 
224
          (> (abs (- (f2cl-lib:fref rlist (maxerr) ((1 *))) area12))
 
225
             (* 1.0e-5 (abs area12)))
 
226
          (< erro12 (* 0.99 errmax)))
 
227
         (go label90))
 
228
        (if extrap (setf iroff2 (f2cl-lib:int-add iroff2 1)))
 
229
        (if (not extrap) (setf iroff1 (f2cl-lib:int-add iroff1 1)))
 
230
       label90
 
231
        (if (and (> last$ 10) (> erro12 errmax))
 
232
            (setf iroff3 (f2cl-lib:int-add iroff3 1)))
 
233
       label95
 
234
        (f2cl-lib:fset (f2cl-lib:fref level (maxerr) ((1 *))) levcur)
 
235
        (f2cl-lib:fset (f2cl-lib:fref level (last$) ((1 *))) levcur)
 
236
        (f2cl-lib:fset (f2cl-lib:fref rlist (maxerr) ((1 *))) area1)
 
237
        (f2cl-lib:fset (f2cl-lib:fref rlist (last$) ((1 *))) area2)
 
238
        (setf errbnd (max epsabs (* epsrel (abs area))))
 
239
        (if (or (>= (f2cl-lib:int-add iroff1 iroff2) 10) (>= iroff3 20))
 
240
            (setf ier 2))
 
241
        (if (>= iroff2 5) (setf ierro 3))
 
242
        (if (= last$ limit) (setf ier 1))
 
243
        (if
 
244
         (<= (max (abs a1) (abs b2))
 
245
             (* (+ 1.0 (* 100.0 epmach)) (+ (abs a2) (* 1000.0 uflow))))
 
246
         (setf ier 4))
 
247
        (if (> error2 error1) (go label100))
 
248
        (f2cl-lib:fset (f2cl-lib:fref alist (last$) ((1 *))) a2)
 
249
        (f2cl-lib:fset (f2cl-lib:fref blist (maxerr) ((1 *))) b1)
 
250
        (f2cl-lib:fset (f2cl-lib:fref blist (last$) ((1 *))) b2)
 
251
        (f2cl-lib:fset (f2cl-lib:fref elist (maxerr) ((1 *))) error1)
 
252
        (f2cl-lib:fset (f2cl-lib:fref elist (last$) ((1 *))) error2)
 
253
        (go label110)
 
254
       label100
 
255
        (f2cl-lib:fset (f2cl-lib:fref alist (maxerr) ((1 *))) a2)
 
256
        (f2cl-lib:fset (f2cl-lib:fref alist (last$) ((1 *))) a1)
 
257
        (f2cl-lib:fset (f2cl-lib:fref blist (last$) ((1 *))) b1)
 
258
        (f2cl-lib:fset (f2cl-lib:fref rlist (maxerr) ((1 *))) area2)
 
259
        (f2cl-lib:fset (f2cl-lib:fref rlist (last$) ((1 *))) area1)
 
260
        (f2cl-lib:fset (f2cl-lib:fref elist (maxerr) ((1 *))) error2)
 
261
        (f2cl-lib:fset (f2cl-lib:fref elist (last$) ((1 *))) error1)
 
262
       label110
 
263
        (multiple-value-bind
 
264
            (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 
265
            (dqpsrt limit last$ maxerr errmax elist iord nrmax)
 
266
          (declare (ignore var-0 var-1 var-4 var-5))
 
267
          (setf maxerr var-2)
 
268
          (setf errmax var-3)
 
269
          (setf nrmax var-6))
 
270
        (if (<= errsum errbnd) (go label190))
 
271
        (if (/= ier 0) (go label170))
 
272
        (if noext (go label160))
 
273
        (setf erlarg (- erlarg erlast))
 
274
        (if (<= (f2cl-lib:int-add levcur 1) levmax)
 
275
            (setf erlarg (+ erlarg erro12)))
 
276
        (if extrap (go label120))
 
277
        (if
 
278
         (<= (f2cl-lib:int-add (f2cl-lib:fref level (maxerr) ((1 *))) 1)
 
279
             levmax)
 
280
         (go label160))
 
281
        (setf extrap f2cl-lib:%true%)
 
282
        (setf nrmax 2)
 
283
       label120
 
284
        (if (or (= ierro 3) (<= erlarg ertest)) (go label140))
 
285
        (setf id nrmax)
 
286
        (setf jupbnd last$)
 
287
        (if (> last$ (+ 2 (the f2cl-lib:integer4 (truncate limit 2))))
 
288
            (setf jupbnd (f2cl-lib:int-sub (f2cl-lib:int-add limit 3) last$)))
 
289
        (f2cl-lib:fdo (k id (f2cl-lib:int-add k 1))
 
290
                      ((> k jupbnd) nil)
 
291
          (tagbody
 
292
            (setf maxerr (f2cl-lib:fref iord (nrmax) ((1 *))))
 
293
            (setf errmax (f2cl-lib:fref elist (maxerr) ((1 *))))
 
294
            (if
 
295
             (<= (f2cl-lib:int-add (f2cl-lib:fref level (maxerr) ((1 *))) 1)
 
296
                 levmax)
 
297
             (go label160))
 
298
            (setf nrmax (f2cl-lib:int-add nrmax 1))
 
299
           label130))
 
300
       label140
 
301
        (setf numrl2 (f2cl-lib:int-add numrl2 1))
 
302
        (f2cl-lib:fset (f2cl-lib:fref rlist2 (numrl2) ((1 52))) area)
 
303
        (if (<= numrl2 2) (go label155))
 
304
        (multiple-value-bind
 
305
            (var-0 var-1 var-2 var-3 var-4 var-5)
 
306
            (dqelg numrl2 rlist2 reseps abseps res3la nres)
 
307
          (declare (ignore var-1 var-4))
 
308
          (setf numrl2 var-0)
 
309
          (setf reseps var-2)
 
310
          (setf abseps var-3)
 
311
          (setf nres var-5))
 
312
        (setf ktmin (f2cl-lib:int-add ktmin 1))
 
313
        (if (and (> ktmin 5) (< abserr (* 0.001 errsum))) (setf ier 5))
 
314
        (if (>= abseps abserr) (go label150))
 
315
        (setf ktmin 0)
 
316
        (setf abserr abseps)
 
317
        (setf result reseps)
 
318
        (setf correc erlarg)
 
319
        (setf ertest (max epsabs (* epsrel (abs reseps))))
 
320
        (if (< abserr ertest) (go label170))
 
321
       label150
 
322
        (if (= numrl2 1) (setf noext f2cl-lib:%true%))
 
323
        (if (>= ier 5) (go label170))
 
324
       label155
 
325
        (setf maxerr (f2cl-lib:fref iord (1) ((1 *))))
 
326
        (setf errmax (f2cl-lib:fref elist (maxerr) ((1 *))))
 
327
        (setf nrmax 1)
 
328
        (setf extrap f2cl-lib:%false%)
 
329
        (setf levmax (f2cl-lib:int-add levmax 1))
 
330
        (setf erlarg errsum)
 
331
       label160))
 
332
   label170
 
333
    (if (= abserr oflow) (go label190))
 
334
    (if (= (f2cl-lib:int-add ier ierro) 0) (go label180))
 
335
    (if (= ierro 3) (setf abserr (+ abserr correc)))
 
336
    (if (= ier 0) (setf ier 3))
 
337
    (if (and (/= result 0.0) (/= area 0.0)) (go label175))
 
338
    (if (> abserr errsum) (go label190))
 
339
    (if (= area 0.0) (go label210))
 
340
    (go label180)
 
341
   label175
 
342
    (if (> (/ abserr (abs result)) (/ errsum (abs area))) (go label190))
 
343
   label180
 
344
    (if
 
345
     (and (= ksgn -1)
 
346
          (<= (max (abs result) (abs area)) (* resabs 0.010000000000000002)))
 
347
     (go label210))
 
348
    (if
 
349
     (or (> 0.010000000000000002 (/ result area))
 
350
         (> (/ result area) 100.0)
 
351
         (> errsum (abs area)))
 
352
     (setf ier 6))
 
353
    (go label210)
 
354
   label190
 
355
    (setf result 0.0)
 
356
    (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
357
                  ((> k last$) nil)
 
358
      (tagbody
 
359
        (setf result (+ result (f2cl-lib:fref rlist (k) ((1 *)))))
 
360
       label200))
 
361
    (setf abserr errsum)
 
362
   label210
 
363
    (if (> ier 2) (setf ier (f2cl-lib:int-sub ier 1)))
 
364
    (setf result (* result f2cl-lib:sign))
 
365
   label999
 
366
    (go end_label)
 
367
   end_label
 
368
    (return
 
369
     (values nil
 
370
             nil
 
371
             nil
 
372
             nil
 
373
             nil
 
374
             nil
 
375
             nil
 
376
             nil
 
377
             result
 
378
             abserr
 
379
             neval
 
380
             ier
 
381
             nil
 
382
             nil
 
383
             nil
 
384
             nil
 
385
             nil
 
386
             nil
 
387
             nil
 
388
             nil
 
389
             last$))))
 
390