← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/cunk1.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/cunk1.f
 
1
*DECK CUNK1
 
2
      SUBROUTINE CUNK1 (Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  CUNK1
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to CBESK
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      ALL (CUNK1-A, ZUNK1-A)
 
8
C***AUTHOR  Amos, D. E., (SNL)
 
9
C***DESCRIPTION
 
10
C
 
11
C     CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
 
12
C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
 
13
C     UNIFORM ASYMPTOTIC EXPANSION.
 
14
C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
 
15
C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
 
16
C
 
17
C***SEE ALSO  CBESK
 
18
C***ROUTINES CALLED  CS1S2, CUCHK, CUNIK, R1MACH
 
19
C***REVISION HISTORY  (YYMMDD)
 
20
C   830501  DATE WRITTEN
 
21
C   910415  Prologue converted to Version 4.0 format.  (BAB)
 
22
C***END PROLOGUE  CUNK1
 
23
      COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS,
 
24
     * CWRK, CY, CZERO, C1, C2, PHI,  RZ, SUM,  S1, S2, Y, Z,
 
25
     * ZETA1,  ZETA2,  ZR, PHID, ZETA1D, ZETA2D, SUMD
 
26
      REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM,
 
27
     * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH
 
28
      INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
 
29
     * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC, M
 
30
      DIMENSION BRY(3), INIT(2), Y(N), SUM(2), PHI(2), ZETA1(2),
 
31
     * ZETA2(2), CY(2), CWRK(16,3), CSS(3), CSR(3)
 
32
      DATA CZERO, CONE / (0.0E0,0.0E0) , (1.0E0,0.0E0) /
 
33
      DATA PI / 3.14159265358979324E0 /
 
34
C***FIRST EXECUTABLE STATEMENT  CUNK1
 
35
      KDFLG = 1
 
36
      NZ = 0
 
37
C-----------------------------------------------------------------------
 
38
C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
 
39
C     THE UNDERFLOW LIMIT
 
40
C-----------------------------------------------------------------------
 
41
      CSCL = CMPLX(1.0E0/TOL,0.0E0)
 
42
      CRSC = CMPLX(TOL,0.0E0)
 
43
      CSS(1) = CSCL
 
44
      CSS(2) = CONE
 
45
      CSS(3) = CRSC
 
46
      CSR(1) = CRSC
 
47
      CSR(2) = CONE
 
48
      CSR(3) = CSCL
 
49
      BRY(1) = 1.0E+3*R1MACH(1)/TOL
 
50
      BRY(2) = 1.0E0/BRY(1)
 
51
      BRY(3) = R1MACH(2)
 
52
      X = REAL(Z)
 
53
      ZR = Z
 
54
      IF (X.LT.0.0E0) ZR = -Z
 
55
      J=2
 
56
      DO 70 I=1,N
 
57
C-----------------------------------------------------------------------
 
58
C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
 
59
C-----------------------------------------------------------------------
 
60
        J = 3 - J
 
61
        FN = FNU + (I-1)
 
62
        INIT(J) = 0
 
63
        CALL CUNIK(ZR, FN, 2, 0, TOL, INIT(J), PHI(J), ZETA1(J),
 
64
     *   ZETA2(J), SUM(J), CWRK(1,J))
 
65
        IF (KODE.EQ.1) GO TO 20
 
66
        CFN = CMPLX(FN,0.0E0)
 
67
        S1 = ZETA1(J) - CFN*(CFN/(ZR+ZETA2(J)))
 
68
        GO TO 30
 
69
   20   CONTINUE
 
70
        S1 = ZETA1(J) - ZETA2(J)
 
71
   30   CONTINUE
 
72
C-----------------------------------------------------------------------
 
73
C     TEST FOR UNDERFLOW AND OVERFLOW
 
74
C-----------------------------------------------------------------------
 
75
        RS1 = REAL(S1)
 
76
        IF (ABS(RS1).GT.ELIM) GO TO 60
 
77
        IF (KDFLG.EQ.1) KFLAG = 2
 
78
        IF (ABS(RS1).LT.ALIM) GO TO 40
 
79
C-----------------------------------------------------------------------
 
80
C     REFINE  TEST AND SCALE
 
81
C-----------------------------------------------------------------------
 
82
        APHI = ABS(PHI(J))
 
83
        RS1 = RS1 + ALOG(APHI)
 
84
        IF (ABS(RS1).GT.ELIM) GO TO 60
 
85
        IF (KDFLG.EQ.1) KFLAG = 1
 
86
        IF (RS1.LT.0.0E0) GO TO 40
 
87
        IF (KDFLG.EQ.1) KFLAG = 3
 
88
   40   CONTINUE
 
89
C-----------------------------------------------------------------------
 
90
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
 
91
C     EXPONENT EXTREMES
 
92
C-----------------------------------------------------------------------
 
93
        S2 = PHI(J)*SUM(J)
 
94
        C2R = REAL(S1)
 
95
        C2I = AIMAG(S1)
 
96
        C2M = EXP(C2R)*REAL(CSS(KFLAG))
 
97
        S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
 
98
        S2 = S2*S1
 
99
        IF (KFLAG.NE.1) GO TO 50
 
100
        CALL CUCHK(S2, NW, BRY(1), TOL)
 
101
        IF (NW.NE.0) GO TO 60
 
102
   50   CONTINUE
 
103
        CY(KDFLG) = S2
 
104
        Y(I) = S2*CSR(KFLAG)
 
105
        IF (KDFLG.EQ.2) GO TO 75
 
106
        KDFLG = 2
 
107
        GO TO 70
 
108
   60   CONTINUE
 
109
        IF (RS1.GT.0.0E0) GO TO 290
 
110
C-----------------------------------------------------------------------
 
111
C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 
112
C-----------------------------------------------------------------------
 
113
        IF (X.LT.0.0E0) GO TO 290
 
114
        KDFLG = 1
 
115
        Y(I) = CZERO
 
116
        NZ=NZ+1
 
117
        IF (I.EQ.1) GO TO 70
 
118
        IF (Y(I-1).EQ.CZERO) GO TO 70
 
119
        Y(I-1) = CZERO
 
120
        NZ=NZ+1
 
121
   70 CONTINUE
 
122
      I=N
 
123
   75 CONTINUE
 
124
      RZ = CMPLX(2.0E0,0.0E0)/ZR
 
125
      CK = CMPLX(FN,0.0E0)*RZ
 
126
      IB = I+1
 
127
      IF (N.LT.IB) GO TO 160
 
128
C-----------------------------------------------------------------------
 
129
C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
 
130
C     ON UNDERFLOW
 
131
C-----------------------------------------------------------------------
 
132
      FN = FNU+(N-1)
 
133
      IPARD = 1
 
134
      IF (MR.NE.0) IPARD = 0
 
135
      INITD = 0
 
136
      CALL CUNIK(ZR,FN,2,IPARD,TOL,INITD,PHID,ZETA1D,ZETA2D,SUMD,
 
137
     *CWRK(1,3))
 
138
      IF (KODE.EQ.1) GO TO 80
 
139
      CFN=CMPLX(FN,0.0E0)
 
140
      S1=ZETA1D-CFN*(CFN/(ZR+ZETA2D))
 
141
      GO TO 90
 
142
   80 CONTINUE
 
143
      S1=ZETA1D-ZETA2D
 
144
   90 CONTINUE
 
145
      RS1=REAL(S1)
 
146
      IF (ABS(RS1).GT.ELIM) GO TO 95
 
147
      IF (ABS(RS1).LT.ALIM) GO TO 100
 
148
C-----------------------------------------------------------------------
 
149
C     REFINE ESTIMATE AND TEST
 
150
C-----------------------------------------------------------------------
 
151
      APHI=ABS(PHID)
 
152
      RS1=RS1+ALOG(APHI)
 
153
      IF (ABS(RS1).LT.ELIM) GO TO 100
 
154
   95 CONTINUE
 
155
      IF (RS1.GT.0.0E0) GO TO 290
 
156
C-----------------------------------------------------------------------
 
157
C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 
158
C-----------------------------------------------------------------------
 
159
      IF (X.LT.0.0E0) GO TO 290
 
160
      NZ=N
 
161
      DO 96 I=1,N
 
162
        Y(I) = CZERO
 
163
   96 CONTINUE
 
164
      RETURN
 
165
  100 CONTINUE
 
166
C-----------------------------------------------------------------------
 
167
C     RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
 
168
C-----------------------------------------------------------------------
 
169
      S1 = CY(1)
 
170
      S2 = CY(2)
 
171
      C1 = CSR(KFLAG)
 
172
      ASCLE = BRY(KFLAG)
 
173
      DO 120 I=IB,N
 
174
        C2 = S2
 
175
        S2 = CK*S2 + S1
 
176
        S1 = C2
 
177
        CK = CK + RZ
 
178
        C2 = S2*C1
 
179
        Y(I) = C2
 
180
        IF (KFLAG.GE.3) GO TO 120
 
181
        C2R = REAL(C2)
 
182
        C2I = AIMAG(C2)
 
183
        C2R = ABS(C2R)
 
184
        C2I = ABS(C2I)
 
185
        C2M = MAX(C2R,C2I)
 
186
        IF (C2M.LE.ASCLE) GO TO 120
 
187
        KFLAG = KFLAG + 1
 
188
        ASCLE = BRY(KFLAG)
 
189
        S1 = S1*C1
 
190
        S2 = C2
 
191
        S1 = S1*CSS(KFLAG)
 
192
        S2 = S2*CSS(KFLAG)
 
193
        C1 = CSR(KFLAG)
 
194
  120 CONTINUE
 
195
  160 CONTINUE
 
196
      IF (MR.EQ.0) RETURN
 
197
C-----------------------------------------------------------------------
 
198
C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
 
199
C-----------------------------------------------------------------------
 
200
      NZ = 0
 
201
      FMR = MR
 
202
      SGN = -SIGN(PI,FMR)
 
203
C-----------------------------------------------------------------------
 
204
C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
 
205
C-----------------------------------------------------------------------
 
206
      CSGN = CMPLX(0.0E0,SGN)
 
207
      INU = FNU
 
208
      FNF = FNU - INU
 
209
      IFN = INU + N - 1
 
210
      ANG = FNF*SGN
 
211
      CPN = COS(ANG)
 
212
      SPN = SIN(ANG)
 
213
      CSPN = CMPLX(CPN,SPN)
 
214
      IF (MOD(IFN,2).EQ.1) CSPN = -CSPN
 
215
      ASC = BRY(1)
 
216
      KK = N
 
217
      IUF = 0
 
218
      KDFLG = 1
 
219
      IB = IB-1
 
220
      IC = IB-1
 
221
      DO 260 K=1,N
 
222
        FN = FNU + (KK-1)
 
223
C-----------------------------------------------------------------------
 
224
C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
 
225
C     FUNCTION ABOVE
 
226
C-----------------------------------------------------------------------
 
227
        M=3
 
228
        IF (N.GT.2) GO TO 175
 
229
  170   CONTINUE
 
230
        INITD = INIT(J)
 
231
        PHID = PHI(J)
 
232
        ZETA1D = ZETA1(J)
 
233
        ZETA2D = ZETA2(J)
 
234
        SUMD = SUM(J)
 
235
        M = J
 
236
        J = 3 - J
 
237
        GO TO 180
 
238
  175   CONTINUE
 
239
        IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
 
240
        IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 170
 
241
        INITD = 0
 
242
  180   CONTINUE
 
243
        CALL CUNIK(ZR, FN, 1, 0, TOL, INITD, PHID, ZETA1D,
 
244
     *   ZETA2D, SUMD, CWRK(1,M))
 
245
        IF (KODE.EQ.1) GO TO 190
 
246
        CFN = CMPLX(FN,0.0E0)
 
247
        S1 = -ZETA1D + CFN*(CFN/(ZR+ZETA2D))
 
248
        GO TO 200
 
249
  190   CONTINUE
 
250
        S1 = -ZETA1D + ZETA2D
 
251
  200   CONTINUE
 
252
C-----------------------------------------------------------------------
 
253
C     TEST FOR UNDERFLOW AND OVERFLOW
 
254
C-----------------------------------------------------------------------
 
255
        RS1 = REAL(S1)
 
256
        IF (ABS(RS1).GT.ELIM) GO TO 250
 
257
        IF (KDFLG.EQ.1) IFLAG = 2
 
258
        IF (ABS(RS1).LT.ALIM) GO TO 210
 
259
C-----------------------------------------------------------------------
 
260
C     REFINE  TEST AND SCALE
 
261
C-----------------------------------------------------------------------
 
262
        APHI = ABS(PHID)
 
263
        RS1 = RS1 + ALOG(APHI)
 
264
        IF (ABS(RS1).GT.ELIM) GO TO 250
 
265
        IF (KDFLG.EQ.1) IFLAG = 1
 
266
        IF (RS1.LT.0.0E0) GO TO 210
 
267
        IF (KDFLG.EQ.1) IFLAG = 3
 
268
  210   CONTINUE
 
269
        S2 = CSGN*PHID*SUMD
 
270
        C2R = REAL(S1)
 
271
        C2I = AIMAG(S1)
 
272
        C2M = EXP(C2R)*REAL(CSS(IFLAG))
 
273
        S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
 
274
        S2 = S2*S1
 
275
        IF (IFLAG.NE.1) GO TO 220
 
276
        CALL CUCHK(S2, NW, BRY(1), TOL)
 
277
        IF (NW.NE.0) S2 = CMPLX(0.0E0,0.0E0)
 
278
  220   CONTINUE
 
279
        CY(KDFLG) = S2
 
280
        C2 = S2
 
281
        S2 = S2*CSR(IFLAG)
 
282
C-----------------------------------------------------------------------
 
283
C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
 
284
C-----------------------------------------------------------------------
 
285
        S1 = Y(KK)
 
286
        IF (KODE.EQ.1) GO TO 240
 
287
        CALL CS1S2(ZR, S1, S2, NW, ASC, ALIM, IUF)
 
288
        NZ = NZ + NW
 
289
  240   CONTINUE
 
290
        Y(KK) = S1*CSPN + S2
 
291
        KK = KK - 1
 
292
        CSPN = -CSPN
 
293
        IF (C2.NE.CZERO) GO TO 245
 
294
        KDFLG = 1
 
295
        GO TO 260
 
296
  245   CONTINUE
 
297
        IF (KDFLG.EQ.2) GO TO 265
 
298
        KDFLG = 2
 
299
        GO TO 260
 
300
  250   CONTINUE
 
301
        IF (RS1.GT.0.0E0) GO TO 290
 
302
        S2 = CZERO
 
303
        GO TO 220
 
304
  260 CONTINUE
 
305
      K = N
 
306
  265 CONTINUE
 
307
      IL = N - K
 
308
      IF (IL.EQ.0) RETURN
 
309
C-----------------------------------------------------------------------
 
310
C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
 
311
C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
 
312
C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
 
313
C-----------------------------------------------------------------------
 
314
      S1 = CY(1)
 
315
      S2 = CY(2)
 
316
      CS = CSR(IFLAG)
 
317
      ASCLE = BRY(IFLAG)
 
318
      FN = (INU+IL)
 
319
      DO 280 I=1,IL
 
320
        C2 = S2
 
321
        S2 = S1 + CMPLX(FN+FNF,0.0E0)*RZ*S2
 
322
        S1 = C2
 
323
        FN = FN - 1.0E0
 
324
        C2 = S2*CS
 
325
        CK = C2
 
326
        C1 = Y(KK)
 
327
        IF (KODE.EQ.1) GO TO 270
 
328
        CALL CS1S2(ZR, C1, C2, NW, ASC, ALIM, IUF)
 
329
        NZ = NZ + NW
 
330
  270   CONTINUE
 
331
        Y(KK) = C1*CSPN + C2
 
332
        KK = KK - 1
 
333
        CSPN = -CSPN
 
334
        IF (IFLAG.GE.3) GO TO 280
 
335
        C2R = REAL(CK)
 
336
        C2I = AIMAG(CK)
 
337
        C2R = ABS(C2R)
 
338
        C2I = ABS(C2I)
 
339
        C2M = MAX(C2R,C2I)
 
340
        IF (C2M.LE.ASCLE) GO TO 280
 
341
        IFLAG = IFLAG + 1
 
342
        ASCLE = BRY(IFLAG)
 
343
        S1 = S1*CS
 
344
        S2 = CK
 
345
        S1 = S1*CSS(IFLAG)
 
346
        S2 = S2*CSS(IFLAG)
 
347
        CS = CSR(IFLAG)
 
348
  280 CONTINUE
 
349
      RETURN
 
350
  290 CONTINUE
 
351
      NZ = -1
 
352
      RETURN
 
353
      END