      SUBROUTINE SQRANK(A,LDA,M,N,TOL,KR,JPVT,QRAUX,WORK)
C***BEGIN PROLOGUE  SQRANK
C***REVISION NOVEMBER 15, 1980
C***CATEGORY NO.  E2G1A
C***KEYWORD(S)  OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS
C***AUTHOR  D. KAHANER (NBS)
C***DATE WRITTEN
C***PURPOSE
C      COMPUTES THE QR FACTORIZATION OF AN  M  BY  N  MATRIX A.  THIS
C      ROUTINE IS USED IN CONJUNCTION WITH SQRLSS TO SOLVE LINEAR
C      SYSTEMS OF EQUATIONS IN A LEAST SQUARE SENSE.
C***DESCRIPTION
C
C     SQRANK IS USED IN CONJUNCTION WITH SQRLSS TO SOLVE
C     OVERDETERMINED, UNDERDETERMINED AND SINGULAR LINEAR SYSTEMS
C     IN A LEAST SQUARES SENSE.
C     SQRANK 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     ON ENTRY
C
C        A     REAL (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   REAL.
C              A RELATIVE TOLERANCE USED TO DETERMINE THE NUMERICAL
C              RANK.  THE PROBLEM SHOULD BE SCALED SO THAT ALL THE ELEMENTS
C              OF  A   HAVE ROUGHLY THE SAME ABSOLUTE ACCURACY, EPS.  THEN A
C              REASONABLE VALUE FOR  TOL  IS ROUGHLY  EPS  DIVIDED BY
C              THE MAGNITUDE OF THE LARGEST ELEMENT.
C
C        JPVT  INTEGER(N)
C
C        QRAUX REAL(N)
C
C        WORK  REAL(N)
C              THREE AUXILLIARY VECTORS USED BY SQRDC .
C
C     ON RETURN
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.....SEE SUBROUTINE SQRLSS
C
C***REFERENCE(S)
C      DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979
C***ROUTINES CALLED  SQRDC
C***END PROLOGUE
      INTEGER LDA,M,N,KR,JPVT(1)
      REAL A(LDA,1),TOL,QRAUX(1),WORK(1)
      INTEGER J,K
C***FIRST EXECUTABLE STATEMENT
      DO 10 J = 1, N
         JPVT(J) = 0
   10 CONTINUE
      CALL SQRDC(A,LDA,M,N,QRAUX,JPVT,WORK,1)
      KR = 0
      K = MIN0(M,N)
      DO 20 J = 1, K
         IF (ABS(A(J,J)) .LE. TOL*ABS(A(1,1))) GO TO 30
         KR = J
   20 CONTINUE
   30 RETURN
      END
      SUBROUTINE SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
C***BEGIN PROLOGUE  SQRDC
C***DATE WRITTEN   780814   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  D5
C***KEYWORDS  DECOMPOSITION,LINEAR ALGEBRA,LINPACK,MATRIX,
C             ORTHOGONAL TRIANGULAR
C***AUTHOR  STEWART, G. W., (U. OF MARYLAND)
C***PURPOSE  Uses Householder transformations to compute the QR
C            factorization of an N by P matrix X.  Column pivoting is
C            a users option.
C***DESCRIPTION
C
C     SQRDC uses Householder transformations to compute the QR
C     factorization of an N by P matrix X.  Column pivoting
C     based on the 2-norms of the reduced columns may be
C     performed at the user's option.
C
C     On Entry
C
C        X       REAL(LDX,P), where LDX .GE. N.
C                X contains the matrix whose decomposition is to be
C                computed.
C
C        LDX     INTEGER.
C                LDX is the leading dimension of the array X.
C
C        N       INTEGER.
C                N is the number of rows of the matrix X.
C
C        P       INTEGER.
C                P is the number of columns of the matrix X.
C
C        JPVT    INTEGER(P).
C                JPVT contains integers that control the selection
C                of the pivot columns.  The K-th column X(K) of X
C                is placed in one of three classes according to the
C                value of JPVT(K).
C
C                   If JPVT(K) .GT. 0, then X(K) is an initial
C                                      column.
C
C                   If JPVT(K) .EQ. 0, then X(K) is a free column.
C
C                   If JPVT(K) .LT. 0, then X(K) is a final column.
C
C                Before the decomposition is computed, initial columns
C                are moved to the beginning of the array X and final
C                columns to the end.  Both initial and final columns
C                are frozen in place during the computation and only
C                free columns are moved.  At the K-th stage of the
C                reduction, if X(K) is occupied by a free column,
C                it is interchanged with the free column of largest
C                reduced norm.  JPVT is not referenced if
C                JOB .EQ. 0.
C
C        WORK    REAL(P).
C                WORK is a work array.  WORK is not referenced if
C                JOB .EQ. 0.
C
C        JOB     INTEGER.
C                JOB is an integer that 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
C        X       X contains in its upper triangle the upper
C                triangular matrix R of the QR factorization.
C                Below its diagonal X contains information from
C                which the orthogonal part of the decomposition
C                can be recovered.  Note that if pivoting has
C                been requested, the decomposition is not that
C                of the original matrix X but that of X
C                with its columns permuted as described by JPVT.
C
C        QRAUX   REAL(P).
C                QRAUX contains further information required to recover
C                the orthogonal part of the decomposition.
C
C        JPVT    JPVT(K) contains the index of the column of the
C                original matrix that has been interchanged into
C                the K-th column, if pivoting was requested.
C
C     LINPACK.  This version dated 08/14/78 .
C     G. W. Stewart, University of Maryland, Argonne National Lab.
C
C     SQRDC uses the following functions and subprograms.
C
C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2
C     Fortran ABS,AMAX1,MIN0,SQRT
C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
C***ROUTINES CALLED  SAXPY,SDOT,SNRM2,SSCAL,SSWAP
C***END PROLOGUE  SQRDC
      INTEGER LDX,N,P,JOB
      INTEGER JPVT(1)
      REAL X(LDX,1),QRAUX(1),WORK(1)
C
      INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
      REAL MAXNRM,SNRM2,TT
      REAL SDOT,NRMXL,T
      LOGICAL NEGJ,SWAPJ
C
C***FIRST EXECUTABLE STATEMENT  SQRDC
      PL = 1
      PU = 0
      IF (JOB .EQ. 0) GO TO 60
C
C        PIVOTING HAS BEEN REQUESTED.  REARRANGE THE COLUMNS
C        ACCORDING TO JPVT.
C
         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 SSWAP(N,X(1,PL),1,X(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 SSWAP(N,X(1,PU),1,X(1,J),1)
                  JP = JPVT(PU)
                  JPVT(PU) = JPVT(J)
                  JPVT(J) = JP
   30          CONTINUE
               PU = PU - 1
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
C
C     COMPUTE THE NORMS OF THE FREE COLUMNS.
C
      IF (PU .LT. PL) GO TO 80
      DO 70 J = PL, PU
         QRAUX(J) = SNRM2(N,X(1,J),1)
         WORK(J) = QRAUX(J)
   70 CONTINUE
   80 CONTINUE
C
C     PERFORM THE HOUSEHOLDER REDUCTION OF X.
C
      LUP = MIN0(N,P)
      DO 200 L = 1, LUP
         IF (L .LT. PL .OR. L .GE. PU) GO TO 120
C
C           LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
C           INTO THE PIVOT POSITION.
C
            MAXNRM = 0.0E0
            MAXJ = L
            DO 100 J = L, PU
               IF (QRAUX(J) .LE. MAXNRM) GO TO 90
                  MAXNRM = QRAUX(J)
                  MAXJ = J
   90          CONTINUE
  100       CONTINUE
            IF (MAXJ .EQ. L) GO TO 110
               CALL SSWAP(N,X(1,L),1,X(1,MAXJ),1)
               QRAUX(MAXJ) = QRAUX(L)
               WORK(MAXJ) = WORK(L)
               JP = JPVT(MAXJ)
               JPVT(MAXJ) = JPVT(L)
               JPVT(L) = JP
  110       CONTINUE
  120    CONTINUE
         QRAUX(L) = 0.0E0
         IF (L .EQ. N) GO TO 190
C
C           COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
C
            NRMXL = SNRM2(N-L+1,X(L,L),1)
            IF (NRMXL .EQ. 0.0E0) GO TO 180
               IF (X(L,L) .NE. 0.0E0) NRMXL = SIGN(NRMXL,X(L,L))
               CALL SSCAL(N-L+1,1.0E0/NRMXL,X(L,L),1)
               X(L,L) = 1.0E0 + X(L,L)
C
C              APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
C              UPDATING THE NORMS.
C
               LP1 = L + 1
               IF (P .LT. LP1) GO TO 170
               DO 160 J = LP1, P
                  T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
                  CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
                  IF (J .LT. PL .OR. J .GT. PU) GO TO 150
                  IF (QRAUX(J) .EQ. 0.0E0) GO TO 150
                     TT = 1.0E0 - (ABS(X(L,J))/QRAUX(J))**2
                     TT = AMAX1(TT,0.0E0)
                     T = TT
                     TT = 1.0E0 + 0.05E0*TT*(QRAUX(J)/WORK(J))**2
                     IF (TT .EQ. 1.0E0) GO TO 130
                        QRAUX(J) = QRAUX(J)*SQRT(T)
                     GO TO 140
  130                CONTINUE
                        QRAUX(J) = SNRM2(N-L,X(L+1,J),1)
                        WORK(J) = QRAUX(J)
  140                CONTINUE
  150             CONTINUE
  160          CONTINUE
  170          CONTINUE
C
C              SAVE THE TRANSFORMATION.
C
               QRAUX(L) = X(L,L)
               X(L,L) = -NRMXL
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
      RETURN
      END
      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
      SUBROUTINE SQRLSS(A,LDA,M,N,KR,B,X,RSD,JPVT,QRAUX)
C***BEGIN PROLOGUE  SQRLSS
C***REVISION NOVEMBER 15, 1980
C***CATEGORY NO.  E2G1A
C***KEYWORD(S)  OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS
C***AUTHOR  D. KAHANER (NBS)
C***DATE WRITTEN
C***PURPOSE
C      SOLVES AN OVERDETERMINED, UNDERDETERMINED, OR SINGULAR SYSTEM OF
C      LINEAR EQUATIONS IN LEAST SQUARE SENSE.
C***DESCRIPTION
C
C     SQRLSS USES THE LINPACK SUBROUTINE SQRSL TO SOLVE AN OVERDETERMINED,
C     UNDERDETERMINED, OR SINGULAR LINEAR SYSTEM IN A LEAST SQUARES
C     SENSE.  IT MUST BE PRECEEDED BY SQRANK .
C     THE SYSTEM IS  A*X  APPROXIMATES  B  WHERE  A  IS  M  BY  N  WITH
C     RANK  KR  (DETERMINED BY SQRANK),  B  IS A GIVEN  M-VECTOR,
C     AND  X  IS THE  N-VECTOR TO BE COMPUTED.  A SOLUTION  X  WITH
C     AT MOST  KR  NONZERO COMPONENTS IS FOUND WHICH MINIMIMZES
C     THE 2-NORM OF THE RESIDUAL,  A*X - B .
C
C     ON ENTRY
C
C        A,LDA,M,N,KR,JPVT,QRAUX
C              THE OUTPUT FROM SQRANK .
C
C        B     REAL(M) .
C              THE RIGHT HAND SIDE OF THE LINEAR SYSTEM.
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 OVERWITE  B .
C
C      USAGE....
C        ONCE THE MATRIX A HAS BEEN FORMED, SQRANK SHOULD BE
C      CALLED ONCE TO DECOMPOSE IT.  THEN FOR EACH RIGHT HAND
C      SIDE, B, SQRLSS SHOULD BE CALLED ONCE TO OBTAIN THE
C      SOLUTION AND RESIDUAL.
C
C     INTEGER J,K,INFO
C     REAL T
C
C***REFERENCE(S)
C      DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979
C***ROUTINES CALLED  SQRSL
C***END PROLOGUE
      INTEGER LDA,M,N,KR,JPVT(1)
      REAL A(LDA,1),B(1),X(1),RSD(1),QRAUX(1)
C***FIRST EXECUTABLE STATEMENT
      IF (KR .NE. 0)
     1   CALL SQRSL(A,LDA,M,KR,QRAUX,B,RSD,RSD,X,RSD,RSD,110,INFO)
      DO 40 J = 1, N
         JPVT(J) = -JPVT(J)
         IF (J .GT. KR) X(J) = 0.0
   40 CONTINUE
      DO 70 J = 1, N
         IF (JPVT(J) .GT. 0) GO TO 70
         K = -JPVT(J)
         JPVT(J) = K
   50    CONTINUE
            IF (K .EQ. J) GO TO 60
            T = X(J)
            X(J) = X(K)
            X(K) = T
            JPVT(K) = -JPVT(K)
            K = JPVT(K)
         GO TO 50
   60    CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE SQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO)
C***BEGIN PROLOGUE  SQRSL
C***DATE WRITTEN   780814   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  D9,D2A1
C***KEYWORDS  LINEAR ALGEBRA,LINPACK,MATRIX,ORTHOGONAL TRIANGULAR,SOLVE
C***AUTHOR  STEWART, G. W., (U. OF MARYLAND)
C***PURPOSE  Applies the output of SQRDC to compute coordinate trans-
C            formations projections, and least squares solutions.
C***DESCRIPTION
C
C     SQRSL applies the output of SQRDC to compute coordinate
C     transformations, projections, and least squares solutions.
C     For K .LE. MIN(N,P), let XK be the matrix
C
C            XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))
C
C     formed from columnns JPVT(1), ... ,JPVT(K) of the original
C     N x P matrix X that was input to SQRDC (if no pivoting was
C     done, XK consists of the first K columns of X in their
C     original order).  SQRDC produces a factored orthogonal matrix Q
C     and an upper triangular matrix R such that
C
C              XK = Q * (R)
C                       (0)
C
C     This information is contained in coded form in the arrays
C     X and QRAUX.
C
C     On Entry
C
C        X      REAL(LDX,P)
C               X contains the output of SQRDC.
C
C        LDX    INTEGER
C               LDX is the leading dimension of the array X.
C
C        N      INTEGER
C               N is the number of rows of the matrix XK.  It must
C               have the same value as N in SQRDC.
C
C        K      INTEGER
C               K is the number of columns of the matrix XK.  K
C               must not be greater than MIN(N,P), where P is the
C               same as in the calling sequence to SQRDC.
C
C        QRAUX  REAL(P)
C               QRAUX contains the auxiliary output from SQRDC.
C
C        Y      REAL(N)
C               Y contains an N-vector that is to be manipulated
C               by SQRSL.
C
C        JOB    INTEGER
C               JOB specifies what is to be computed.  JOB has
C               the decimal expansion ABCDE, with the following
C               meaning.
C
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
C               Note that a request to compute B, RSD, or XB
C               automatically triggers the computation of QTY, for
C               which an array must be provided in the calling
C               sequence.
C
C     On Return
C
C        QY     REAL(N).
C               QY contains Q*Y, if its computation has been
C               requested.
C
C        QTY    REAL(N).
C               QTY contains TRANS(Q)*Y, if its computation has
C               been requested.  Here TRANS(Q) is the
C               transpose of the matrix Q.
C
C        B      REAL(K)
C               B contains the solution of the least squares problem
C
C                    minimize norm2(Y - XK*B),
C
C               if its computation has been requested.  (Note that
C               if pivoting was requested in SQRDC, the J-th
C               component of B will be associated with column JPVT(J)
C               of the original matrix X that was input into SQRDC.)
C
C        RSD    REAL(N).
C               RSD contains the least squares residual Y - XK*B,
C               if its computation has been requested.  RSD is
C               also the orthogonal projection of Y onto the
C               orthogonal complement of the column space of XK.
C
C        XB     REAL(N).
C               XB contains the least squares approximation XK*B,
C               if its computation has been requested.  XB is also
C               the orthogonal projection of Y onto the column space
C               of X.
C
C        INFO   INTEGER.
C               INFO is zero unless the computation of B has
C               been requested and R is exactly singular.  In
C               this case, INFO is the index of the first zero
C               diagonal element of R and B is left unaltered.
C
C     The parameters QY, QTY, B, RSD, and XB are not referenced
C     if their computation is not requested and in this case
C     can be replaced by dummy variables in the calling program.
C     To save storage, the user may in some cases use the same
C     array for different parameters in the calling sequence.  A
C     frequently occuring example is when one wishes to compute
C     any of B, RSD, or XB and does not need Y or QTY.  In this
C     case one may identify Y, QTY, and one of B, RSD, or XB, while
C     providing separate arrays for anything else that is to be
C     computed.  Thus the calling sequence
C
C          CALL SQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
C
C     will result in the computation of B and RSD, with RSD
C     overwriting Y.  More generally, each item in the following
C     list contains groups of permissible identifications for
C     a single callinng sequence.
C
C          1. (Y,QTY,B) (RSD) (XB) (QY)
C
C          2. (Y,QTY,RSD) (B) (XB) (QY)
C
C          3. (Y,QTY,XB) (B) (RSD) (QY)
C
C          4. (Y,QY) (QTY,B) (RSD) (XB)
C
C          5. (Y,QY) (QTY,RSD) (B) (XB)
C
C          6. (Y,QY) (QTY,XB) (B) (RSD)
C
C     In any group the value returned in the array allocated to
C     the group corresponds to the last member of the group.
C
C     LINPACK.  This version dated 08/14/78 .
C     G. W. Stewart, University of Maryland, Argonne National Lab.
C
C     SQRSL uses the following functions and subprograms.
C
C     BLAS SAXPY,SCOPY,SDOT
C     Fortran ABS,MIN0,MOD
C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
C***ROUTINES CALLED  SAXPY,SCOPY,SDOT
C***END PROLOGUE  SQRSL
      INTEGER LDX,N,K,JOB,INFO
      REAL X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1)
C
      INTEGER I,J,JJ,JU,KP1
      REAL SDOT,T,TEMP
      LOGICAL CB,CQY,CQTY,CR,CXB
C
C     SET INFO FLAG.
C
C***FIRST EXECUTABLE STATEMENT  SQRSL
      INFO = 0
C
C     DETERMINE WHAT IS TO BE COMPUTED.
C
      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.
C
      IF (JU .NE. 0) GO TO 40
         IF (CQY) QY(1) = Y(1)
         IF (CQTY) QTY(1) = Y(1)
         IF (CXB) XB(1) = Y(1)
         IF (.NOT.CB) GO TO 30
            IF (X(1,1) .NE. 0.0E0) GO TO 10
               INFO = 1
            GO TO 20
   10       CONTINUE
               B(1) = Y(1)/X(1,1)
   20       CONTINUE
   30    CONTINUE
         IF (CR) RSD(1) = 0.0E0
      GO TO 250
   40 CONTINUE
C
C        SET UP TO COMPUTE QY OR QTY.
C
         IF (CQY) CALL SCOPY(N,Y,1,QY,1)
         IF (CQTY) CALL SCOPY(N,Y,1,QTY,1)
         IF (.NOT.CQY) GO TO 70
C
C           COMPUTE QY.
C
            DO 60 JJ = 1, JU
               J = JU - JJ + 1
               IF (QRAUX(J) .EQ. 0.0E0) GO TO 50
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  T = -SDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J)
                  CALL SAXPY(N-J+1,T,X(J,J),1,QY(J),1)
                  X(J,J) = TEMP
   50          CONTINUE
   60       CONTINUE
   70    CONTINUE
         IF (.NOT.CQTY) GO TO 100
C
C           COMPUTE TRANS(Q)*Y.
C
            DO 90 J = 1, JU
               IF (QRAUX(J) .EQ. 0.0E0) GO TO 80
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  T = -SDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J)
                  CALL SAXPY(N-J+1,T,X(J,J),1,QTY(J),1)
                  X(J,J) = TEMP
   80          CONTINUE
   90       CONTINUE
  100    CONTINUE
C
C        SET UP TO COMPUTE B, RSD, OR XB.
C
         IF (CB) CALL SCOPY(K,QTY,1,B,1)
         KP1 = K + 1
         IF (CXB) CALL SCOPY(K,QTY,1,XB,1)
         IF (CR .AND. K .LT. N) CALL SCOPY(N-K,QTY(KP1),1,RSD(KP1),1)
         IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120
            DO 110 I = KP1, N
               XB(I) = 0.0E0
  110       CONTINUE
  120    CONTINUE
         IF (.NOT.CR) GO TO 140
            DO 130 I = 1, K
               RSD(I) = 0.0E0
  130       CONTINUE
  140    CONTINUE
         IF (.NOT.CB) GO TO 190
C
C           COMPUTE B.
C
            DO 170 JJ = 1, K
               J = K - JJ + 1
               IF (X(J,J) .NE. 0.0E0) GO TO 150
                  INFO = J
C           ......EXIT
                  GO TO 180
  150          CONTINUE
               B(J) = B(J)/X(J,J)
               IF (J .EQ. 1) GO TO 160
                  T = -B(J)
                  CALL SAXPY(J-1,T,X(1,J),1,B,1)
  160          CONTINUE
  170       CONTINUE
  180       CONTINUE
  190    CONTINUE
         IF (.NOT.CR .AND. .NOT.CXB) GO TO 240
C
C           COMPUTE RSD OR XB AS REQUIRED.
C
            DO 230 JJ = 1, JU
               J = JU - JJ + 1
               IF (QRAUX(J) .EQ. 0.0E0) GO TO 220
                  TEMP = X(J,J)
                  X(J,J) = QRAUX(J)
                  IF (.NOT.CR) GO TO 200
                     T = -SDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J)
                     CALL SAXPY(N-J+1,T,X(J,J),1,RSD(J),1)
  200             CONTINUE
                  IF (.NOT.CXB) GO TO 210
                     T = -SDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J)
                     CALL SAXPY(N-J+1,T,X(J,J),1,XB(J),1)
  210             CONTINUE
                  X(J,J) = TEMP
  220          CONTINUE
  230       CONTINUE
  240    CONTINUE
  250 CONTINUE
      RETURN
      END
