GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ckscl.f
Go to the documentation of this file.
1  SUBROUTINE ckscl(ZR, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
2 C***BEGIN PROLOGUE CKSCL
3 C***REFER TO CBKNU,CUNK1,CUNK2
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 CUCHK
10 C***END PROLOGUE CKSCL
11  COMPLEX ck, cs, cy, czero, rz, s1, s2, y, zr, zd, celm
12  REAL aa, ascle, acs, as, csi, csr, elim, fn, fnu, tol, xx, zri,
13  * elm, alas, helim
14  INTEGER i, ic, k, kk, n, nn, nw, nz
15  dimension y(n), cy(2)
16  DATA czero / (0.0e0,0.0e0) /
17 C
18  nz = 0
19  ic = 0
20  xx = REAL(zr)
21  nn = min0(2,n)
22  DO 10 i=1,nn
23  s1 = y(i)
24  cy(i) = s1
25  as = cabs(s1)
26  acs = -xx + alog(as)
27  nz = nz + 1
28  y(i) = czero
29  IF (acs.LT.(-elim)) go to 10
30  cs = -zr + clog(s1)
31  csr = REAL(cs)
32  csi = aimag(cs)
33  aa = exp(csr)/tol
34  cs = cmplx(aa,0.0e0)*cmplx(cos(csi),sin(csi))
35  CALL cuchk(cs, nw, ascle, tol)
36  IF (nw.NE.0) go to 10
37  y(i) = cs
38  nz = nz - 1
39  ic = i
40  10 CONTINUE
41  IF (n.EQ.1) RETURN
42  IF (ic.GT.1) go to 20
43  y(1) = czero
44  nz = 2
45  20 CONTINUE
46  IF (n.EQ.2) RETURN
47  IF (nz.EQ.0) RETURN
48  fn = fnu + 1.0e0
49  ck = cmplx(fn,0.0e0)*rz
50  s1 = cy(1)
51  s2 = cy(2)
52  helim = 0.5e0*elim
53  elm = exp(-elim)
54  celm = cmplx(elm,0.0e0)
55  zri =aimag(zr)
56  zd = zr
57 C
58 C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
59 C S2 GETS LARGER THAN EXP(ELIM/2)
60 C
61  DO 30 i=3,n
62  kk = i
63  cs = s2
64  s2 = ck*s2 + s1
65  s1 = cs
66  ck = ck + rz
67  as = cabs(s2)
68  alas = alog(as)
69  acs = -xx + alas
70  nz = nz + 1
71  y(i) = czero
72  IF (acs.LT.(-elim)) go to 25
73  cs = -zd + clog(s2)
74  csr = REAL(cs)
75  csi = aimag(cs)
76  aa = exp(csr)/tol
77  cs = cmplx(aa,0.0e0)*cmplx(cos(csi),sin(csi))
78  CALL cuchk(cs, nw, ascle, tol)
79  IF (nw.NE.0) go to 25
80  y(i) = cs
81  nz = nz - 1
82  IF (ic.EQ.(kk-1)) go to 40
83  ic = kk
84  go to 30
85  25 CONTINUE
86  IF(alas.LT.helim) go to 30
87  xx = xx-elim
88  s1 = s1*celm
89  s2 = s2*celm
90  zd = cmplx(xx,zri)
91  30 CONTINUE
92  nz = n
93  IF(ic.EQ.n) nz=n-1
94  go to 45
95  40 CONTINUE
96  nz = kk - 2
97  45 CONTINUE
98  DO 50 k=1,nz
99  y(k) = czero
100  50 CONTINUE
101  RETURN
102  END