      SUBROUTINE DQRLS (A,LDA,M,N,TOL,KR,B,X,RSD,JPVT,QRAUX,WORK,
     *                  ITASK,IND)
C***BEGIN PROLOGUE  DQRLS
C***REVISION SEPTEMBER 11, 1987
C***CATEGORY NO.  E2G1A
C***KEYWORD(S)  OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS
C***AUTHOR  STEPHEN NASH (GEORGE MASON UNIVERSITY)
C***DATE WRITTEN SEPTEMBER 11, 1987
C***PURPOSE
C      SOLVES AN OVERDETERMINED, UNDERDETERMINED, OR SINGULAR SYSTEM OF
C      LINEAR EQUATIONS IN LEAST SQUARE SENSE.  THE SOLUTION IS OBTAINED
C      USING A QR FACTORIZATION OF THE  M  BY  N  COEFFICIENT MATRIX  A.
C
C***DESCRIPTION
C     DQRLS IS USED TO SOLVE OVERDETERMINED, UNDERDETERMINED AND 
C     SINGULAR LINEAR SYSTEMS IN A LEAST SQUARES SENSE.
C     THE SYSTEM IS  A*X  APPROXIMATES  B  WHERE  A  IS  M  BY  N.
C     B  IS A GIVEN  M-VECTOR, AND  X  IS THE  N-VECTOR TO BE COMPUTED.  
C     A SOLUTION  X  IS FOUND WHICH MINIMIZES THE SUM OF SQUARES (2-NORM)
C     OF THE RESIDUAL,  A*X - B .
C     DQRLS USES THE LINPACK SUBROUTINE SQRDC TO COMPUTE THE QR
C     FACTORIZATION, WITH COLUMN PIVOTING, OF AN  M  BY  N  MATRIX  A .
C     THE NUMERICAL RANK IS DETERMINED USING THE TOLERANCE TOL.
C
C     FOR MORE INFORMATION, SEE CHAPTER 9 OF LINPACK USERS GUIDE,
C     J. DONGARRA ET ALL, SIAM, 1979.
C
C
C     ON ENTRY
C
C        A     DOUBLE PRECISION (LDA,N) . 
C              THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED.
C
C        LDA   INTEGER.
C              THE LEADING DIMENSION OF A .
C
C        M     INTEGER.
C              THE NUMBER OF ROWS OF A .
C
C        N     INTEGER.
C              THE NUMBER OF COLUMNS OF  A .
C
C        TOL   DOUBLE PRECISION.
C              A RELATIVE TOLERANCE USED TO DETERMINE THE NUMERICAL
C              RANK.  THE PROBLEM SHOULD BE SCALED SO THAT ALL THE 
C              ELEMENTS OF  A   HAVE ROUGHLY THE SAME ABSOLUTE ACCURACY
C              EPS.  THEN A REASONABLE VALUE FOR  TOL  IS ROUGHLY  EPS  
C              DIVIDED BY THE MAGNITUDE OF THE LARGEST ELEMENT.
C
C        JPVT  INTEGER(N)
C
C        QRAUX DOUBLE PRECISION(N)
C
C        WORK  DOUBLE PRECISION(N)
C              THREE AUXILIARY VECTORS USED TO FACTOR THE MATRIX A.
C              (ARRAY WORK IS NOT REQUIRED IF ITASK .GT. 1)
C
C        B     DOUBLE PRECISION(M)
C              THE RIGHT HAND SIDE OF THE LINEAR SYSTEM.
C
C        ITASK INTEGER.
C              IF ITASK=1, THEN DQRLS FACTORS THE MATRIX A AND
C                          SOLVES THE LEAST SQUARES PROBLEM.
C              IF ITASK>1, THEN DQRLS ASSUMES THAT THE MATRIX A
C                          WAS FACTORED WITH AN EARLIER CALL TO
C                          DQRLS, AND ONLY SOLVES THE LEAST SQUARES
C                          PROBLEM.
C
C     ON RETURN
C
C        X     DOUBLE PRECISION(N) .
C              A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM.
C
C        RSD   DOUBLE PRECISION(M) .
C              THE RESIDUAL, B - A*X .  RSD MAY OVERWRITE  B .
C
C        IND   INTEGER
C              ERROR CODE:  IND =  0:  NO ERROR
C                           IND = -1: N .GT. LDA   (FATAL ERROR)
C                           IND = -2: N .LT. 1     (FATAL ERROR)
C                           IND = -3: ITASK .LT. 1 (FATAL ERROR)
C
C        A     CONTAINS THE OUTPUT FROM SQRDC.
C              THE TRIANGULAR MATRIX  R  OF THE QR FACTORIZATION IS
C              CONTAINED IN THE UPPER TRIANGLE AND INFORMATION NEEDED 
C              TO RECOVER THE ORTHOGONAL MATRIX  Q  IS STORED BELOW
C              THE DIAGONAL IN  A  AND IN THE VECTOR QRAUX .
C
C        KR    INTEGER.
C              THE NUMERICAL RANK.
C
C        JPVT  THE PIVOT INFORMATION FROM SQRDC.
C
C     COLUMNS JPVT(1),...,JPVT(KR) OF THE ORIGINAL MATRIX ARE LINEARLY
C     INDEPENDENT TO WITHIN THE TOLERANCE TOL AND THE REMAINING COLUMNS
C     ARE LINEARLY DEPENDENT.  ABS(A(1,1))/ABS(A(KR,KR))  IS AN ESTIMATE
C     OF THE CONDITION NUMBER OF THE MATRIX OF INDEPENDENT COLUMNS,
C     AND OF  R .  THIS ESTIMATE WILL BE .LE. 1/TOL .
C
C      USAGE....
C        DQRLS CAN BE EFFICIENTLY USED TO SOLVE SEVERAL LEAST SQUARES
C      PROBLEMS WITH THE SAME MATRIX A.  THE FIRST SYSTEM IS SOLVED
C      WITH ITASK = 1.  THE SUBSEQUENT SYSTEMS ARE SOLVED WITH
C      ITASK = 2, TO AVOID THE RECOMPUTATION OF THE MATRIX FACTORS.
C      THE PARAMETERS  KR, JPVT, AND QRAUX MUST NOT BE MODIFIED
C      BETWEEN CALLS TO DQRLS.
C
C***REFERENCE(S)
C      DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979
C***ROUTINES CALLED  DQRANK, DQRLSS, XERROR
C***END PROLOGUE
      INTEGER  LDA, M, N, ITASK, JPVT(1), KR, IND
      DOUBLE PRECISION  A(LDA,1), B(1), X(1), RSD(1), QRAUX(1), 
     *     WORK(1), TOL
      CHARACTER  MSG*54
C***FIRST EXECUTABLE STATEMENT
      IF (LDA   .LT. N) GO TO 10
      IF (N     .LE. 0) GO TO 20
      IF (ITASK .LT. 1) GO TO 30
      IND = 0
C
      IF (ITASK .EQ. 1)
C
C FACTOR MATRIX
C
     *    CALL DQRANK (A, LDA, M, N, TOL, KR, JPVT, QRAUX, WORK)
C
C SOLVE LEAST-SQUARES PROBLEM
C
      CALL DQRLSS (A, LDA, M, N, KR, B, X, RSD, JPVT, QRAUX)
      RETURN
C
C ERROR IN CALLING SEQUENCE
C
C LDA .LT. N
10    IND = -1
      WRITE (MSG, '(
     *  ''DQRLS ERROR (IND=-1) -- LDA='', I5, '' IS LESS THAN N='',
     *      I5       )'   ) LDA, N
      CALL XERROR(MSG(1:54), 54, -1, 0)
      RETURN
C
C N .LT. 1
20    IND = -2
      WRITE (MSG, '(
     *  ''DQRLS ERROR (IND=-2) -- N='', I5, '' IS LESS THAN 1.'') ')N
      CALL XERROR(MSG(1:47), 47, -2, 0)
      RETURN
C
C ITASK .LT. 1
30    IND = -3
      WRITE (MSG, '(
     *  ''DQRLS ERROR (IND=-3) -- ITASK='', I5, '' IS LESS THAN 1.'')
     *               ') ITASK
      CALL XERROR(MSG(1:51), 51, -3, 0)
      RETURN
      END

