← Back to branch summary

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

« back to all changes in this revision

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