GNU Octave: libcruft/amos/zkscl.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZKSCL(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
00002 C***BEGIN PROLOGUE  ZKSCL
00003 C***REFER TO  ZBESK
00004 C
00005 C     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
00006 C     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
00007 C     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
00008 C
00009 C***ROUTINES CALLED  ZUCHK,XZABS,XZLOG
00010 C***END PROLOGUE  ZKSCL
00011 C     COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
00012       DOUBLE PRECISION ACS, AS, ASCLE, CKI, CKR, CSI, CSR, CYI,
00013      * CYR, ELIM, FN, FNU, RZI, RZR, STR, S1I, S1R, S2I,
00014      * S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, XZABS,
00015      * ZDR, ZDI, CELMR, ELM, HELIM, ALAS
00016       INTEGER I, IC, IDUM, KK, N, NN, NW, NZ
00017       DIMENSION YR(N), YI(N), CYR(2), CYI(2)
00018       DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 /
00019 C
00020       NZ = 0
00021       IC = 0
00022       NN = MIN0(2,N)
00023       DO 10 I=1,NN
00024         S1R = YR(I)
00025         S1I = YI(I)
00026         CYR(I) = S1R
00027         CYI(I) = S1I
00028         AS = XZABS(S1R,S1I)
00029         ACS = -ZRR + DLOG(AS)
00030         NZ = NZ + 1
00031         YR(I) = ZEROR
00032         YI(I) = ZEROI
00033         IF (ACS.LT.(-ELIM)) GO TO 10
00034         CALL XZLOG(S1R, S1I, CSR, CSI, IDUM)
00035         CSR = CSR - ZRR
00036         CSI = CSI - ZRI
00037         STR = DEXP(CSR)/TOL
00038         CSR = STR*DCOS(CSI)
00039         CSI = STR*DSIN(CSI)
00040         CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
00041         IF (NW.NE.0) GO TO 10
00042         YR(I) = CSR
00043         YI(I) = CSI
00044         IC = I
00045         NZ = NZ - 1
00046    10 CONTINUE
00047       IF (N.EQ.1) RETURN
00048       IF (IC.GT.1) GO TO 20
00049       YR(1) = ZEROR
00050       YI(1) = ZEROI
00051       NZ = 2
00052    20 CONTINUE
00053       IF (N.EQ.2) RETURN
00054       IF (NZ.EQ.0) RETURN
00055       FN = FNU + 1.0D0
00056       CKR = FN*RZR
00057       CKI = FN*RZI
00058       S1R = CYR(1)
00059       S1I = CYI(1)
00060       S2R = CYR(2)
00061       S2I = CYI(2)
00062       HELIM = 0.5D0*ELIM
00063       ELM = DEXP(-ELIM)
00064       CELMR = ELM
00065       ZDR = ZRR
00066       ZDI = ZRI
00067 C
00068 C     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
00069 C     S2 GETS LARGER THAN EXP(ELIM/2)
00070 C
00071       DO 30 I=3,N
00072         KK = I
00073         CSR = S2R
00074         CSI = S2I
00075         S2R = CKR*CSR - CKI*CSI + S1R
00076         S2I = CKI*CSR + CKR*CSI + S1I
00077         S1R = CSR
00078         S1I = CSI
00079         CKR = CKR + RZR
00080         CKI = CKI + RZI
00081         AS = XZABS(S2R,S2I)
00082         ALAS = DLOG(AS)
00083         ACS = -ZDR + ALAS
00084         NZ = NZ + 1
00085         YR(I) = ZEROR
00086         YI(I) = ZEROI
00087         IF (ACS.LT.(-ELIM)) GO TO 25
00088         CALL XZLOG(S2R, S2I, CSR, CSI, IDUM)
00089         CSR = CSR - ZDR
00090         CSI = CSI - ZDI
00091         STR = DEXP(CSR)/TOL
00092         CSR = STR*DCOS(CSI)
00093         CSI = STR*DSIN(CSI)
00094         CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
00095         IF (NW.NE.0) GO TO 25
00096         YR(I) = CSR
00097         YI(I) = CSI
00098         NZ = NZ - 1
00099         IF (IC.EQ.KK-1) GO TO 40
00100         IC = KK
00101         GO TO 30
00102    25   CONTINUE
00103         IF(ALAS.LT.HELIM) GO TO 30
00104         ZDR = ZDR - ELIM
00105         S1R = S1R*CELMR
00106         S1I = S1I*CELMR
00107         S2R = S2R*CELMR
00108         S2I = S2I*CELMR
00109    30 CONTINUE
00110       NZ = N
00111       IF(IC.EQ.N) NZ=N-1
00112       GO TO 45
00113    40 CONTINUE
00114       NZ = KK - 2
00115    45 CONTINUE
00116       DO 50 I=1,NZ
00117         YR(I) = ZEROR
00118         YI(I) = ZEROI
00119    50 CONTINUE
00120       RETURN
00121       END