GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zkscl.f
Go to the documentation of this file.
1  SUBROUTINE zkscl(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
2 C***BEGIN PROLOGUE ZKSCL
3 C***REFER TO ZBESK
4 C
5 C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
6 C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
7 C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
8 C
9 C***ROUTINES CALLED ZUCHK,XZABS,XZLOG
10 C***END PROLOGUE ZKSCL
11 C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
12  DOUBLE PRECISION acs, as, ascle, cki, ckr, csi, csr, cyi,
13  * cyr, elim, fn, fnu, rzi, rzr, str, s1i, s1r, s2i,
14  * s2r, tol, yi, yr, zeroi, zeror, zri, zrr, xzabs,
15  * zdr, zdi, celmr, elm, helim, alas
16  INTEGER i, ic, idum, kk, n, nn, nw, nz
17  dimension yr(n), yi(n), cyr(2), cyi(2)
18  DATA zeror,zeroi / 0.0d0 , 0.0d0 /
19 C
20  nz = 0
21  ic = 0
22  nn = min0(2,n)
23  DO 10 i=1,nn
24  s1r = yr(i)
25  s1i = yi(i)
26  cyr(i) = s1r
27  cyi(i) = s1i
28  as = xzabs(s1r,s1i)
29  acs = -zrr + dlog(as)
30  nz = nz + 1
31  yr(i) = zeror
32  yi(i) = zeroi
33  IF (acs.LT.(-elim)) go to 10
34  CALL xzlog(s1r, s1i, csr, csi, idum)
35  csr = csr - zrr
36  csi = csi - zri
37  str = dexp(csr)/tol
38  csr = str*dcos(csi)
39  csi = str*dsin(csi)
40  CALL zuchk(csr, csi, nw, ascle, tol)
41  IF (nw.NE.0) go to 10
42  yr(i) = csr
43  yi(i) = csi
44  ic = i
45  nz = nz - 1
46  10 CONTINUE
47  IF (n.EQ.1) RETURN
48  IF (ic.GT.1) go to 20
49  yr(1) = zeror
50  yi(1) = zeroi
51  nz = 2
52  20 CONTINUE
53  IF (n.EQ.2) RETURN
54  IF (nz.EQ.0) RETURN
55  fn = fnu + 1.0d0
56  ckr = fn*rzr
57  cki = fn*rzi
58  s1r = cyr(1)
59  s1i = cyi(1)
60  s2r = cyr(2)
61  s2i = cyi(2)
62  helim = 0.5d0*elim
63  elm = dexp(-elim)
64  celmr = elm
65  zdr = zrr
66  zdi = zri
67 C
68 C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
69 C S2 GETS LARGER THAN EXP(ELIM/2)
70 C
71  DO 30 i=3,n
72  kk = i
73  csr = s2r
74  csi = s2i
75  s2r = ckr*csr - cki*csi + s1r
76  s2i = cki*csr + ckr*csi + s1i
77  s1r = csr
78  s1i = csi
79  ckr = ckr + rzr
80  cki = cki + rzi
81  as = xzabs(s2r,s2i)
82  alas = dlog(as)
83  acs = -zdr + alas
84  nz = nz + 1
85  yr(i) = zeror
86  yi(i) = zeroi
87  IF (acs.LT.(-elim)) go to 25
88  CALL xzlog(s2r, s2i, csr, csi, idum)
89  csr = csr - zdr
90  csi = csi - zdi
91  str = dexp(csr)/tol
92  csr = str*dcos(csi)
93  csi = str*dsin(csi)
94  CALL zuchk(csr, csi, nw, ascle, tol)
95  IF (nw.NE.0) go to 25
96  yr(i) = csr
97  yi(i) = csi
98  nz = nz - 1
99  IF (ic.EQ.kk-1) go to 40
100  ic = kk
101  go to 30
102  25 CONTINUE
103  IF(alas.LT.helim) go to 30
104  zdr = zdr - elim
105  s1r = s1r*celmr
106  s1i = s1i*celmr
107  s2r = s2r*celmr
108  s2i = s2i*celmr
109  30 CONTINUE
110  nz = n
111  IF(ic.EQ.n) nz=n-1
112  go to 45
113  40 CONTINUE
114  nz = kk - 2
115  45 CONTINUE
116  DO 50 i=1,nz
117  yr(i) = zeror
118  yi(i) = zeroi
119  50 CONTINUE
120  RETURN
121  END