← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/besknu.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/besknu.f
 
1
*DECK BESKNU
 
2
      SUBROUTINE BESKNU (X, FNU, KODE, N, Y, NZ)
 
3
C***BEGIN PROLOGUE  BESKNU
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to BESK
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      SINGLE PRECISION (BESKNU-S, DBSKNU-D)
 
8
C***AUTHOR  Amos, D. E., (SNLA)
 
9
C***DESCRIPTION
 
10
C
 
11
C     Abstract
 
12
C         BESKNU computes N member sequences of K Bessel functions
 
13
C         K/SUB(FNU+I-1)/(X), I=1,N for non-negative orders FNU and
 
14
C         positive X. Equations of the references are implemented on
 
15
C         small orders DNU for K/SUB(DNU)/(X) and K/SUB(DNU+1)/(X).
 
16
C         Forward recursion with the three term recursion relation
 
17
C         generates higher orders FNU+I-1, I=1,...,N. The parameter
 
18
C         KODE permits K/SUB(FNU+I-1)/(X) values or scaled values
 
19
C         EXP(X)*K/SUB(FNU+I-1)/(X), I=1,N to be returned.
 
20
C
 
21
C         To start the recursion FNU is normalized to the interval
 
22
C         -0.5.LE.DNU.LT.0.5. A special form of the power series is
 
23
C         implemented on 0.LT.X.LE.X1 while the Miller algorithm for the
 
24
C         K Bessel function in terms of the confluent hypergeometric
 
25
C         function U(FNU+0.5,2*FNU+1,X) is implemented on X1.LT.X.LE.X2.
 
26
C         For X.GT.X2, the asymptotic expansion for large X is used.
 
27
C         When FNU is a half odd integer, a special formula for
 
28
C         DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion.
 
29
C
 
30
C         BESKNU assumes that a significant digit SINH(X) function is
 
31
C         available.
 
32
C
 
33
C     Description of Arguments
 
34
C
 
35
C         Input
 
36
C           X      - X.GT.0.0E0
 
37
C           FNU    - Order of initial K function, FNU.GE.0.0E0
 
38
C           N      - Number of members of the sequence, N.GE.1
 
39
C           KODE   - A parameter to indicate the scaling option
 
40
C                    KODE= 1  returns
 
41
C                             Y(I)=       K/SUB(FNU+I-1)/(X)
 
42
C                                  I=1,...,N
 
43
C                        = 2  returns
 
44
C                             Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X)
 
45
C                                  I=1,...,N
 
46
C
 
47
C         Output
 
48
C           Y      - A vector whose first N components contain values
 
49
C                    for the sequence
 
50
C                    Y(I)=       K/SUB(FNU+I-1)/(X), I=1,...,N or
 
51
C                    Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N
 
52
C                    depending on KODE
 
53
C           NZ     - Number of components set to zero due to
 
54
C                    underflow,
 
55
C                    NZ= 0   , Normal return
 
56
C                    NZ.NE.0 , First NZ components of Y set to zero
 
57
C                              due to underflow, Y(I)=0.0E0,I=1,...,NZ
 
58
C
 
59
C     Error Conditions
 
60
C         Improper input arguments - a fatal error
 
61
C         Overflow - a fatal error
 
62
C         Underflow with KODE=1 - a non-fatal error (NZ.NE.0)
 
63
C
 
64
C***SEE ALSO  BESK
 
65
C***REFERENCES  N. M. Temme, On the numerical evaluation of the modified
 
66
C                 Bessel function of the third kind, Journal of
 
67
C                 Computational Physics 19, (1975), pp. 324-337.
 
68
C***ROUTINES CALLED  GAMMA, I1MACH, R1MACH, XERMSG
 
69
C***REVISION HISTORY  (YYMMDD)
 
70
C   790201  DATE WRITTEN
 
71
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
72
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
73
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
74
C   900326  Removed duplicate information from DESCRIPTION section.
 
75
C           (WRB)
 
76
C   900328  Added TYPE section.  (WRB)
 
77
C   900727  Added EXTERNAL statement.  (WRB)
 
78
C   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
 
79
C   920501  Reformatted the REFERENCES section.  (WRB)
 
80
C***END PROLOGUE  BESKNU
 
81
C
 
82
      INTEGER I, IFLAG, INU, J, K, KK, KODE, KODED, N, NN, NZ
 
83
      INTEGER I1MACH
 
84
      REAL A, AK, A1, A2, B, BK, CC, CK, COEF, CX, DK, DNU, DNU2, ELIM,
 
85
     1 ETEST, EX, F, FC, FHS, FK, FKS, FLRX, FMU, FNU, G1, G2, P, PI,
 
86
     2 PT, P1, P2, Q, RTHPI, RX, S, SMU, SQK, ST, S1, S2, TM, TOL, T1,
 
87
     3 T2, X, X1, X2, Y
 
88
      REAL GAMMA, R1MACH
 
89
      DIMENSION A(160), B(160), Y(*), CC(8)
 
90
      EXTERNAL GAMMA
 
91
      SAVE X1, X2, PI, RTHPI, CC
 
92
      DATA X1, X2 / 2.0E0, 17.0E0 /
 
93
      DATA PI,RTHPI        / 3.14159265358979E+00, 1.25331413731550E+00/
 
94
      DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)
 
95
     1                     / 5.77215664901533E-01,-4.20026350340952E-02,
 
96
     2-4.21977345555443E-02, 7.21894324666300E-03,-2.15241674114900E-04,
 
97
     3-2.01348547807000E-05, 1.13302723200000E-06, 6.11609500000000E-09/
 
98
C***FIRST EXECUTABLE STATEMENT  BESKNU
 
99
      KK = -I1MACH(12)
 
100
      ELIM = 2.303E0*(KK*R1MACH(5)-3.0E0)
 
101
      AK = R1MACH(3)
 
102
      TOL = MAX(AK,1.0E-15)
 
103
      IF (X.LE.0.0E0) GO TO 350
 
104
      IF (FNU.LT.0.0E0) GO TO 360
 
105
      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 370
 
106
      IF (N.LT.1) GO TO 380
 
107
      NZ = 0
 
108
      IFLAG = 0
 
109
      KODED = KODE
 
110
      RX = 2.0E0/X
 
111
      INU = INT(FNU+0.5E0)
 
112
      DNU = FNU - INU
 
113
      IF (ABS(DNU).EQ.0.5E0) GO TO 120
 
114
      DNU2 = 0.0E0
 
115
      IF (ABS(DNU).LT.TOL) GO TO 10
 
116
      DNU2 = DNU*DNU
 
117
   10 CONTINUE
 
118
      IF (X.GT.X1) GO TO 120
 
119
C
 
120
C     SERIES FOR X.LE.X1
 
121
C
 
122
      A1 = 1.0E0 - DNU
 
123
      A2 = 1.0E0 + DNU
 
124
      T1 = 1.0E0/GAMMA(A1)
 
125
      T2 = 1.0E0/GAMMA(A2)
 
126
      IF (ABS(DNU).GT.0.1E0) GO TO 40
 
127
C     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
 
128
      S = CC(1)
 
129
      AK = 1.0E0
 
130
      DO 20 K=2,8
 
131
        AK = AK*DNU2
 
132
        TM = CC(K)*AK
 
133
        S = S + TM
 
134
        IF (ABS(TM).LT.TOL) GO TO 30
 
135
   20 CONTINUE
 
136
   30 G1 = -S
 
137
      GO TO 50
 
138
   40 CONTINUE
 
139
      G1 = (T1-T2)/(DNU+DNU)
 
140
   50 CONTINUE
 
141
      G2 = (T1+T2)*0.5E0
 
142
      SMU = 1.0E0
 
143
      FC = 1.0E0
 
144
      FLRX = LOG(RX)
 
145
      FMU = DNU*FLRX
 
146
      IF (DNU.EQ.0.0E0) GO TO 60
 
147
      FC = DNU*PI
 
148
      FC = FC/SIN(FC)
 
149
      IF (FMU.NE.0.0E0) SMU = SINH(FMU)/FMU
 
150
   60 CONTINUE
 
151
      F = FC*(G1*COSH(FMU)+G2*FLRX*SMU)
 
152
      FC = EXP(FMU)
 
153
      P = 0.5E0*FC/T2
 
154
      Q = 0.5E0/(FC*T1)
 
155
      AK = 1.0E0
 
156
      CK = 1.0E0
 
157
      BK = 1.0E0
 
158
      S1 = F
 
159
      S2 = P
 
160
      IF (INU.GT.0 .OR. N.GT.1) GO TO 90
 
161
      IF (X.LT.TOL) GO TO 80
 
162
      CX = X*X*0.25E0
 
163
   70 CONTINUE
 
164
      F = (AK*F+P+Q)/(BK-DNU2)
 
165
      P = P/(AK-DNU)
 
166
      Q = Q/(AK+DNU)
 
167
      CK = CK*CX/AK
 
168
      T1 = CK*F
 
169
      S1 = S1 + T1
 
170
      BK = BK + AK + AK + 1.0E0
 
171
      AK = AK + 1.0E0
 
172
      S = ABS(T1)/(1.0E0+ABS(S1))
 
173
      IF (S.GT.TOL) GO TO 70
 
174
   80 CONTINUE
 
175
      Y(1) = S1
 
176
      IF (KODED.EQ.1) RETURN
 
177
      Y(1) = S1*EXP(X)
 
178
      RETURN
 
179
   90 CONTINUE
 
180
      IF (X.LT.TOL) GO TO 110
 
181
      CX = X*X*0.25E0
 
182
  100 CONTINUE
 
183
      F = (AK*F+P+Q)/(BK-DNU2)
 
184
      P = P/(AK-DNU)
 
185
      Q = Q/(AK+DNU)
 
186
      CK = CK*CX/AK
 
187
      T1 = CK*F
 
188
      S1 = S1 + T1
 
189
      T2 = CK*(P-AK*F)
 
190
      S2 = S2 + T2
 
191
      BK = BK + AK + AK + 1.0E0
 
192
      AK = AK + 1.0E0
 
193
      S = ABS(T1)/(1.0E0+ABS(S1)) + ABS(T2)/(1.0E0+ABS(S2))
 
194
      IF (S.GT.TOL) GO TO 100
 
195
  110 CONTINUE
 
196
      S2 = S2*RX
 
197
      IF (KODED.EQ.1) GO TO 170
 
198
      F = EXP(X)
 
199
      S1 = S1*F
 
200
      S2 = S2*F
 
201
      GO TO 170
 
202
  120 CONTINUE
 
203
      COEF = RTHPI/SQRT(X)
 
204
      IF (KODED.EQ.2) GO TO 130
 
205
      IF (X.GT.ELIM) GO TO 330
 
206
      COEF = COEF*EXP(-X)
 
207
  130 CONTINUE
 
208
      IF (ABS(DNU).EQ.0.5E0) GO TO 340
 
209
      IF (X.GT.X2) GO TO 280
 
210
C
 
211
C     MILLER ALGORITHM FOR X1.LT.X.LE.X2
 
212
C
 
213
      ETEST = COS(PI*DNU)/(PI*X*TOL)
 
214
      FKS = 1.0E0
 
215
      FHS = 0.25E0
 
216
      FK = 0.0E0
 
217
      CK = X + X + 2.0E0
 
218
      P1 = 0.0E0
 
219
      P2 = 1.0E0
 
220
      K = 0
 
221
  140 CONTINUE
 
222
      K = K + 1
 
223
      FK = FK + 1.0E0
 
224
      AK = (FHS-DNU2)/(FKS+FK)
 
225
      BK = CK/(FK+1.0E0)
 
226
      PT = P2
 
227
      P2 = BK*P2 - AK*P1
 
228
      P1 = PT
 
229
      A(K) = AK
 
230
      B(K) = BK
 
231
      CK = CK + 2.0E0
 
232
      FKS = FKS + FK + FK + 1.0E0
 
233
      FHS = FHS + FK + FK
 
234
      IF (ETEST.GT.FK*P1) GO TO 140
 
235
      KK = K
 
236
      S = 1.0E0
 
237
      P1 = 0.0E0
 
238
      P2 = 1.0E0
 
239
      DO 150 I=1,K
 
240
        PT = P2
 
241
        P2 = (B(KK)*P2-P1)/A(KK)
 
242
        P1 = PT
 
243
        S = S + P2
 
244
        KK = KK - 1
 
245
  150 CONTINUE
 
246
      S1 = COEF*(P2/S)
 
247
      IF (INU.GT.0 .OR. N.GT.1) GO TO 160
 
248
      GO TO 200
 
249
  160 CONTINUE
 
250
      S2 = S1*(X+DNU+0.5E0-P1/P2)/X
 
251
C
 
252
C     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
 
253
C
 
254
  170 CONTINUE
 
255
      CK = (DNU+DNU+2.0E0)/X
 
256
      IF (N.EQ.1) INU = INU - 1
 
257
      IF (INU.GT.0) GO TO 180
 
258
      IF (N.GT.1) GO TO 200
 
259
      S1 = S2
 
260
      GO TO 200
 
261
  180 CONTINUE
 
262
      DO 190 I=1,INU
 
263
        ST = S2
 
264
        S2 = CK*S2 + S1
 
265
        S1 = ST
 
266
        CK = CK + RX
 
267
  190 CONTINUE
 
268
      IF (N.EQ.1) S1 = S2
 
269
  200 CONTINUE
 
270
      IF (IFLAG.EQ.1) GO TO 220
 
271
      Y(1) = S1
 
272
      IF (N.EQ.1) RETURN
 
273
      Y(2) = S2
 
274
      IF (N.EQ.2) RETURN
 
275
      DO 210 I=3,N
 
276
        Y(I) = CK*Y(I-1) + Y(I-2)
 
277
        CK = CK + RX
 
278
  210 CONTINUE
 
279
      RETURN
 
280
C     IFLAG=1 CASES
 
281
  220 CONTINUE
 
282
      S = -X + LOG(S1)
 
283
      Y(1) = 0.0E0
 
284
      NZ = 1
 
285
      IF (S.LT.-ELIM) GO TO 230
 
286
      Y(1) = EXP(S)
 
287
      NZ = 0
 
288
  230 CONTINUE
 
289
      IF (N.EQ.1) RETURN
 
290
      S = -X + LOG(S2)
 
291
      Y(2) = 0.0E0
 
292
      NZ = NZ + 1
 
293
      IF (S.LT.-ELIM) GO TO 240
 
294
      NZ = NZ - 1
 
295
      Y(2) = EXP(S)
 
296
  240 CONTINUE
 
297
      IF (N.EQ.2) RETURN
 
298
      KK = 2
 
299
      IF (NZ.LT.2) GO TO 260
 
300
      DO 250 I=3,N
 
301
        KK = I
 
302
        ST = S2
 
303
        S2 = CK*S2 + S1
 
304
        S1 = ST
 
305
        CK = CK + RX
 
306
        S = -X + LOG(S2)
 
307
        NZ = NZ + 1
 
308
        Y(I) = 0.0E0
 
309
        IF (S.LT.-ELIM) GO TO 250
 
310
        Y(I) = EXP(S)
 
311
        NZ = NZ - 1
 
312
        GO TO 260
 
313
  250 CONTINUE
 
314
      RETURN
 
315
  260 CONTINUE
 
316
      IF (KK.EQ.N) RETURN
 
317
      S2 = S2*CK + S1
 
318
      CK = CK + RX
 
319
      KK = KK + 1
 
320
      Y(KK) = EXP(-X+LOG(S2))
 
321
      IF (KK.EQ.N) RETURN
 
322
      KK = KK + 1
 
323
      DO 270 I=KK,N
 
324
        Y(I) = CK*Y(I-1) + Y(I-2)
 
325
        CK = CK + RX
 
326
  270 CONTINUE
 
327
      RETURN
 
328
C
 
329
C     ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2
 
330
C
 
331
C     IFLAG=0 MEANS NO UNDERFLOW OCCURRED
 
332
C     IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
 
333
C     KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
 
334
C     RECURSION
 
335
  280 CONTINUE
 
336
      NN = 2
 
337
      IF (INU.EQ.0 .AND. N.EQ.1) NN = 1
 
338
      DNU2 = DNU + DNU
 
339
      FMU = 0.0E0
 
340
      IF (ABS(DNU2).LT.TOL) GO TO 290
 
341
      FMU = DNU2*DNU2
 
342
  290 CONTINUE
 
343
      EX = X*8.0E0
 
344
      S2 = 0.0E0
 
345
      DO 320 K=1,NN
 
346
        S1 = S2
 
347
        S = 1.0E0
 
348
        AK = 0.0E0
 
349
        CK = 1.0E0
 
350
        SQK = 1.0E0
 
351
        DK = EX
 
352
        DO 300 J=1,30
 
353
          CK = CK*(FMU-SQK)/DK
 
354
          S = S + CK
 
355
          DK = DK + EX
 
356
          AK = AK + 8.0E0
 
357
          SQK = SQK + AK
 
358
          IF (ABS(CK).LT.TOL) GO TO 310
 
359
  300   CONTINUE
 
360
  310   S2 = S*COEF
 
361
        FMU = FMU + 8.0E0*DNU + 4.0E0
 
362
  320 CONTINUE
 
363
      IF (NN.GT.1) GO TO 170
 
364
      S1 = S2
 
365
      GO TO 200
 
366
  330 CONTINUE
 
367
      KODED = 2
 
368
      IFLAG = 1
 
369
      GO TO 120
 
370
C
 
371
C     FNU=HALF ODD INTEGER CASE
 
372
C
 
373
  340 CONTINUE
 
374
      S1 = COEF
 
375
      S2 = COEF
 
376
      GO TO 170
 
377
C
 
378
C
 
379
  350 CALL XERMSG ('SLATEC', 'BESKNU', 'X NOT GREATER THAN ZERO', 2, 1)
 
380
      RETURN
 
381
  360 CALL XERMSG ('SLATEC', 'BESKNU', 'FNU NOT ZERO OR POSITIVE', 2,
 
382
     +   1)
 
383
      RETURN
 
384
  370 CALL XERMSG ('SLATEC', 'BESKNU', 'KODE NOT 1 OR 2', 2, 1)
 
385
      RETURN
 
386
  380 CALL XERMSG ('SLATEC', 'BESKNU', 'N NOT GREATER THAN 0', 2, 1)
 
387
      RETURN
 
388
      END