1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
*DECK ZKSCL
SUBROUTINE ZKSCL (ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE,
+ TOL, ELIM)
C***BEGIN PROLOGUE ZKSCL
C***SUBSIDIARY
C***PURPOSE Subsidiary to ZBESK
C***LIBRARY SLATEC
C***TYPE ALL (CKSCL-A, ZKSCL-A)
C***AUTHOR Amos, D. E., (SNL)
C***DESCRIPTION
C
C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
C
C***SEE ALSO ZBESK
C***ROUTINES CALLED ZABS, ZLOG, ZUCHK
C***REVISION HISTORY (YYMMDD)
C 830501 DATE WRITTEN
C 910415 Prologue converted to Version 4.0 format. (BAB)
C 930122 Added ZLOG to EXTERNAL statement. (RWC)
C***END PROLOGUE ZKSCL
C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
DOUBLE PRECISION ACS, AS, ASCLE, CKI, CKR, CSI, CSR, CYI,
* CYR, ELIM, FN, FNU, RZI, RZR, STR, S1I, S1R, S2I,
* S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, ZABS,
* ZDR, ZDI, CELMR, ELM, HELIM, ALAS
INTEGER I, IC, IDUM, KK, N, NN, NW, NZ
DIMENSION YR(N), YI(N), CYR(2), CYI(2)
EXTERNAL ZABS, ZLOG
DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 /
C***FIRST EXECUTABLE STATEMENT ZKSCL
NZ = 0
IC = 0
NN = MIN(2,N)
DO 10 I=1,NN
S1R = YR(I)
S1I = YI(I)
CYR(I) = S1R
CYI(I) = S1I
AS = ZABS(S1R,S1I)
ACS = -ZRR + LOG(AS)
NZ = NZ + 1
YR(I) = ZEROR
YI(I) = ZEROI
IF (ACS.LT.(-ELIM)) GO TO 10
CALL ZLOG(S1R, S1I, CSR, CSI, IDUM)
CSR = CSR - ZRR
CSI = CSI - ZRI
STR = EXP(CSR)/TOL
CSR = STR*COS(CSI)
CSI = STR*SIN(CSI)
CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 10
YR(I) = CSR
YI(I) = CSI
IC = I
NZ = NZ - 1
10 CONTINUE
IF (N.EQ.1) RETURN
IF (IC.GT.1) GO TO 20
YR(1) = ZEROR
YI(1) = ZEROI
NZ = 2
20 CONTINUE
IF (N.EQ.2) RETURN
IF (NZ.EQ.0) RETURN
FN = FNU + 1.0D0
CKR = FN*RZR
CKI = FN*RZI
S1R = CYR(1)
S1I = CYI(1)
S2R = CYR(2)
S2I = CYI(2)
HELIM = 0.5D0*ELIM
ELM = EXP(-ELIM)
CELMR = ELM
ZDR = ZRR
ZDI = ZRI
C
C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
C S2 GETS LARGER THAN EXP(ELIM/2)
C
DO 30 I=3,N
KK = I
CSR = S2R
CSI = S2I
S2R = CKR*CSR - CKI*CSI + S1R
S2I = CKI*CSR + CKR*CSI + S1I
S1R = CSR
S1I = CSI
CKR = CKR + RZR
CKI = CKI + RZI
AS = ZABS(S2R,S2I)
ALAS = LOG(AS)
ACS = -ZDR + ALAS
NZ = NZ + 1
YR(I) = ZEROR
YI(I) = ZEROI
IF (ACS.LT.(-ELIM)) GO TO 25
CALL ZLOG(S2R, S2I, CSR, CSI, IDUM)
CSR = CSR - ZDR
CSI = CSI - ZDI
STR = EXP(CSR)/TOL
CSR = STR*COS(CSI)
CSI = STR*SIN(CSI)
CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 25
YR(I) = CSR
YI(I) = CSI
NZ = NZ - 1
IF (IC.EQ.KK-1) GO TO 40
IC = KK
GO TO 30
25 CONTINUE
IF(ALAS.LT.HELIM) GO TO 30
ZDR = ZDR - ELIM
S1R = S1R*CELMR
S1I = S1I*CELMR
S2R = S2R*CELMR
S2I = S2I*CELMR
30 CONTINUE
NZ = N
IF(IC.EQ.N) NZ=N-1
GO TO 45
40 CONTINUE
NZ = KK - 2
45 CONTINUE
DO 50 I=1,NZ
YR(I) = ZEROR
YI(I) = ZEROI
50 CONTINUE
RETURN
END
|
