      SUBROUTINE CPSI (Z, W)
C     ******************************************************************
C     EVALUATION OF THE COMPLEX DIGAMMA FUNCTION
C     ******************************************************************
      COMPLEX Z, W
      REAL COEFF(18)
      DOUBLE PRECISION DS1, DS2
C     -------------------
C     DPI = 2*PI
C     -------------------
      DATA PI/3.14159265358979324/
      DATA DPI/6.28318530717958648/
C     -------------------
C     COEFFICIENTS OF A RATIONAL APPROXIMATION FOR 1.0 - EXP(U)
C     WHEN -0.1054 .LE. U .LE. 0.0
C     -------------------
      DATA F0/840.07385296052619/,F1/20.001230821894200/
      DATA D0/1680.1477059210524/,D1/180.01477047052042/
C     -------------------
C     COEFFICIENTS FOR THE ASYMPTOTIC EXPANSION
C     -------------------
      DATA COEFF(1) /+8.333333333333333E-2/,
     1     COEFF(2) /-8.333333333333333E-3/,
     2     COEFF(3) /+3.968253968253968E-3/,
     3     COEFF(4) /-4.166666666666667E-3/,
     4     COEFF(5) /+7.575757575757576E-3/,
     5     COEFF(6) /-2.109279609279609E-2/,
     6     COEFF(7) /+8.333333333333333E-2/,
     7     COEFF(8) /-4.432598039215686E-1/,
     8     COEFF(9) /+3.053954330270120E00/,
     9     COEFF(10)/-2.645621212121212E+1/,
     1     COEFF(11)/+2.814601449275362E+2/,
     2     COEFF(12)/-3.607510546398046E+3/,
     3     COEFF(13)/+5.482758333333333E+4/,
     4     COEFF(14)/-9.749368238505747E+5/,
     5     COEFF(15)/+2.005269579668808E+7/,
     6     COEFF(16)/-4.723848677216299E+8/,
     7     COEFF(17)/+1.263572479591667E10/,
     8     COEFF(18)/-3.808793112524537E11/
C     -------------------
      X = REAL(Z)
      Y = AIMAG(Z)
      IF (X .GE. 0.0) GO TO 40
C     ------------------------------------------------------------------
C     CASE WHEN THE REAL PART OF Z IS NEGATIVE
C     ------------------------------------------------------------------
      Y = ABS(Y)
      T = -DPI*Y
      ET = EXP(T)
C
C     SET  A1 = (1 + ET)/2  AND  A2 = (1 - ET)/2
C
      A1 = 0.5*(1.0 + ET)
      IF (T .GE. -0.1054) GO TO 10
      A2 = 0.5*(1.0 - ET)
      GO TO 20
   10 T2 = T*T
      C = -T*(F0 + F1*T2)
      A2 = C/((T2 + D1)*T2 + D0 + C)
C
C     COMPUTE SIN(PI*X) AND COS(PI*X), OR -SIN(PI*X) AND -COS(PI*X)
C
   20 K = -X
      U = X + K
      IF (U .LE. -0.5) U = U + 1.0
      U = PI*U
      SN = SIN(U)
      CN = COS(U)
C
C     SET H1 + H2*I = PI*COT(PI*Z)
C
      S1 = A1*SN
      S2 = A2*CN
      C1 = A1*CN
      C2 = -A2*SN
      S = S1*S1 + S2*S2
      H1 = PI*(S1*C1 + S2*C2)/S
      H2 = PI*(S1*C2 - S2*C1)/S
C
      IF (AIMAG(Z) .LT. 0.0) GO TO 30
      X = 1.0 - X
      Y = -Y
      GO TO 40
   30 H2 = -H2
      X = 1.0 - X
C     ------------------------------------------------------------------
C     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
C     ------------------------------------------------------------------
   40 T = X
      Y2 = Y*Y
C
C     LET S1 + S2*I BE THE SUM OF THE TERMS 1/(Z+J) FOR J = 0,1,...,N-1
C
      DS1 = 0.D0
      DS2 = 0.D0
   50 A = T*T + Y2
      IF (A .GE. 64.0) GO TO 51
      DS1 = DS1 + DBLE(T/A)
      DS2 = DS2 - DBLE(Y/A)
      T = T + 1.0
      GO TO 50
   51 S1 = DS1
      S2 = DS2
C
C     SET W1 + W2*I = LOG(Z+N)
C
      W1 = 0.5*ALOG(A)
      W2 = ATAN2(Y,T)
C
C     LET A1 + A2*I BE THE ASYMPTOTIC SUM
C
      U1 = T/A
      U2 = -Y/A
      Q1 = U1*U1 - U2*U2
      Q2 = 2.0*U1*U2
      V1 = Q1
      V2 = Q2
      A1 = 0.0
      A2 = 0.0
      SM1 = 0.0
      SM2 = 0.0
      DO 61 J = 1,18
      T1 = COEFF(J) * V1
      T2 = COEFF(J) * V2
      SM1 = SM1 + T1
      SM2 = SM2 + T2
      IF (A1 .NE. SM1) GO TO 60
      IF (A2 .EQ. SM2) GO TO 70
   60 A1 = SM1
      A2 = SM2
      T1 = V1*Q1 - V2*Q2
      T2 = V1*Q2 + V2*Q1
      V1 = T1
   61 V2 = T2
C     ------------------------------------------------------------------
C     FINAL ASSEMBLY
C     ------------------------------------------------------------------
   70 A1 = 0.5*U1 + A1
      A2 = 0.5*U2 + A2
      W1 = (W1 - S1) - A1
      W2 = (W2 - A2) - S2
      W  = CMPLX(W1,W2)
      IF (REAL(Z) .GE. 0.0) RETURN
      W = CMPLX(W1 - H1, W2 - H2)
      RETURN
      END
