← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dbesk.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/dbesk.f
 
1
*DECK DBESK
 
2
      SUBROUTINE DBESK (X, FNU, KODE, N, Y, NZ)
 
3
C***BEGIN PROLOGUE  DBESK
 
4
C***PURPOSE  Implement forward recursion on the three term recursion
 
5
C            relation for a sequence of non-negative order Bessel
 
6
C            functions K/SUB(FNU+I-1)/(X), or scaled Bessel functions
 
7
C            EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N for real, positive
 
8
C            X and non-negative orders FNU.
 
9
C***LIBRARY   SLATEC
 
10
C***CATEGORY  C10B3
 
11
C***TYPE      DOUBLE PRECISION (BESK-S, DBESK-D)
 
12
C***KEYWORDS  K BESSEL FUNCTION, SPECIAL FUNCTIONS
 
13
C***AUTHOR  Amos, D. E., (SNLA)
 
14
C***DESCRIPTION
 
15
C
 
16
C     Abstract  **** a double precision routine ****
 
17
C         DBESK implements forward recursion on the three term
 
18
C         recursion relation for a sequence of non-negative order Bessel
 
19
C         functions K/sub(FNU+I-1)/(X), or scaled Bessel functions
 
20
C         EXP(X)*K/sub(FNU+I-1)/(X), I=1,..,N for real X .GT. 0.0D0 and
 
21
C         non-negative orders FNU.  If FNU .LT. NULIM, orders FNU and
 
22
C         FNU+1 are obtained from DBSKNU to start the recursion.  If
 
23
C         FNU .GE. NULIM, the uniform asymptotic expansion is used for
 
24
C         orders FNU and FNU+1 to start the recursion.  NULIM is 35 or
 
25
C         70 depending on whether N=1 or N .GE. 2.  Under and overflow
 
26
C         tests are made on the leading term of the asymptotic expansion
 
27
C         before any extensive computation is done.
 
28
C
 
29
C         The maximum number of significant digits obtainable
 
30
C         is the smaller of 14 and the number of digits carried in
 
31
C         double precision arithmetic.
 
32
C
 
33
C     Description of Arguments
 
34
C
 
35
C         Input      X,FNU are double precision
 
36
C           X      - X .GT. 0.0D0
 
37
C           FNU    - order of the initial K function, FNU .GE. 0.0D0
 
38
C           KODE   - a parameter to indicate the scaling option
 
39
C                    KODE=1 returns Y(I)=       K/sub(FNU+I-1)/(X),
 
40
C                                        I=1,...,N
 
41
C                    KODE=2 returns Y(I)=EXP(X)*K/sub(FNU+I-1)/(X),
 
42
C                                        I=1,...,N
 
43
C           N      - number of members in the sequence, N .GE. 1
 
44
C
 
45
C         Output     Y is double precision
 
46
C           Y      - a vector whose first N components contain values
 
47
C                    for the sequence
 
48
C                    Y(I)=       k/sub(FNU+I-1)/(X), I=1,...,N  or
 
49
C                    Y(I)=EXP(X)*K/sub(FNU+I-1)/(X), I=1,...,N
 
50
C                    depending on KODE
 
51
C           NZ     - number of components of Y set to zero due to
 
52
C                    underflow with KODE=1,
 
53
C                    NZ=0   , normal return, computation completed
 
54
C                    NZ .NE. 0, first NZ components of Y set to zero
 
55
C                             due to underflow, Y(I)=0.0D0, I=1,...,NZ
 
56
C
 
57
C     Error Conditions
 
58
C         Improper input arguments - a fatal error
 
59
C         Overflow - a fatal error
 
60
C         Underflow with KODE=1 -  a non-fatal error (NZ .NE. 0)
 
61
C
 
62
C***REFERENCES  F. W. J. Olver, Tables of Bessel Functions of Moderate
 
63
C                 or Large Orders, NPL Mathematical Tables 6, Her
 
64
C                 Majesty's Stationery Office, London, 1962.
 
65
C               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  D1MACH, DASYIK, DBESK0, DBESK1, DBSK0E, DBSK1E,
 
69
C                    DBSKNU, I1MACH, XERMSG
 
70
C***REVISION HISTORY  (YYMMDD)
 
71
C   790201  DATE WRITTEN
 
72
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
73
C   890911  Removed unnecessary intrinsics.  (WRB)
 
74
C   890911  REVISION DATE from Version 3.2
 
75
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
76
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
77
C   920501  Reformatted the REFERENCES section.  (WRB)
 
78
C***END PROLOGUE  DBESK
 
79
C
 
80
      INTEGER I, J, K, KODE, MZ, N, NB, ND, NN, NUD, NULIM, NZ
 
81
      INTEGER I1MACH
 
82
      DOUBLE PRECISION CN,DNU,ELIM,ETX,FLGIK,FN,FNN,FNU,GLN,GNU,RTZ,
 
83
     1 S, S1, S2, T, TM, TRX, W, X, XLIM, Y, ZN
 
84
      DOUBLE PRECISION DBESK0, DBESK1, DBSK1E, DBSK0E, D1MACH
 
85
      DIMENSION W(2), NULIM(2), Y(*)
 
86
      SAVE NULIM
 
87
      DATA NULIM(1),NULIM(2) / 35 , 70 /
 
88
C***FIRST EXECUTABLE STATEMENT  DBESK
 
89
      NN = -I1MACH(15)
 
90
      ELIM = 2.303D0*(NN*D1MACH(5)-3.0D0)
 
91
      XLIM = D1MACH(1)*1.0D+3
 
92
      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 280
 
93
      IF (FNU.LT.0.0D0) GO TO 290
 
94
      IF (X.LE.0.0D0) GO TO 300
 
95
      IF (X.LT.XLIM) GO TO 320
 
96
      IF (N.LT.1) GO TO 310
 
97
      ETX = KODE - 1
 
98
C
 
99
C     ND IS A DUMMY VARIABLE FOR N
 
100
C     GNU IS A DUMMY VARIABLE FOR FNU
 
101
C     NZ = NUMBER OF UNDERFLOWS ON KODE=1
 
102
C
 
103
      ND = N
 
104
      NZ = 0
 
105
      NUD = INT(FNU)
 
106
      DNU = FNU - NUD
 
107
      GNU = FNU
 
108
      NN = MIN(2,ND)
 
109
      FN = FNU + N - 1
 
110
      FNN = FN
 
111
      IF (FN.LT.2.0D0) GO TO 150
 
112
C
 
113
C     OVERFLOW TEST  (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
 
114
C     FOR THE LAST ORDER, FNU+N-1.GE.NULIM
 
115
C
 
116
      ZN = X/FN
 
117
      IF (ZN.EQ.0.0D0) GO TO 320
 
118
      RTZ = SQRT(1.0D0+ZN*ZN)
 
119
      GLN = LOG((1.0D0+RTZ)/ZN)
 
120
      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)
 
121
      CN = -FN*(T-GLN)
 
122
      IF (CN.GT.ELIM) GO TO 320
 
123
      IF (NUD.LT.NULIM(NN)) GO TO 30
 
124
      IF (NN.EQ.1) GO TO 20
 
125
   10 CONTINUE
 
126
C
 
127
C     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
 
128
C     FOR THE FIRST ORDER, FNU.GE.NULIM
 
129
C
 
130
      FN = GNU
 
131
      ZN = X/FN
 
132
      RTZ = SQRT(1.0D0+ZN*ZN)
 
133
      GLN = LOG((1.0D0+RTZ)/ZN)
 
134
      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)
 
135
      CN = -FN*(T-GLN)
 
136
   20 CONTINUE
 
137
      IF (CN.LT.-ELIM) GO TO 230
 
138
C
 
139
C     ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
 
140
C
 
141
      FLGIK = -1.0D0
 
142
      CALL DASYIK(X,GNU,KODE,FLGIK,RTZ,CN,NN,Y)
 
143
      IF (NN.EQ.1) GO TO 240
 
144
      TRX = 2.0D0/X
 
145
      TM = (GNU+GNU+2.0D0)/X
 
146
      GO TO 130
 
147
C
 
148
   30 CONTINUE
 
149
      IF (KODE.EQ.2) GO TO 40
 
150
C
 
151
C     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X)
 
152
C     FOR ORDER DNU
 
153
C
 
154
      IF (X.GT.ELIM) GO TO 230
 
155
   40 CONTINUE
 
156
      IF (DNU.NE.0.0D0) GO TO 80
 
157
      IF (KODE.EQ.2) GO TO 50
 
158
      S1 = DBESK0(X)
 
159
      GO TO 60
 
160
   50 S1 = DBSK0E(X)
 
161
   60 CONTINUE
 
162
      IF (NUD.EQ.0 .AND. ND.EQ.1) GO TO 120
 
163
      IF (KODE.EQ.2) GO TO 70
 
164
      S2 = DBESK1(X)
 
165
      GO TO 90
 
166
   70 S2 = DBSK1E(X)
 
167
      GO TO 90
 
168
   80 CONTINUE
 
169
      NB = 2
 
170
      IF (NUD.EQ.0 .AND. ND.EQ.1) NB = 1
 
171
      CALL DBSKNU(X, DNU, KODE, NB, W, NZ)
 
172
      S1 = W(1)
 
173
      IF (NB.EQ.1) GO TO 120
 
174
      S2 = W(2)
 
175
   90 CONTINUE
 
176
      TRX = 2.0D0/X
 
177
      TM = (DNU+DNU+2.0D0)/X
 
178
C     FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
 
179
      IF (ND.EQ.1) NUD = NUD - 1
 
180
      IF (NUD.GT.0) GO TO 100
 
181
      IF (ND.GT.1) GO TO 120
 
182
      S1 = S2
 
183
      GO TO 120
 
184
  100 CONTINUE
 
185
      DO 110 I=1,NUD
 
186
        S = S2
 
187
        S2 = TM*S2 + S1
 
188
        S1 = S
 
189
        TM = TM + TRX
 
190
  110 CONTINUE
 
191
      IF (ND.EQ.1) S1 = S2
 
192
  120 CONTINUE
 
193
      Y(1) = S1
 
194
      IF (ND.EQ.1) GO TO 240
 
195
      Y(2) = S2
 
196
  130 CONTINUE
 
197
      IF (ND.EQ.2) GO TO 240
 
198
C     FORWARD RECUR FROM FNU+2 TO FNU+N-1
 
199
      DO 140 I=3,ND
 
200
        Y(I) = TM*Y(I-1) + Y(I-2)
 
201
        TM = TM + TRX
 
202
  140 CONTINUE
 
203
      GO TO 240
 
204
C
 
205
  150 CONTINUE
 
206
C     UNDERFLOW TEST FOR KODE=1
 
207
      IF (KODE.EQ.2) GO TO 160
 
208
      IF (X.GT.ELIM) GO TO 230
 
209
  160 CONTINUE
 
210
C     OVERFLOW TEST
 
211
      IF (FN.LE.1.0D0) GO TO 170
 
212
      IF (-FN*(LOG(X)-0.693D0).GT.ELIM) GO TO 320
 
213
  170 CONTINUE
 
214
      IF (DNU.EQ.0.0D0) GO TO 180
 
215
      CALL DBSKNU(X, FNU, KODE, ND, Y, MZ)
 
216
      GO TO 240
 
217
  180 CONTINUE
 
218
      J = NUD
 
219
      IF (J.EQ.1) GO TO 210
 
220
      J = J + 1
 
221
      IF (KODE.EQ.2) GO TO 190
 
222
      Y(J) = DBESK0(X)
 
223
      GO TO 200
 
224
  190 Y(J) = DBSK0E(X)
 
225
  200 IF (ND.EQ.1) GO TO 240
 
226
      J = J + 1
 
227
  210 IF (KODE.EQ.2) GO TO 220
 
228
      Y(J) = DBESK1(X)
 
229
      GO TO 240
 
230
  220 Y(J) = DBSK1E(X)
 
231
      GO TO 240
 
232
C
 
233
C     UPDATE PARAMETERS ON UNDERFLOW
 
234
C
 
235
  230 CONTINUE
 
236
      NUD = NUD + 1
 
237
      ND = ND - 1
 
238
      IF (ND.EQ.0) GO TO 240
 
239
      NN = MIN(2,ND)
 
240
      GNU = GNU + 1.0D0
 
241
      IF (FNN.LT.2.0D0) GO TO 230
 
242
      IF (NUD.LT.NULIM(NN)) GO TO 230
 
243
      GO TO 10
 
244
  240 CONTINUE
 
245
      NZ = N - ND
 
246
      IF (NZ.EQ.0) RETURN
 
247
      IF (ND.EQ.0) GO TO 260
 
248
      DO 250 I=1,ND
 
249
        J = N - I + 1
 
250
        K = ND - I + 1
 
251
        Y(J) = Y(K)
 
252
  250 CONTINUE
 
253
  260 CONTINUE
 
254
      DO 270 I=1,NZ
 
255
        Y(I) = 0.0D0
 
256
  270 CONTINUE
 
257
      RETURN
 
258
C
 
259
C
 
260
C
 
261
  280 CONTINUE
 
262
      CALL XERMSG ('SLATEC', 'DBESK',
 
263
     +   'SCALING OPTION, KODE, NOT 1 OR 2', 2, 1)
 
264
      RETURN
 
265
  290 CONTINUE
 
266
      CALL XERMSG ('SLATEC', 'DBESK', 'ORDER, FNU, LESS THAN ZERO', 2,
 
267
     +   1)
 
268
      RETURN
 
269
  300 CONTINUE
 
270
      CALL XERMSG ('SLATEC', 'DBESK', 'X LESS THAN OR EQUAL TO ZERO',
 
271
     +   2, 1)
 
272
      RETURN
 
273
  310 CONTINUE
 
274
      CALL XERMSG ('SLATEC', 'DBESK', 'N LESS THAN ONE', 2, 1)
 
275
      RETURN
 
276
  320 CONTINUE
 
277
      CALL XERMSG ('SLATEC', 'DBESK',
 
278
     +   'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL', 6, 1)
 
279
      RETURN
 
280
      END