github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zkscl.f (about)

     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,ZABS,ZLOG
    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, ZABS,
    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 = ZABS(CMPLX(S1R,S1I,kind=KIND(1.0D0)))
    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 ZLOG(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 = ZABS(CMPLX(S2R,S2I,kind=KIND(1.0D0)))
    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 ZLOG(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