← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/qawfe.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/qawfe.f
 
1
*DECK QAWFE
 
2
      SUBROUTINE QAWFE (F, A, OMEGA, INTEGR, EPSABS, LIMLST, LIMIT,
 
3
     +   MAXP1, RESULT, ABSERR, NEVAL, IER, RSLST, ERLST, IERLST, LST,
 
4
     +   ALIST, BLIST, RLIST, ELIST, IORD, NNLOG, CHEBMO)
 
5
C***BEGIN PROLOGUE  QAWFE
 
6
C***PURPOSE  The routine calculates an approximation result to a
 
7
C            given Fourier integral
 
8
C            I = Integral of F(X)*W(X) over (A,INFINITY)
 
9
C             where W(X) = COS(OMEGA*X) or W(X) = SIN(OMEGA*X),
 
10
C            hopefully satisfying following claim for accuracy
 
11
C            ABS(I-RESULT).LE.EPSABS.
 
12
C***LIBRARY   SLATEC (QUADPACK)
 
13
C***CATEGORY  H2A3A1
 
14
C***TYPE      SINGLE PRECISION (QAWFE-S, DQAWFE-D)
 
15
C***KEYWORDS  AUTOMATIC INTEGRATOR, CONVERGENCE ACCELERATION,
 
16
C             FOURIER INTEGRALS, INTEGRATION BETWEEN ZEROS, QUADPACK,
 
17
C             QUADRATURE, SPECIAL-PURPOSE INTEGRAL
 
18
C***AUTHOR  Piessens, Robert
 
19
C             Applied Mathematics and Programming Division
 
20
C             K. U. Leuven
 
21
C           de Doncker, Elise
 
22
C             Applied Mathematics and Programming Division
 
23
C             K. U. Leuven
 
24
C***DESCRIPTION
 
25
C
 
26
C        Computation of Fourier integrals
 
27
C        Standard fortran subroutine
 
28
C        Real version
 
29
C
 
30
C        PARAMETERS
 
31
C         ON ENTRY
 
32
C            F      - Real
 
33
C                     Function subprogram defining the integrand
 
34
C                     Function F(X). The actual name for F needs to
 
35
C                     be declared E X T E R N A L in the driver program.
 
36
C
 
37
C            A      - Real
 
38
C                     Lower limit of integration
 
39
C
 
40
C            OMEGA  - Real
 
41
C                     Parameter in the WEIGHT function
 
42
C
 
43
C            INTEGR - Integer
 
44
C                     Indicates which WEIGHT function is used
 
45
C                     INTEGR = 1      W(X) = COS(OMEGA*X)
 
46
C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
 
47
C                     If INTEGR.NE.1.AND.INTEGR.NE.2, the routine will
 
48
C                     end with IER = 6.
 
49
C
 
50
C            EPSABS - Real
 
51
C                     absolute accuracy requested, EPSABS.GT.0
 
52
C                     If EPSABS.LE.0, the routine will end with IER = 6.
 
53
C
 
54
C            LIMLST - Integer
 
55
C                     LIMLST gives an upper bound on the number of
 
56
C                     cycles, LIMLST.GE.1.
 
57
C                     If LIMLST.LT.3, the routine will end with IER = 6.
 
58
C
 
59
C            LIMIT  - Integer
 
60
C                     Gives an upper bound on the number of subintervals
 
61
C                     allowed in the partition of each cycle, LIMIT.GE.1
 
62
C                     each cycle, LIMIT.GE.1.
 
63
C
 
64
C            MAXP1  - Integer
 
65
C                     Gives an upper bound on the number of
 
66
C                     Chebyshev moments which can be stored, I.E.
 
67
C                     for the intervals of lengths ABS(B-A)*2**(-L),
 
68
C                     L=0,1, ..., MAXP1-2, MAXP1.GE.1
 
69
C
 
70
C         ON RETURN
 
71
C            RESULT - Real
 
72
C                     Approximation to the integral X
 
73
C
 
74
C            ABSERR - Real
 
75
C                     Estimate of the modulus of the absolute error,
 
76
C                     which should equal or exceed ABS(I-RESULT)
 
77
C
 
78
C            NEVAL  - Integer
 
79
C                     Number of integrand evaluations
 
80
C
 
81
C            IER    - IER = 0 Normal and reliable termination of
 
82
C                             the routine. It is assumed that the
 
83
C                             requested accuracy has been achieved.
 
84
C                     IER.GT.0 Abnormal termination of the routine. The
 
85
C                             estimates for integral and error are less
 
86
C                             reliable. It is assumed that the requested
 
87
C                             accuracy has not been achieved.
 
88
C            ERROR MESSAGES
 
89
C                    If OMEGA.NE.0
 
90
C                     IER = 1 Maximum number of  cycles  allowed
 
91
C                             Has been achieved., i.e. of subintervals
 
92
C                             (A+(K-1)C,A+KC) where
 
93
C                             C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
 
94
C                             for K = 1, 2, ..., LST.
 
95
C                             One can allow more cycles by increasing
 
96
C                             the value of LIMLST (and taking the
 
97
C                             according dimension adjustments into
 
98
C                             account).
 
99
C                             Examine the array IWORK which contains
 
100
C                             the error flags on the cycles, in order to
 
101
C                             look for eventual local integration
 
102
C                             difficulties. If the position of a local
 
103
C                             difficulty can be determined (e.g.
 
104
C                             SINGULARITY, DISCONTINUITY within the
 
105
C                             interval) one will probably gain from
 
106
C                             splitting up the interval at this point
 
107
C                             and calling appropriate integrators on
 
108
C                             the subranges.
 
109
C                         = 4 The extrapolation table constructed for
 
110
C                             convergence acceleration of the series
 
111
C                             formed by the integral contributions over
 
112
C                             the cycles, does not converge to within
 
113
C                             the requested accuracy. As in the case of
 
114
C                             IER = 1, it is advised to examine the
 
115
C                             array IWORK which contains the error
 
116
C                             flags on the cycles.
 
117
C                         = 6 The input is invalid because
 
118
C                             (INTEGR.NE.1 AND INTEGR.NE.2) or
 
119
C                              EPSABS.LE.0 or LIMLST.LT.3.
 
120
C                              RESULT, ABSERR, NEVAL, LST are set
 
121
C                              to zero.
 
122
C                         = 7 Bad integrand behaviour occurs within one
 
123
C                             or more of the cycles. Location and type
 
124
C                             of the difficulty involved can be
 
125
C                             determined from the vector IERLST. Here
 
126
C                             LST is the number of cycles actually
 
127
C                             needed (see below).
 
128
C                             IERLST(K) = 1 The maximum number of
 
129
C                                           subdivisions (= LIMIT) has
 
130
C                                           been achieved on the K th
 
131
C                                           cycle.
 
132
C                                       = 2 Occurrence of roundoff error
 
133
C                                           is detected and prevents the
 
134
C                                           tolerance imposed on the
 
135
C                                           K th cycle, from being
 
136
C                                           achieved.
 
137
C                                       = 3 Extremely bad integrand
 
138
C                                           behaviour occurs at some
 
139
C                                           points of the K th cycle.
 
140
C                                       = 4 The integration procedure
 
141
C                                           over the K th cycle does
 
142
C                                           not converge (to within the
 
143
C                                           required accuracy) due to
 
144
C                                           roundoff in the
 
145
C                                           extrapolation procedure
 
146
C                                           invoked on this cycle. It
 
147
C                                           is assumed that the result
 
148
C                                           on this interval is the
 
149
C                                           best which can be obtained.
 
150
C                                       = 5 The integral over the K th
 
151
C                                           cycle is probably divergent
 
152
C                                           or slowly convergent. It
 
153
C                                           must be noted that
 
154
C                                           divergence can occur with
 
155
C                                           any other value of
 
156
C                                           IERLST(K).
 
157
C                    If OMEGA = 0 and INTEGR = 1,
 
158
C                    The integral is calculated by means of DQAGIE
 
159
C                    and IER = IERLST(1) (with meaning as described
 
160
C                    for IERLST(K), K = 1).
 
161
C
 
162
C            RSLST  - Real
 
163
C                     Vector of dimension at least LIMLST
 
164
C                     RSLST(K) contains the integral contribution
 
165
C                     over the interval (A+(K-1)C,A+KC) where
 
166
C                     C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
 
167
C                     K = 1, 2, ..., LST.
 
168
C                     Note that, if OMEGA = 0, RSLST(1) contains
 
169
C                     the value of the integral over (A,INFINITY).
 
170
C
 
171
C            ERLST  - Real
 
172
C                     Vector of dimension at least LIMLST
 
173
C                     ERLST(K) contains the error estimate corresponding
 
174
C                     with RSLST(K).
 
175
C
 
176
C            IERLST - Integer
 
177
C                     Vector of dimension at least LIMLST
 
178
C                     IERLST(K) contains the error flag corresponding
 
179
C                     with RSLST(K). For the meaning of the local error
 
180
C                     flags see description of output parameter IER.
 
181
C
 
182
C            LST    - Integer
 
183
C                     Number of subintervals needed for the integration
 
184
C                     If OMEGA = 0 then LST is set to 1.
 
185
C
 
186
C            ALIST, BLIST, RLIST, ELIST - Real
 
187
C                     vector of dimension at least LIMIT,
 
188
C
 
189
C            IORD, NNLOG - Integer
 
190
C                     Vector of dimension at least LIMIT, providing
 
191
C                     space for the quantities needed in the subdivision
 
192
C                     process of each cycle
 
193
C
 
194
C            CHEBMO - Real
 
195
C                     Array of dimension at least (MAXP1,25), providing
 
196
C                     space for the Chebyshev moments needed within the
 
197
C                     cycles
 
198
C
 
199
C***REFERENCES  (NONE)
 
200
C***ROUTINES CALLED  QAGIE, QAWOE, QELG, R1MACH
 
201
C***REVISION HISTORY  (YYMMDD)
 
202
C   800101  DATE WRITTEN
 
203
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
204
C   890831  Modified array declarations.  (WRB)
 
205
C   891009  Removed unreferenced variable.  (WRB)
 
206
C   891009  REVISION DATE from Version 3.2
 
207
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
208
C***END PROLOGUE  QAWFE
 
209
C
 
210
      REAL A,ABSEPS,ABSERR,ALIST,BLIST,CHEBMO,CORREC,CYCLE,
 
211
     1  C1,C2,DL,DRL,ELIST,EP,EPS,EPSA,EPSABS,ERLST,
 
212
     2  ERRSUM,FACT,OMEGA,P,PI,P1,PSUM,RESEPS,RESULT,RES3LA,RLIST,RSLST
 
213
     3  ,R1MACH,UFLOW
 
214
      INTEGER IER,IERLST,INTEGR,IORD,KTMIN,L,LST,LIMIT,LL,MAXP1,
 
215
     1    NEV,NEVAL,NNLOG,NRES,NUMRL2
 
216
C
 
217
      DIMENSION ALIST(*),BLIST(*),CHEBMO(MAXP1,25),ELIST(*),
 
218
     1  ERLST(*),IERLST(*),IORD(*),NNLOG(*),PSUM(52),
 
219
     2  RES3LA(3),RLIST(*),RSLST(*)
 
220
C
 
221
      EXTERNAL F
 
222
C
 
223
C
 
224
C            THE DIMENSION OF  PSUM  IS DETERMINED BY THE VALUE OF
 
225
C            LIMEXP IN SUBROUTINE QELG (PSUM MUST BE
 
226
C            OF DIMENSION (LIMEXP+2) AT LEAST).
 
227
C
 
228
C           LIST OF MAJOR VARIABLES
 
229
C           -----------------------
 
230
C
 
231
C           C1, C2    - END POINTS OF SUBINTERVAL (OF LENGTH
 
232
C                       CYCLE)
 
233
C           CYCLE     - (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA)
 
234
C           PSUM      - VECTOR OF DIMENSION AT LEAST (LIMEXP+2)
 
235
C                       (SEE ROUTINE QELG)
 
236
C                       PSUM CONTAINS THE PART OF THE EPSILON
 
237
C                       TABLE WHICH IS STILL NEEDED FOR FURTHER
 
238
C                       COMPUTATIONS.
 
239
C                       EACH ELEMENT OF PSUM IS A PARTIAL SUM OF
 
240
C                       THE SERIES WHICH SHOULD SUM TO THE VALUE OF
 
241
C                       THE INTEGRAL.
 
242
C           ERRSUM    - SUM OF ERROR ESTIMATES OVER THE
 
243
C                       SUBINTERVALS, CALCULATED CUMULATIVELY
 
244
C           EPSA      - ABSOLUTE TOLERANCE REQUESTED OVER CURRENT
 
245
C                       SUBINTERVAL
 
246
C           CHEBMO    - ARRAY CONTAINING THE MODIFIED CHEBYSHEV
 
247
C                       MOMENTS (SEE ALSO ROUTINE QC25F)
 
248
C
 
249
      SAVE P, PI
 
250
      DATA P/0.9E+00/,PI/0.31415926535897932E+01/
 
251
C
 
252
C           TEST ON VALIDITY OF PARAMETERS
 
253
C           ------------------------------
 
254
C
 
255
C***FIRST EXECUTABLE STATEMENT  QAWFE
 
256
      RESULT = 0.0E+00
 
257
      ABSERR = 0.0E+00
 
258
      NEVAL = 0
 
259
      LST = 0
 
260
      IER = 0
 
261
      IF((INTEGR.NE.1.AND.INTEGR.NE.2).OR.EPSABS.LE.0.0E+00.OR.
 
262
     1  LIMLST.LT.3) IER = 6
 
263
      IF(IER.EQ.6) GO TO 999
 
264
      IF(OMEGA.NE.0.0E+00) GO TO 10
 
265
C
 
266
C           INTEGRATION BY QAGIE IF OMEGA IS ZERO
 
267
C           --------------------------------------
 
268
C
 
269
      IF(INTEGR.EQ.1) CALL QAGIE(F,A,1,EPSABS,0.0E+00,LIMIT,
 
270
     1  RESULT,ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST)
 
271
      RSLST(1) = RESULT
 
272
      ERLST(1) = ABSERR
 
273
      IERLST(1) = IER
 
274
      LST = 1
 
275
      GO TO 999
 
276
C
 
277
C           INITIALIZATIONS
 
278
C           ---------------
 
279
C
 
280
   10 L = ABS(OMEGA)
 
281
      DL = 2*L+1
 
282
      CYCLE = DL*PI/ABS(OMEGA)
 
283
      IER = 0
 
284
      KTMIN = 0
 
285
      NEVAL = 0
 
286
      NUMRL2 = 0
 
287
      NRES = 0
 
288
      C1 = A
 
289
      C2 = CYCLE+A
 
290
      P1 = 0.1E+01-P
 
291
      EPS = EPSABS
 
292
      UFLOW = R1MACH(1)
 
293
      IF(EPSABS.GT.UFLOW/P1) EPS = EPSABS*P1
 
294
      EP = EPS
 
295
      FACT = 0.1E+01
 
296
      CORREC = 0.0E+00
 
297
      ABSERR = 0.0E+00
 
298
      ERRSUM = 0.0E+00
 
299
C
 
300
C           MAIN DO-LOOP
 
301
C           ------------
 
302
C
 
303
      DO 50 LST = 1,LIMLST
 
304
C
 
305
C           INTEGRATE OVER CURRENT SUBINTERVAL.
 
306
C
 
307
        EPSA = EPS*FACT
 
308
        CALL QAWOE(F,C1,C2,OMEGA,INTEGR,EPSA,0.0E+00,LIMIT,LST,MAXP1,
 
309
     1  RSLST(LST),ERLST(LST),NEV,IERLST(LST),LAST,ALIST,BLIST,RLIST,
 
310
     2  ELIST,IORD,NNLOG,MOMCOM,CHEBMO)
 
311
        NEVAL = NEVAL+NEV
 
312
        FACT = FACT*P
 
313
        ERRSUM = ERRSUM+ERLST(LST)
 
314
        DRL = 0.5E+02*ABS(RSLST(LST))
 
315
C
 
316
C           TEST ON ACCURACY WITH PARTIAL SUM
 
317
C
 
318
        IF(ERRSUM+DRL.LE.EPSABS.AND.LST.GE.6) GO TO 80
 
319
        CORREC = MAX(CORREC,ERLST(LST))
 
320
        IF(IERLST(LST).NE.0) EPS = MAX(EP,CORREC*P1)
 
321
        IF(IERLST(LST).NE.0) IER = 7
 
322
        IF(IER.EQ.7.AND.(ERRSUM+DRL).LE.CORREC*0.1E+02.AND.
 
323
     1  LST.GT.5) GO TO 80
 
324
        NUMRL2 = NUMRL2+1
 
325
        IF(LST.GT.1) GO TO 20
 
326
        PSUM(1) = RSLST(1)
 
327
        GO TO 40
 
328
   20   PSUM(NUMRL2) = PSUM(LL)+RSLST(LST)
 
329
        IF(LST.EQ.2) GO TO 40
 
330
C
 
331
C           TEST ON MAXIMUM NUMBER OF SUBINTERVALS
 
332
C
 
333
        IF(LST.EQ.LIMLST) IER = 1
 
334
C
 
335
C           PERFORM NEW EXTRAPOLATION
 
336
C
 
337
        CALL QELG(NUMRL2,PSUM,RESEPS,ABSEPS,RES3LA,NRES)
 
338
C
 
339
C           TEST WHETHER EXTRAPOLATED RESULT IS INFLUENCED BY
 
340
C           ROUNDOFF
 
341
C
 
342
        KTMIN = KTMIN+1
 
343
        IF(KTMIN.GE.15.AND.ABSERR.LE.0.1E-02*(ERRSUM+DRL)) IER = 4
 
344
        IF(ABSEPS.GT.ABSERR.AND.LST.NE.3) GO TO 30
 
345
        ABSERR = ABSEPS
 
346
        RESULT = RESEPS
 
347
        KTMIN = 0
 
348
C
 
349
C           IF IER IS NOT 0, CHECK WHETHER DIRECT RESULT (PARTIAL
 
350
C           SUM) OR EXTRAPOLATED RESULT YIELDS THE BEST INTEGRAL
 
351
C           APPROXIMATION
 
352
C
 
353
        IF((ABSERR+0.1E+02*CORREC).LE.EPSABS.OR.
 
354
     1  (ABSERR.LE.EPSABS.AND.0.1E+02*CORREC.GE.EPSABS)) GO TO 60
 
355
   30   IF(IER.NE.0.AND.IER.NE.7) GO TO 60
 
356
   40   LL = NUMRL2
 
357
        C1 = C2
 
358
        C2 = C2+CYCLE
 
359
   50 CONTINUE
 
360
C
 
361
C         SET FINAL RESULT AND ERROR ESTIMATE
 
362
C         -----------------------------------
 
363
C
 
364
   60 ABSERR = ABSERR+0.1E+02*CORREC
 
365
      IF(IER.EQ.0) GO TO 999
 
366
      IF(RESULT.NE.0.0E+00.AND.PSUM(NUMRL2).NE.0.0E+00) GO TO 70
 
367
      IF(ABSERR.GT.ERRSUM) GO TO 80
 
368
      IF(PSUM(NUMRL2).EQ.0.0E+00) GO TO 999
 
369
   70 IF(ABSERR/ABS(RESULT).GT.(ERRSUM+DRL)/ABS(PSUM(NUMRL2)))
 
370
     1  GO TO 80
 
371
      IF(IER.GE.1.AND.IER.NE.7) ABSERR = ABSERR+DRL
 
372
      GO TO 999
 
373
   80 RESULT = PSUM(NUMRL2)
 
374
      ABSERR = ERRSUM+DRL
 
375
  999 RETURN
 
376
      END