      SUBROUTINE WGECO (AR, AI, LDA, N, IPVT, RCOND, ZR, ZI)
      IMPLICIT NONE
C
C FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION
C   AND ESTIMATES THE CONDITION OF THE MATRIX
C
C IF RCOND IS NOT NEEDED, WGEFA IS SLIGHTLY FASTER
C   TO SOLVE    A*X = B , FOLLOW WGECO BY WGESL
C   TO COMPUTE  INVERSE(A)*C , FOLLOW WGECO BY WGESL
C   TO COMPUTE  DETERMINANT(A) , FOLLOW WGECO BY WGEDI
C   TO COMPUTE  INVERSE(A) , FOLLOW WGECO BY WGEDI
C
C ON ENTRY:
C    A    DOUBLE-COMPLEX(LDA, N),  MATRIX TO BE FACTORED
C    LDA  INTEGER,  LEADING DIMENSION OF THE ARRAY A
C    N    INTEGER,  ORDER OF THE MATRIX A
C
C ON RETURN:
C    A     AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE USED TO
C            OBTAIN IT.  THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C            L IS A PRODUCT OF PERMUTATION AND UNIT LOWER TRIANGULAR MATRICES
C            AND U IS UPPER TRIANGULAR.
C    IPVT  INTEGER(N),  AN INTEGER VECTOR OF PIVOT INDICES.
C    RCOND DOUBLE PRECISION,  AN ESTIMATE OF THE RECIPROCAL CONDITION OF A.
C            FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS IN A AND B OF
C            SIZE EPSILON MAY CAUSE RELATIVE PERTURBATIONS IN X OF SIZE
C            EPSILON/RCOND.  IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION
C            1.0 + RCOND .EQ. 1.0  IS TRUE, THEN A MAY BE SINGULAR TO WORKING
C            PRECISION.  IN PARTICULAR, RCOND IS ZERO IF EXACT SINGULARITY IS
C            DETECTED OR THE ESTIMATE UNDERFLOWS.
C    Z     DOUBLE-COMPLEX(N),  WORK VECTOR WHOSE CONTENTS ARE USUALLY
C            UNIMPORTANT.  IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS AN
C            APPROXIMATE NULL VECTOR IN THE SENSE THAT NORM(A*Z) =
C            RCOND*NORM(A)*NORM(Z).
C
C VERSION 07/01/79. CLEVE MOLER, UNIV OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
      INTEGER LDA, N, IPVT(*)
      DOUBLE PRECISION AR(LDA,*), AI(LDA,*), RCOND, ZR(*), ZI(*)
C
      DOUBLE PRECISION EKR, EKI, TR, TI, WKR, WKI, WKMR, WKMI
      DOUBLE PRECISION ANORM, S, SM, YNORM
      INTEGER INFO, J, K, KB, KP1, L
      DOUBLE PRECISION ZDUMR, ZDUMI
C
      DOUBLE PRECISION CABS1, WDOTCR, WDOTCI, WASUM, FLOP
C
C
      CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
C
C
C ***      COMPUTE 1-NORM OF A
      ANORM = 0.0D0
      DO 10 J = 1, N
        ANORM = DMAX1 (ANORM, WASUM (N, AR(1,J), AI(1,J), 1))
10    CONTINUE
C
C ***      FACTOR
      CALL WGEFA (AR, AI, LDA, N, IPVT, INFO)
C
C *** RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A))))
C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  CTRANS(A)*Y = E
C     CTRANS(A)  IS THE CONJUGATE TRANSPOSE OF A
C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN
C       THE ELEMENTS OF W  WHERE  CTRANS(U)*W = E
C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW
C
C ***      SOLVE CTRANS(U)*W = E
      EKR = 1.0D0
      EKI = 0.0D0
      DO 20 J = 1, N
        ZR(J) = 0.0D0
        ZI(J) = 0.0D0
20    CONTINUE
      DO 110 K = 1, N
        CALL WSIGN (EKR, EKI, -ZR(K), -ZI(K), EKR, EKI)
        IF (CABS1 (EKR-ZR(K), EKI-ZI(K)).LE.CABS1 (AR(K,K), AI(K,K)))
     .      GO TO 40
        S = CABS1 (AR(K,K), AI(K,K))/CABS1 (EKR-ZR(K), EKI-ZI(K))
        CALL WRSCAL (N, S, ZR, ZI, 1)
        EKR = S*EKR
        EKI = S*EKI
40      CONTINUE
        WKR = EKR-ZR(K)
        WKI = EKI-ZI(K)
        WKMR = -EKR-ZR(K)
        WKMI = -EKI-ZI(K)
        S = CABS1 (WKR, WKI)
        SM = CABS1 (WKMR, WKMI)
        IF (CABS1 (AR(K,K), AI(K,K)).EQ.0.0D0) GO TO 50
        CALL WDIV (WKR, WKI, AR(K,K), -AI(K,K), WKR, WKI)
        CALL WDIV (WKMR, WKMI, AR(K,K), -AI(K,K), WKMR, WKMI)
        GO TO 60
C
50      CONTINUE
        WKR = 1.0D0
        WKI = 0.0D0
        WKMR = 1.0D0
        WKMI = 0.0D0
60      CONTINUE
        KP1 = K+1
        IF (KP1.GT.N) GO TO 100
        DO 70 J = KP1, N
          CALL WMUL (WKMR, WKMI, AR(K,J), -AI(K,J), TR, TI)
          SM = FLOP (SM+CABS1 (ZR(J)+TR, ZI(J)+TI))
          CALL WAXPY (1, WKR, WKI, AR(K,J), -AI(K,J),
     .                1, ZR(J), ZI(J), 1)
          S = FLOP (S+CABS1 (ZR(J), ZI(J)))
70      CONTINUE
        IF (S.GE.SM) GO TO 90
        TR = WKMR-WKR
        TI = WKMI-WKI
        WKR = WKMR
        WKI = WKMI
        DO 80 J = KP1, N
          CALL WAXPY (1, TR, TI, AR(K,J), -AI(K,J), 1, ZR(J), ZI(J), 1)
80      CONTINUE
90      CONTINUE
100     CONTINUE
        ZR(K) = WKR
        ZI(K) = WKI
110   CONTINUE
      S = 1.0D0/WASUM (N, ZR, ZI, 1)
      CALL WRSCAL (N, S, ZR, ZI, 1)
C
C ***      SOLVE CTRANS(L)*Y = W
      DO 140 KB = 1, N
        K = N+1-KB
        IF (K.GE.N) GO TO 120
        ZR(K) = ZR(K)+WDOTCR (N-K, AR(K+1,K), AI(K+1,K),
     .                        1, ZR(K+1), ZI(K+1), 1)
        ZI(K) = ZI(K)+WDOTCI (N-K, AR(K+1,K), AI(K+1,K),
     .                        1, ZR(K+1), ZI(K+1), 1)
120     CONTINUE
        IF (CABS1 (ZR(K), ZI(K)).LE.1.0D0) GO TO 130
        S = 1.0D0/CABS1 (ZR(K), ZI(K))
        CALL WRSCAL (N, S, ZR, ZI, 1)
130     CONTINUE
        L = IPVT(K)
        TR = ZR(L)
        TI = ZI(L)
        ZR(L) = ZR(K)
        ZI(L) = ZI(K)
        ZR(K) = TR
        ZI(K) = TI
140   CONTINUE
      S = 1.0D0/WASUM (N, ZR, ZI, 1)
      CALL WRSCAL (N, S, ZR, ZI, 1)
      YNORM = 1.0D0
C
C ***      SOLVE L*V = Y
      DO 160 K = 1, N
        L = IPVT(K)
        TR = ZR(L)
        TI = ZI(L)
        ZR(L) = ZR(K)
        ZI(L) = ZI(K)
        ZR(K) = TR
        ZI(K) = TI
        IF (K.LT.N) CALL WAXPY (N-K, TR, TI, AR(K+1,K), AI(K+1,K),
     .                          1, ZR(K+1), ZI(K+1), 1)
        IF (CABS1 (ZR(K), ZI(K)).LE.1.0D0) GO TO 150
        S = 1.0D0/CABS1 (ZR(K), ZI(K))
        CALL WRSCAL (N, S, ZR, ZI, 1)
        YNORM = S*YNORM
150     CONTINUE
160   CONTINUE
      S = 1.0D0/WASUM (N, ZR, ZI, 1)
      CALL WRSCAL (N, S, ZR, ZI, 1)
      YNORM = S*YNORM
C
C ***      SOLVE  U*Z = V
      DO 200 KB = 1, N
        K = N+1-KB
        IF (CABS1 (ZR(K), ZI(K)).LE.CABS1 (AR(K,K), AI(K,K))) GO TO 170
        S = CABS1 (AR(K,K), AI(K,K))/CABS1 (ZR(K), ZI(K))
        CALL WRSCAL (N, S, ZR, ZI, 1)
        YNORM = S*YNORM
170     CONTINUE
        IF (CABS1 (AR(K,K), AI(K,K)).EQ.0.0D0) GO TO 180
        CALL WDIV (ZR(K), ZI(K), AR(K,K), AI(K,K), ZR(K), ZI(K))
180     CONTINUE
        IF (CABS1 (AR(K,K), AI(K,K)).NE.0.0D0) GO TO 190
        ZR(K) = 1.0D0
        ZI(K) = 0.0D0
190     CONTINUE
        TR = -ZR(K)
        TI = -ZI(K)
        CALL WAXPY (K-1, TR, TI, AR(1,K), AI(1,K), 1, ZR(1), ZI(1), 1)
200   CONTINUE
C
C ***      MAKE ZNORM = 1.0
      S = 1.0D0/WASUM (N, ZR, ZI, 1)
      CALL WRSCAL (N, S, ZR, ZI, 1)
      YNORM = S*YNORM
      IF (ANORM.NE.0.0D0) RCOND = YNORM/ANORM
      IF (ANORM.EQ.0.0D0) RCOND = 0.0D0
C
      RETURN
      END
