      SUBROUTINE SQRLS (A,LDA,M,N,TOL,KR,B,X,RSD,JPVT,QRAUX,WORK,
     *                  ITASK,IND)
C***BEGIN PROLOGUE  SQRLS
C***DATE WRITTEN   870911   (YYMMDD)
C***REVISION DATE  871016   (YYMMDD)
C***CATEGORY NO.  D9
C***KEYWORDS  LEAST SQUARES,OVERDETERMINED,LINEAR EQUATIONS
C***AUTHOR  STEPHEN NASH (GEORGE MASON UNIVERSITY)
C***PURPOSE
C***PURPOSE  SQRLS solves an overdetermined, underdetermined or singular
C            system of linear equations in least square sense.  The
C            solution is obtained using a QR factorization of the
C            M  by  N  coefficient matrix  A.
C***DESCRIPTION
C    From the book "Numerical Methods and Software"
C       by  D. Kahaner, C. Moler, S. Nash
C           Prentice Hall 1988
C
C     SQRLS IS USED TO SOLVE IN A LEAST SQUARES SENSE
C     OVERDETERMINED, UNDERDETERMINED AND SINGULAR LINEAR SYSTEMS .
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 MINIMIMZES THE SUM OF SQUARES (2-NORM)
C     OF THE RESIDUAL,  A*X - B .
C     THE NUMERICAL RANK OF A IS DETERMINED USING THE TOLERANCE TOL.
C
C     SQRLS USES THE LINPACK SUBROUTINE SQRDC TO COMPUTE THE QR
C     FACTORIZATION, WITH COLUMN PIVOTING, OF AN  M  BY  N  MATRIX  A .
C     FOR MORE INFORMATION, SEE CHAPTER 9 OF THE REFERENCE BELOW.
C
C
C     ON ENTRY
C
C        A     REAL (LDA,N) .
C              THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED.
C              IN A LEAST SQUARES DATA FITTING PROBLEM, A(I,J) IS THE
C              VALUE OF THE J-TH BASIS (MODEL) FUNCTION AT THE I-TH
C              DATA POINT.
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   REAL.
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        QRAUX REAL(N)
C        WORK  REAL(N)
C                    THREE AUXILIARY ARRAYS USED TO FACTOR THE MATRIX A.
C                    (NOT REQUIRED IF ITASK .GT. 1)
C
C        B     REAL(M)
C              THE RIGHT HAND SIDE OF THE LINEAR SYSTEM.
C              IN A LEAST SQUARES DATA FITTING PROBLEM B(I) CONTAINS THE
C              VALUE OF I-TH OBSERVATION.
C
C        ITASK INTEGER.
C              IF ITASK=1, THEN SQRLS FACTORS THE MATRIX A AND
C                          SOLVES THE LEAST SQUARES PROBLEM.
C              IF ITASK=2, THEN SQRLS ASSUMES THAT THE MATRIX A
C                          WAS FACTORED WITH AN EARLIER CALL TO
C                          SQRLS, AND ONLY SOLVES THE LEAST SQUARES
C                          PROBLEM.
C
C     ON RETURN
C
C        X     REAL(N) .
C              A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM.
C
C        RSD   REAL(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        SQRLS 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 SQRLS.
C
C***REFERENCES  DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979
C***ROUTINES CALLED  SQRANK,SQRLSS,XERROR
C***END PROLOGUE  SQRLS
      INTEGER   LDA, M, N, ITASK, JPVT(1), KR, IND
      REAL      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 SQRANK (A, LDA, M, N, TOL, KR, JPVT, QRAUX, WORK)
C
C SOLVE LEAST-SQUARES PROBLEM
C
      CALL SQRLSS (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, '(
     *  ''SQRLS 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, '(
     *  ''SQRLS 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, '(
     *  ''SQRLS ERROR (IND=-3) -- ITASK='', I5, '' IS LESS THAN 1.'')
     *               ') ITASK
      CALL XERROR(MSG(1:51), 51, -3, 0)
      RETURN
      END
