      SUBROUTINE ZHECO (H,LDH,N,IPVT,RCOND,Z)
      INTEGER LDH,N,IPVT(N)
      DOUBLE PRECISION RCOND
      DOUBLE COMPLEX H(LDH,N),Z(N)
C
C     *** PURPOSE:
C
C     ZHECO FACTORS A COMPLEX UPPER HESSENBERG MATRIX BY GAUSSIAN
C     ELIMINATION AND ESTIMATES THE CONDITION OF THE MATRIX.  THE
C     LINPACK ANALOGUE OF ZHECO IS ZGECO.
C
C     IF RCOND IS NOT NEEDED, ZHEFA MAY BE SUBSTANTIALLY FASTER.  TO
C     SOLVE THE LINEAR SYSTEM  H*X = B  (I.E., TO COMPUTE
C     (INVERSE(H)*B)), FOLLOW ZHECO BY ZHESL.
C
C     ON ENTRY:
C
C        H     DOUBLE COMPLEX(LDH,N)
C              THE MATRIX TO BE FACTORED.
C
C        LDH   INTEGER
C              THE LEADING OR ROW DIMENSION OF THE ARRAY H AS DECLARED
C              IN THE MAIN CALLING PROGRAM.
C
C        N     INTEGER
C              THE ORDER OF THE MATRIX H.
C
C     ON RETURN:
C
C        H     AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE
C              USED TO OBTAIN IT.  THE FACTORIZATION CAN BE WRITTEN
C              H = L*U  WHERE L IS A PRODUCT OF PERMUTATION AND UNIT
C              LOWER BIDIAGONAL MATRICES AND U IS UPPER TRIANGULAR.
C
C        IPVT  INTEGER(N)
C              AN INTEGER VECTOR OF PIVOT INDICES.
C
C        RCOND DOUBLE PRECISION
C              AN ESTIMATE OF THE RECIPROCAL OF THE CONDITION OF H.
C              FOR THE SYSTEM  H*X = B , RELATIVE PERTURBATIONS
C              IN H AND B OF SIZE EPSILON MAY CAUSE RELATIVE PERTURBA-
C              TIONS IN X OF SIZE  EPSILON/RCOND.  IF RCOND IS SO SMALL
C              THAT THE LOGICAL EXPRESSION    (1.0+RCOND) .EQ. 1.0
C              IS TRUE, THEN H MAY BE SINGULAR TO WORKING PRECISION.
C              IN PARTICULAR, RCOND IS ZERO IF EXACT SINGULARITY IS
C              DETECTED OR THE ESTIMATE UNDERFLOWS.
C
C        Z     DOUBLE COMPLEX(N)
C              A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
C              IF H IS CLOSE TO A SINGULAR MATRIX, THEN Z IS AN
C              APPROXIMATE NULL VECTOR IN THE SENSE THAT
C              NORM(H*Z) = RCOND*NORM(H)*NORM(Z)
C
C     THIS VERSION DATED MAY 1980.
C     ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA.
C
C     SUBROUTINES AND FUNCTIONS CALLED:
C
C     DCL1,ZHEFA,ZL1NRM
C
      DOUBLE PRECISION DCL1
C
C     INTERNAL VARIABLES:
C
      INTEGER I,INFO,J,K,KB,KM1,KP1,L
      DOUBLE PRECISION HNORM,S,SM,YNORM,ZNORM
      DOUBLE COMPLEX EK,T,WK,WKM
C
C     FORTRAN FUNCTIONS CALLED:
C
      DOUBLE COMPLEX DCMPLX,DCONJG
C
C     INTERNAL FUNCTION:
C
      DOUBLE COMPLEX Z1,Z2,CSIGN
      CSIGN(Z1,Z2) = DCL1(Z1)*(Z2/DCL1(Z2))
C
      IF (N.GT.1) GO TO 10
C
C     SPECIAL CASE OF SCALAR (N = 1) SYSTEM
C
      HNORM = DCL1(H(1,1))
      YNORM = HNORM
      GO TO 250
   10 CONTINUE
C
C     COMPUTE COMPLEX 1-NORM OF H
C
      HNORM = 0.0D0
      DO 30 J = 1,N
         S = 0.0D0
         DO 20 I = 1,N
              S = S+DCL1(H(I,J))
   20    CONTINUE
         IF (S .GT. HNORM) HNORM = S
   30 CONTINUE
C
C     FACTOR
C
      CALL ZHEFA (H,LDH,N,IPVT,INFO)
C
C     RCOND = 1/(NORM(H)*(ESTIMATE OF NORM(INVERSE(H))).
C     ESTIMATE = NORM(Z)/NORM(Y)  WHERE  H*Z =Y
C     AND  TRANSPOSE(H)*Y = E.  THE COMPONENTS OF E ARE
C     CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS
C     OF W WHERE  TRANSPOSE(U)*W = E.  THE VECTORS ARE
C     FREQUENTLY RESCALED TO AVOID OVERFLOW.
C
C     SOLVE  TRANSPOSE(U)*W = E
C
      EK = (1.0D0,0.0D0)
      DO 40 J = 1,N
         Z(J) = (0.0D0,0.0D0)
   40 CONTINUE
      DO 130 K = 1,N
         IF (DCL1(Z(K)) .NE. 0.0D0) EK = CSIGN(EK,-Z(K))
         IF (DCL1(EK-Z(K)) .LE. DCL1(H(K,K))) GO TO 60
         S = DCL1(H(K,K))/DCL1(EK-Z(K))
         DO 50 I = 1,N
              Z(I) = DCMPLX(S,0.0D0)*Z(I)
   50    CONTINUE
         EK = DCMPLX(S,0.0D0)*EK
   60    CONTINUE
         WK = EK-Z(K)
         WKM = -EK-Z(K)
         S = DCL1(WK)
         SM = DCL1(WKM)
         IF (DCL1(H(K,K)) .EQ. 0.0D0) GO TO 70
         WK = WK/DCONJG(H(K,K))
         WKM = WKM/DCONJG(H(K,K))
         GO TO 80
   70    CONTINUE
         WK = (1.0D0,0.0D0)
         WKM = (1.0D0,0.0D0)
   80    CONTINUE
         KP1 = K+1
         IF (KP1 .GT. N) GO TO 120
         DO 90 J = KP1,N
              SM = SM+DCL1(Z(J)+WKM*DCONJG(H(K,J)))
              Z(J) = Z(J)+WK*DCONJG(H(K,J))
              S = S+DCL1(Z(J))
   90    CONTINUE
         IF (S .GE. SM) GO TO 110
         T = WKM-WK
         WK = WKM
         DO 100 J = KP1,N
              Z(J) = Z(J)+T*DCONJG(H(K,J))
  100    CONTINUE
  110    CONTINUE
  120    CONTINUE
         Z(K) = WK
  130 CONTINUE
      CALL ZL1NRM (N,Z,ZNORM)
C
C     SOLVE  TRANSPOSE(L)*Y = W
C
      DO 160 KB = 1,N
         K = N+1-KB
         IF (K .LT. N) Z(K) = Z(K)-DCONJG(H(K+1,K))*Z(K+1)
         IF (DCL1(Z(K)) .LE. 1.0D0) GO TO 150
         S = 1.0D0/DCL1(Z(K))
         DO 140 I = 1,N
              Z(I) = DCMPLX(S,0.0D0)*Z(I)
  140    CONTINUE
  150    CONTINUE
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
  160 CONTINUE
      CALL ZL1NRM (N,Z,ZNORM)
C
      YNORM = 1.0D0
C
C     SOLVE  L*V = Y
C
      DO 190 K = 1,N
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
         IF (K .LT. N) Z(K+1) = Z(K+1)-T*H(K+1,K)
         IF (DCL1(Z(K)) .LE. 1.0D0) GO TO 180
         S = 1.0D0/DCL1(Z(K))
         DO 170 I = 1,N
              Z(I) = DCMPLX(S,0.0D0)*Z(I)
  170    CONTINUE
         YNORM = S*YNORM
  180    CONTINUE
  190 CONTINUE
      CALL ZL1NRM (N,Z,ZNORM)
      YNORM = YNORM/ZNORM
C
C     SOLVE  U*Z = V
C
      DO 240 KB = 1,N
         K = N+1-KB
         IF (DCL1(Z(K)) .LE. DCL1(H(K,K))) GO TO 210
         S = DCL1(H(K,K))/DCL1(Z(K))
         DO 200 I = 1,N
              Z(I) = DCMPLX(S,0.0D0)*Z(I)
  200    CONTINUE
         YNORM = S*YNORM
  210    CONTINUE
         IF (DCL1(H(K,K)) .NE. 0.0D0) Z(K) = Z(K)/H(K,K)
         IF (DCL1(H(K,K)) .EQ. 0.0D0) Z(K) = (1.0D0,0.0D0)
         T = -Z(K)
         IF (K .LE. 1) GO TO 230
         KM1 = K-1
         DO 220 I = 1,KM1
              Z(I) = Z(I)+T*H(I,K)
  220    CONTINUE
  230    CONTINUE
  240 CONTINUE
      CALL ZL1NRM (N,Z,ZNORM)
      YNORM = YNORM/ZNORM
  250 CONTINUE
      IF (HNORM .NE. 0.0D0) RCOND = YNORM/HNORM
      IF (HNORM .EQ. 0.0D0) RCOND = 0.0D0
      RETURN
      END
