      SUBROUTINE WQRDC (XR, XI, LDX, N, P, QRAUXR, QRAUXI,
     .                  JPVT, WORKR, WORKI, JOB)
      IMPLICIT NONE
C
C USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR FACTORIZATION OF
C  AN N BY P MATRIX X.  COLUMN PIVOTING BASED ON THE 2-NORMS OF THE REDUCED
C  COLUMNS MAY BE PERFORMED AT THE USERS OPTION
C
C ON ENTRY:
C    X     DOUBLE-COMPLEX(LDX,P), WHERE LDX .GE. N.   CONTAINS THE MATRIX
C           WHOSE DECOMPOSITION IS TO BE COMPUTED
C    LDX   INTEGER,  LEADING DIMENSION OF THE ARRAY X
C    N     INTEGER,  NUMBER OF ROWS OF THE MATRIX X
C    P     INTEGER,  NUMBER OF COLUMNS OF THE MATRIX X
C    JPVT  INTEGER(P),  CONTAINS INTEGERS THAT CONTROL THE SELECTION OF THE
C           PIVOT COLUMNS.  THE K-TH COLUMN X(K) OF X IS PLACED IN ONE OF
C           THREE CLASSES ACCORDING TO THE VALUE OF JPVT(K):
C             IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL COLUMN
C             IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN
C             IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN
C           BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS ARE MOVED
C           TO THE BEGINNING OF THE ARRAY X AND FINAL COLUMNS TO THE END.
C           BOTH INITIAL AND FINAL COLUMNS ARE FROZEN IN PLACE DURING THE
C           COMPUTATION AND ONLY FREE COLUMNS ARE MOVED.  AT THE K-TH STAGE
C           OF THE REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN IT IS
C           INTERCHANGED WITH THE FREE COLUMN OF LARGEST REDUCED NORM.
C           JPVT IS NOT REFERENCED IF JOB .EQ. 0
C    WORK  DOUBLE-COMPLEX(P),  WORK ARRAY.  NOT REFERENCED IF JOB .EQ. 0
C    JOB   INTEGER,  INITIATES COLUMN PIVOTING:
C           IF JOB .EQ. 0, NO PIVOTING IS DONE
C           IF JOB .NE. 0, PIVOTING IS DONE
C
C ON RETURN:
C    X      CONTAINS IN ITS UPPER TRIANGLE THE UPPER TRIANGULAR MATRIX R OF
C            THE QR FACTORIZATION.  BELOW ITS DIAGONAL X CONTAINS INFORMATION
C            FROM WHICH THE UNITARY PART OF THE DECOMPOSITION CAN BE RECOVERED.
C            NOTE THAT IF PIVOTING HAS BEEN REQUESTED, THE DECOMPOSITION IS
C            NOT THAT OF THE ORIGINAL MATRIX X BUT THAT OF X WITH ITS COLUMNS
C            PERMUTED AS DESCRIBED BY JPVT
C    QRAUX  DOUBLE-COMPLEX(P),  CONTAINS FURTHER INFORMATION REQUIRED TO
C            RECOVER THE UNITARY PART OF THE DECOMPOSITION
C    JPVT   JPVT(K) =  INDEX OF THE COLUMN OF THE ORIGINAL MATRIX THAT HAS
C            BEEN INTERCHANGED INTO THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
C
C VERSION 07/03/79 C.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
      INTEGER LDX, N, P, JOB, JPVT(*)
      DOUBLE PRECISION XR(LDX,*), XI(LDX,*), QRAUXR(*),
     .                 QRAUXI(*), WORKR(*), WORKI(*)
C
      INTEGER J, JP, L, LP1, LUP, MAXJ, PL, PU, JJ
      DOUBLE PRECISION MAXNRM, TT, NRMXLR, NRMXLI, TR, TI
      LOGICAL NEGJ, SWAPJ
      DOUBLE PRECISION ZDUMR, ZDUMI
C
      DOUBLE PRECISION CABS1, PYTHAG, WDOTCR, WDOTCI, FLOP, WNRM2
C
C
      CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
C
C
      PL = 1
      PU = 0
      IF (JOB.EQ.0) GO TO 60
C
C ***      PIVOTING HAS BEEN REQUESTED.  REARRANGE COLUMNS ACCORDING TO JPVT
      DO 20 J = 1, P
        SWAPJ = JPVT(J).GT.0
        NEGJ = JPVT(J).LT.0
        JPVT(J) = J
        IF (NEGJ) JPVT(J) = -J
        IF (.NOT.SWAPJ) GO TO 10
        IF (J.NE.PL) CALL WSWAP (N, XR(1,PL), XI(1,PL),
     .                           1, XR(1,J), XI(1,J), 1)
        JPVT(J) = JPVT(PL)
        JPVT(PL) = J
        PL = PL+1
10      CONTINUE
20    CONTINUE
      PU = P
      DO 50 JJ = 1, P
        J = P-JJ+1
        IF (JPVT(J).GE.0) GO TO 40
        JPVT(J) = -JPVT(J)
        IF (J.EQ.PU) GO TO 30
        CALL WSWAP (N, XR(1,PU), XI(1,PU), 1, XR(1,J), XI(1,J), 1)
        JP = JPVT(PU)
        JPVT(PU) = JPVT(J)
        JPVT(J) = JP
30      CONTINUE
        PU = PU-1
40      CONTINUE
50    CONTINUE
C
C ***      COMPUTE THE NORMS OF THE FREE COLUMNS
60    CONTINUE
      IF (PU.LT.PL) GO TO 80
      DO 70 J = PL, PU
        QRAUXR(J) = WNRM2 (N, XR(1,J), XI(1,J), 1)
        QRAUXI(J) = 0.0D0
        WORKR(J) = QRAUXR(J)
        WORKI(J) = QRAUXI(J)
70    CONTINUE
C
C ***      PERFORM THE HOUSEHOLDER REDUCTION OF X
80    CONTINUE
      LUP = MIN0 (N, P)
      DO 210 L = 1, LUP
        IF (L.LT.PL .OR. L.GE.PU) GO TO 120
C
C ***      LOCATE COLUMN OF LARGEST NORM AND BRING IT INTO THE PIVOT POSITION
        MAXNRM = 0.0D0
        MAXJ = L
        DO 100 J = L, PU
          IF (QRAUXR(J).LE.MAXNRM) GO TO 90
          MAXNRM = QRAUXR(J)
          MAXJ = J
90        CONTINUE
100     CONTINUE
        IF (MAXJ.EQ.L) GO TO 110
        CALL WSWAP (N, XR(1,L), XI(1,L), 1, XR(1,MAXJ), XI(1,MAXJ), 1)
        QRAUXR(MAXJ) = QRAUXR(L)
        QRAUXI(MAXJ) = QRAUXI(L)
        WORKR(MAXJ) = WORKR(L)
        WORKI(MAXJ) = WORKI(L)
        JP = JPVT(MAXJ)
        JPVT(MAXJ) = JPVT(L)
        JPVT(L) = JP
110     CONTINUE
120     CONTINUE
        QRAUXR(L) = 0.0D0
        QRAUXI(L) = 0.0D0
        IF (L.EQ.N) GO TO 200
C
C ***      COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L
        NRMXLR = WNRM2 (N-L+1, XR(L,L), XI(L,L), 1)
        NRMXLI = 0.0D0
        IF (CABS1 (NRMXLR, NRMXLI).EQ.0.0D0) GO TO 190
        IF (CABS1 (XR(L,L), XI(L,L)).EQ.0.0D0) GO TO 130
        CALL WSIGN (NRMXLR, NRMXLI, XR(L,L), XI(L,L), NRMXLR, NRMXLI)
130     CONTINUE
        CALL WDIV (1.0D0, 0.0D0, NRMXLR, NRMXLI, TR, TI)
        CALL WSCAL (N-L+1, TR, TI, XR(L,L), XI(L,L), 1)
        XR(L,L) = FLOP (1.0D0+XR(L,L))
C
C ***      APPLY TRANSFORMATION TO THE REMAINING COLUMNS, UPDATING THE NORMS
        LP1 = L+1
        IF (P.LT.LP1) GO TO 180
        DO 170 J = LP1, P
          TR = -WDOTCR (N-L+1, XR(L,L), XI(L,L), 1, XR(L,J), XI(L,J), 1)
          TI = -WDOTCI (N-L+1, XR(L,L), XI(L,L), 1, XR(L,J), XI(L,J), 1)
          CALL WDIV (TR, TI, XR(L,L), XI(L,L), TR, TI)
          CALL WAXPY (N-L+1, TR, TI, XR(L,L), XI(L,L),
     .                1, XR(L,J), XI(L,J), 1)
          IF (J.LT.PL .OR. J.GT.PU) GO TO 160
          IF (CABS1 (QRAUXR(J), QRAUXI(J)).EQ.0.0D0) GO TO 160
          TT = 1.0D0-(PYTHAG (XR(L,J), XI(L,J))/QRAUXR(J))**2
          TT = DMAX1 (TT, 0.0D0)
          TR = FLOP (TT)
          TT = FLOP (1.0D0+0.05D0*TT*(QRAUXR(J)/WORKR(J))**2)
          IF (TT.EQ.1.0D0) GO TO 140
          QRAUXR(J) = QRAUXR(J)*DSQRT (TR)
          QRAUXI(J) = QRAUXI(J)*DSQRT (TR)
          GO TO 150
C
140       CONTINUE
          QRAUXR(J) = WNRM2 (N-L, XR(L+1,J), XI(L+1,J), 1)
          QRAUXI(J) = 0.0D0
          WORKR(J) = QRAUXR(J)
          WORKI(J) = QRAUXI(J)
150       CONTINUE
160       CONTINUE
170     CONTINUE
C
C ***      SAVE THE TRANSFORMATION
180     CONTINUE
        QRAUXR(L) = XR(L,L)
        QRAUXI(L) = XI(L,L)
        XR(L,L) = -NRMXLR
        XI(L,L) = -NRMXLI
190     CONTINUE
200     CONTINUE
210   CONTINUE
C
      RETURN
      END
