122 lines
3.0 KiB
Fortran
122 lines
3.0 KiB
Fortran
SUBROUTINE ZKSCL(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
|
|
C***BEGIN PROLOGUE ZKSCL
|
|
C***REFER TO ZBESK
|
|
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***ROUTINES CALLED ZUCHK,AZABS,AZLOG
|
|
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, AZABS,
|
|
* ZDR, ZDI, CELMR, ELM, HELIM, ALAS
|
|
INTEGER I, IC, IDUM, KK, N, NN, NW, NZ
|
|
DIMENSION YR(N), YI(N), CYR(2), CYI(2)
|
|
DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 /
|
|
C
|
|
NZ = 0
|
|
IC = 0
|
|
NN = MIN0(2,N)
|
|
DO 10 I=1,NN
|
|
S1R = YR(I)
|
|
S1I = YI(I)
|
|
CYR(I) = S1R
|
|
CYI(I) = S1I
|
|
AS = AZABS(S1R,S1I)
|
|
ACS = -ZRR + DLOG(AS)
|
|
NZ = NZ + 1
|
|
YR(I) = ZEROR
|
|
YI(I) = ZEROI
|
|
IF (ACS.LT.(-ELIM)) GO TO 10
|
|
CALL AZLOG(S1R, S1I, CSR, CSI, IDUM)
|
|
CSR = CSR - ZRR
|
|
CSI = CSI - ZRI
|
|
STR = DEXP(CSR)/TOL
|
|
CSR = STR*DCOS(CSI)
|
|
CSI = STR*DSIN(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 = DEXP(-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 = AZABS(S2R,S2I)
|
|
ALAS = DLOG(AS)
|
|
ACS = -ZDR + ALAS
|
|
NZ = NZ + 1
|
|
YR(I) = ZEROR
|
|
YI(I) = ZEROI
|
|
IF (ACS.LT.(-ELIM)) GO TO 25
|
|
CALL AZLOG(S2R, S2I, CSR, CSI, IDUM)
|
|
CSR = CSR - ZDR
|
|
CSI = CSI - ZDI
|
|
STR = DEXP(CSR)/TOL
|
|
CSR = STR*DCOS(CSI)
|
|
CSI = STR*DSIN(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
|