      SUBROUTINE WQRSL (XR, XI, LDX, N, K, QRAUXR, QRAUXI, YR, YI,
     .                  QYR, QYI, QTYR, QTYI, BR, BI, RSDR, RSDI,
     .                  XBR, XBI, JOB, INFO)
      IMPLICIT NONE
C
C APPLIES THE OUTPUT OF WQRDC TO COMPUTE COORDINATE TRANSFORMATIONS,
C  PROJECTIONS, AND LEAST SQUARES SOLUTIONS.
C  FOR K .LE. MIN(N,P), LET XK BE THE MATRIX:
C      XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))
C  FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL N X P MATRIX X
C  THAT WAS INPUT TO WQRDC (IF NO PIVOTING WAS DONE, XK CONSISTS OF THE FIRST
C  K COLUMNS OF X IN THEIR ORIGINAL ORDER).  WQRDC PRODUCES A FACTORED UNITARY
C  MATRIX Q AND AN UPPER TRIANGULAR MATRIX R SUCH THAT
C      XK = Q * (R)
C  (0)
C  THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS X AND QRAUX.
C
C ON ENTRY:
C    X      DOUBLE-COMPLEX(LDX,P),  CONTAINS THE OUTPUT OF WQRDC
C    LDX    INTEGER,  LEADING DIMENSION OF THE ARRAY X
C    N      INTEGER,  NUMBER OF ROWS OF THE MATRIX XK.  IT MUST
C            HAVE THE SAME VALUE AS N IN WQRDC
C    K      INTEGER,  NUMBER OF COLUMNS OF THE MATRIX XK.  K MUST NOT BE
C            GREATER THAN MIN(N,P), WHERE P IS THE SAME AS IN THE CALLING
C            SEQUENCE TO WQRDC
C    QRAUX  DOUBLE-COMPLEX(P),  CONTAINS THE AUXILIARY OUTPUT FROM WQRDC
C    Y      DOUBLE-COMPLEX(N),  CONTAINS AN N-VECTOR TO BE MANIPULATED BY WQRSL
C    JOB    INTEGER,  SPECIFIES WHAT IS TO BE COMPUTED.  HAS DECIMAL EXPANSION
C            ABCDE, WITH THE FOLLOWING MEANING:
C            IF A.NE.0, COMPUTE QY
C            IF B,C,D, OR E .NE. 0, COMPUTE QTY
C            IF C.NE.0, COMPUTE B
C            IF D.NE.0, COMPUTE RSD
C            IF E.NE.0, COMPUTE XB
C            NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB AUTOMATICALLY
C             TRIGGERS THE COMPUTATION OF QTY, FOR WHICH AN ARRAY MUST BE
C             PROVIDED IN THE CALLING SEQUENCE
C
C ON RETURN:
C    QY    DOUBLE-COMPLEX(N),  CONTAINS Q*Y, IF ITS COMPUTATION HAS BEEN
C               REQUESTED
C    QTY   DOUBLE-COMPLEX(N),  CONTAINS CTRANS(Q)*Y, IF ITS COMPUTATION HAS
C           BEEN REQUESTED.  HERE CTRANS(Q) IS THE CONJUGATE TRANSPOSE OF THE
C           MATRIX Q
C    B     DOUBLE-COMPLEX(K),  CONTAINS THE SOLUTION OF THE LEAST SQUARES
C           PROBLEM:
C              MINIMIZE NORM2(Y - XK*B),
C           IF ITS COMPUTATION HAS BEEN REQUESTED.  (NOTE THAT IF PIVOTING WAS
C           REQUESTED IN WQRDC, THE J-TH COMPONENT OF B WILL BE ASSOCIATED
C           WITH COLUMN JPVT(J) OF ORIGINAL MATRIX X INPUT INTO WQRDC.)
C    RSD   DOUBLE-COMPLEX(N),  CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B,
C           IF ITS COMPUTATION HAS BEEN REQUESTED.  RSD IS ALSO THE ORTHOGONAL
C           PROJECTION OF Y ONTO ORTHOGONAL COMPLEMENT OF COLUMN SPACE OF XK
C    XB    DOUBLE-COMPLEX(N),  CONTAINS THE LEAST SQUARES APPROXIMATION XK*B,
C           IF ITS COMPUTATION HAS BEEN REQUESTED.  XB IS ALSO THE ORTHOGONAL
C           PROJECTION OF Y ONTO THE COLUMN SPACE OF X
C    INFO  INTEGER,  ZERO UNLESS THE COMPUTATION OF B HAS BEEN REQUESTED AND
C           R IS EXACTLY SINGULAR.  IN THIS CASE, INFO IS THE INDEX OF THE
C           FIRST ZERO DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED
C
C THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED IF THEIR
C  COMPUTATION IS NOT REQUESTED AND IN THIS CASE CAN BE REPLACED BY DUMMY
C  VARIABLES IN THE CALLING PROGRAM.  TO SAVE STORAGE, THE USER MAY IN SOME
C  CASES USE THE SAME ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE.
C  A FREQUENTLY OCCURRING EXAMPLE IS WHEN ONE WISHES TO COMPUTE ANY OF B, RSD,
C  OR XB AND DOES NOT NEED Y OR QTY.  IN THIS CASE ONE MAY IDENTIFY Y, QTY,
C  AND ONE OF B, RSD, OR XB, WHILE PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE
C  THAT IS TO BE COMPUTED.  THUS THE CALLING SEQUENCE:
C     CALL WQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
C  WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD OVERWRITING Y.
C  MORE GENERALLY, EACH ITEM IN THE FOLLOWING LIST CONTAINS GROUPS OF
C  PERMISSIBLE IDENTIFICATIONS FOR A SINGLE CALLINNG SEQUENCE:
C     1. (Y,QTY,B) (RSD) (XB) (QY)
C     2. (Y,QTY,RSD) (B) (XB) (QY)
C     3. (Y,QTY,XB) (B) (RSD) (QY)
C     4. (Y,QY) (QTY,B) (RSD) (XB)
C     5. (Y,QY) (QTY,RSD) (B) (XB)
C     6. (Y,QY) (QTY,XB) (B) (RSD)
C  IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO THE GROUP
C  CORRESPONDS TO THE LAST MEMBER OF THE GROUP.
C
C VERSION 07/03/79 G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
      INTEGER LDX, N, K, JOB, INFO
      DOUBLE PRECISION XR(LDX,*), XI(LDX,*), QRAUXR(*), QRAUXI(*),
     .                 YR(*), YI(*), QYR(*), QYI(*), QTYR(*), QTYI(*),
     .                 BR(*), BI(*), RSDR(*), RSDI(*), XBR(*), XBI(*)
C
      INTEGER I, J, JJ, JU, KP1
      DOUBLE PRECISION TR, TI, TEMPR, TEMPI
      LOGICAL CB, CQY, CQTY, CR, CXB
      DOUBLE PRECISION ZDUMR, ZDUMI
C
      DOUBLE PRECISION CABS1, WDOTCR, WDOTCI
C
C
      CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
C
C
C ***      SET INFO FLAG
      INFO = 0
C
C ***      DETERMINE WHAT IS TO BE COMPUTED
      CQY = JOB/10000.NE.0
      CQTY = MOD (JOB, 10000).NE.0
      CB = MOD (JOB, 1000)/100.NE.0
      CR = MOD (JOB, 100)/10.NE.0
      CXB = MOD (JOB, 10).NE.0
      JU = MIN0 (K, N-1)
C
C ***      SPECIAL ACTION WHEN N = 1
      IF (JU.NE.0) GO TO 80
      IF (.NOT.CQY) GO TO 10
      QYR(1) = YR(1)
      QYI(1) = YI(1)
10    CONTINUE
      IF (.NOT.CQTY) GO TO 20
      QTYR(1) = YR(1)
      QTYI(1) = YI(1)
20    CONTINUE
      IF (.NOT.CXB) GO TO 30
      XBR(1) = YR(1)
      XBI(1) = YI(1)
30    CONTINUE
      IF (.NOT.CB) GO TO 60
      IF (CABS1 (XR(1,1), XI(1,1)).NE.0.0D0) GO TO 40
      INFO = 1
      GO TO 50
C
40    CONTINUE
      CALL WDIV (YR(1), YI(1), XR(1,1), XI(1,1), BR(1), BI(1))
50    CONTINUE
60    CONTINUE
      IF (.NOT.CR) GO TO 70
      RSDR(1) = 0.0D0
      RSDI(1) = 0.0D0
70    CONTINUE
      GO TO 290
C
C ***      SET UP TO COMPUTE QY OR QTY
80    CONTINUE
      IF (CQY) CALL WCOPY (N, YR, YI, 1, QYR, QYI, 1)
      IF (CQTY) CALL WCOPY (N, YR, YI, 1, QTYR, QTYI, 1)
      IF (.NOT.CQY) GO TO 110
C
C ***      COMPUTE QY
      DO 100 JJ = 1, JU
        J = JU-JJ+1
        IF (CABS1 (QRAUXR(J), QRAUXI(J)).EQ.0.0D0) GO TO 90
        TEMPR = XR(J,J)
        TEMPI = XI(J,J)
        XR(J,J) = QRAUXR(J)
        XI(J,J) = QRAUXI(J)
        TR = -WDOTCR (N-J+1, XR(J,J), XI(J,J), 1, QYR(J), QYI(J), 1)
        TI = -WDOTCI (N-J+1, XR(J,J), XI(J,J), 1, QYR(J), QYI(J), 1)
        CALL WDIV (TR, TI, XR(J,J), XI(J,J), TR, TI)
        CALL WAXPY (N-J+1, TR, TI, XR(J,J), XI(J,J),
     .              1, QYR(J), QYI(J), 1)
        XR(J,J) = TEMPR
        XI(J,J) = TEMPI
90      CONTINUE
100   CONTINUE
110   CONTINUE
      IF (.NOT.CQTY) GO TO 140
C
C ***      COMPUTE CTRANS(Q)*Y
      DO 130 J = 1, JU
        IF (CABS1 (QRAUXR(J), QRAUXI(J)).EQ.0.0D0) GO TO 120
        TEMPR = XR(J,J)
        TEMPI = XI(J,J)
        XR(J,J) = QRAUXR(J)
        XI(J,J) = QRAUXI(J)
        TR = -WDOTCR (N-J+1, XR(J,J), XI(J,J), 1, QTYR(J), QTYI(J), 1)
        TI = -WDOTCI (N-J+1, XR(J,J), XI(J,J), 1, QTYR(J), QTYI(J), 1)
        CALL WDIV (TR, TI, XR(J,J), XI(J,J), TR, TI)
        CALL WAXPY (N-J+1, TR, TI, XR(J,J), XI(J,J),
     .              1, QTYR(J), QTYI(J), 1)
        XR(J,J) = TEMPR
        XI(J,J) = TEMPI
120     CONTINUE
130   CONTINUE
C
C ***      SET UP TO COMPUTE B, RSD, OR XB
140   CONTINUE
      IF (CB) CALL WCOPY (K, QTYR, QTYI, 1, BR, BI, 1)
      KP1 = K+1
      IF (CXB) CALL WCOPY (K, QTYR, QTYI, 1, XBR, XBI, 1)
      IF (CR .AND. K.LT.N) CALL WCOPY (N-K, QTYR(KP1), QTYI(KP1),
     .                                 1, RSDR(KP1), RSDI(KP1), 1)
      IF (.NOT.CXB .OR. KP1.GT.N) GO TO 160
      DO 150 I = KP1, N
        XBR(I) = 0.0D0
        XBI(I) = 0.0D0
150   CONTINUE
160   CONTINUE
      IF (.NOT.CR) GO TO 180
      DO 170 I = 1, K
        RSDR(I) = 0.0D0
        RSDI(I) = 0.0D0
170   CONTINUE
180   CONTINUE
      IF (.NOT.CB) GO TO 230
C
C ***      COMPUTE B
      DO 210 JJ = 1, K
        J = K-JJ+1
        IF (CABS1 (XR(J,J), XI(J,J)).NE.0.0D0) GO TO 190
        INFO = J
        GO TO 220
190     CONTINUE
        CALL WDIV (BR(J), BI(J), XR(J,J), XI(J,J), BR(J), BI(J))
        IF (J.EQ.1) GO TO 200
        TR = -BR(J)
        TI = -BI(J)
        CALL WAXPY (J-1, TR, TI, XR(1,J), XI(1,J), 1, BR, BI, 1)
200     CONTINUE
210   CONTINUE
220   CONTINUE
230   CONTINUE
      IF (.NOT.CR .AND. .NOT.CXB) GO TO 280
C
C ***      COMPUTE RSD OR XB AS REQUIRED
      DO 270 JJ = 1, JU
        J = JU-JJ+1
        IF (CABS1 (QRAUXR(J), QRAUXI(J)).EQ.0.0D0) GO TO 260
        TEMPR = XR(J,J)
        TEMPI = XI(J,J)
        XR(J,J) = QRAUXR(J)
        XI(J,J) = QRAUXI(J)
        IF (.NOT.CR) GO TO 240
        TR = -WDOTCR (N-J+1, XR(J,J), XI(J,J), 1, RSDR(J), RSDI(J), 1)
        TI = -WDOTCI (N-J+1, XR(J,J), XI(J,J), 1, RSDR(J), RSDI(J), 1)
        CALL WDIV (TR, TI, XR(J,J), XI(J,J), TR, TI)
        CALL WAXPY (N-J+1, TR, TI, XR(J,J), XI(J,J),
     .              1, RSDR(J), RSDI(J), 1)
240     CONTINUE
        IF (.NOT.CXB) GO TO 250
        TR = -WDOTCR (N-J+1, XR(J,J), XI(J,J), 1, XBR(J), XBI(J), 1)
        TI = -WDOTCI (N-J+1, XR(J,J), XI(J,J), 1, XBR(J), XBI(J), 1)
        CALL WDIV (TR, TI, XR(J,J), XI(J,J), TR, TI)
        CALL WAXPY (N-J+1, TR, TI, XR(J,J), XI(J,J),
     .              1, XBR(J), XBI(J), 1)
250     CONTINUE
        XR(J,J) = TEMPR
        XI(J,J) = TEMPI
260     CONTINUE
270   CONTINUE
C
280   CONTINUE
290   CONTINUE
      RETURN
      END
