      SUBROUTINE L2SLV(M, N, M1, L, A, MM, B, MB, W, TOL, N1, IPIVOT,
     *                 X, NN, RES, MR, QR, MMPNN, C, IFAULT)
C ** PURPOSE **
C SUBROUTINE L2SLV COMPUTES LEAST SQUARES SOLUTIONS TO OVERDETERMINED
C AND UNDERDETERMINED SYSTEMS OF LINEAR EQUATIONS.  THE METHOD USED IS
C A MODIFIED GRAM-SCHMIDT ORTHOGONAL DECOMPOSITION WITH ITERATIVE
C REFINEMENT OF THE SOLUTION.  THE SOLUTION MAY BE SUBJECT TO LINEAR
C EQUALITY CONSTRAINTS.  OUTPUT INCLUDES THE LEAST SQUARES
C COEFFICIENTS, RESIDUALS, UNSCALED COVARIANCE MATRIX, AND INFORMATION
C ON THE BEHAVIOR OF THE ITERATIVE REFINEMENT PROCEDURE.
C MATRIX A IS THE GIVEN MATRIX OF A SYSTEM OF M LINEAR EQUATIONS IN N
C UNKNOWNS, AND MATRIX W IS A GIVEN DIAGONAL MATRIX OF WEIGHTS WITH ALL
C DIAGONAL ELEMENTS NONNEGATIVE.  LET H = W*A.
C IN THE EVENT THAT N1 (THE COMPUTED RANK OF MATRIX H) IS LESS THAN N
C (THE NUMBER OF UNKNOWN COEFFICIENTS), A UNIQUE SOLUTION VECTOR HAVING
C N ELEMENTS CAN BE OBTAINED BY IMPOSING THE CONDITION THAT THE
C SOLUTION BE OF MINIMAL EUCLIDEAN NORM.  SUCH A SOLUTION IS SOUGHT IN
C THE CASE OF UNDERDETERMINED OR RANK-DEFICIENT PROBLEMS.
C
C ** INPUT VARIABLES **
C M      TOTAL NUMBER OF EQUATIONS.
C N      NUMBER OF UNKNOWN COEFFICIENTS.
C M1     NUMBER OF LINEAR CONSTRAINTS (0.LE.M1.LE.M AND M1.LE.N).
C L      NUMBER OF RIGHT-HAND SIDES (VECTORS OF OBSERVATIONS).
C A      TWO-DIMENSIONAL ARRAY OF SIZE (MM,N).  ON ENTRY, THE ARRAY A
C        CONTAINS THE GIVEN MATRIX OF A SYSTEM OF M LINEAR EQUATIONS
C        IN N UNKNOWNS, WHERE THE FIRST M1 EQUATIONS ARE TO BE
C        SATISFIED EXACTLY.  A IS LEFT INTACT ON EXIT.
C B      TWO-DIMENSIONAL ARRAY OF SIZE (MB,L).  ON ENTRY, B CONTAINS
C        THE L GIVEN RIGHT-HAND SIDES (VECTORS OF OBSERVATIONS).  B IS
C        LEFT INTACT ON EXIT.
C W      VECTOR OF SIZE M.  ON ENTRY, W CONTAINS THE DIAGONAL ELEMENTS
C        OF A GIVEN DIAGONAL MATRIX OF WEIGHTS, ALL NONNEGATIVE.
C        (THE FIRST M1 ELEMENTS OF W ARE SET EQUAL TO 1.0 BY THE
C        PROGRAM WHEN M1 IS GREATER THAN ZERO.)
C TOL    PARAMETER USED IN DETERMINING THE RANK OF MATRIX H.
C        NOTE --
C        (1) IF TOL EQUALS ZERO, THE TOLERANCE USED IN SUBROUTINE
C            DECOM2 WILL BE BASED ON MACHINE PRECISION.
C        (2) IF TOL IS GREATER THAN ZERO, THIS VALUE OF TOL WILL BE
C            USED IN SETTING AN ABSOLUTE TOLERANCE FOR COMPARISON WITH
C            DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX OBTAINED IN
C            SUBROUTINE DECOM2.  THE VALUE OF TOL CAN BE BASED ON
C            KNOWLEDGE CONCERNING THE ACCURACY OF THE DATA.
C MM     DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN
C        THE ARRAY A.  MM MUST SATISFY MM.GE.M.
C MB     DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN
C        THE ARRAY B.  MB MUST SATISFY MB.GE.M.
C MR     DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN
C        THE ARRAY RES.  MR MUST SATISFY MR.GE.M.
C NN     DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN
C        THE ARRAY X.  NN MUST SATISFY NN.GE.N.
C MMPNN  DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN
C        THE ARRAY QR.  MMPNN MUST SATISFY MMPNN.GE.M+N.
C
C ** OUTPUT VARIABLES AND INTERNAL VARIABLES **
C N1     COMPUTED RANK OF MATRIX H, WHERE H = W*A.
C IPIVOT VECTOR OF SIZE N.  ON EXIT, THIS ARRAY RECORDS THE ORDER
C        IN WHICH THE COLUMNS OF H WERE SELECTED BY THE PIVOTING
C        SCHEME IN THE COURSE OF THE ORTHOGONAL DECOMPOSITION.
C        WHENEVER N1.LT.N, THE FIRST N1 ELEMENTS OF IPIVOT INDICATE
C        WHICH COLUMNS OF H WERE FOUND TO BE LINEARLY INDEPENDENT.
C X      TWO-DIMENSIONAL ARRAY OF SIZE (NN,L).  ON EXIT, X CONTAINS
C        THE SOLUTION VECTORS.
C RES    TWO-DIMENSIONAL ARRAY OF SIZE (MR,L).  ON EXIT, RES CONTAINS
C        THE RESIDUAL VECTORS.
C QR     TWO-DIMENSIONAL ARRAY OF SIZE (MMPNN,N).  ON EXIT, IF N1 = N
C        THEN QR CONTAINS THE UNSCALED COVARIANCE MATRIX. (QR IS USED
C        INTERNALLY TO STORE THE RESULTS FROM THE SUBROUTINE DECOM2.
C        THE RESULTS FROM DECOM2 ARE DESTROYED WHEN THE COVARIANCE
C        MATRIX IS COMPUTED.)
C C      VECTOR HAVING AT LEAST 6*(M+N)+2*L ELEMENTS USED (1) FOR
C        INTERNAL WORK SPACE AND (2) FOR RETURNING INFORMATION ON THE
C        BEHAVIOR OF THE ITERATIVE REFINEMENT PROCEDURE.
C        (A) NUMIT IS THE NUMBER OF ITERATIONS CARRIED OUT DURING THE
C            ITERATIVE REFINEMENT IN ATTEMPTING TO OBTAIN A SOLUTION
C            FOR THE K-TH RIGHT-HAND SIDE.
C            ON EXIT, C(K) = +NUMIT IF THE SOLUTION CONVERGED, AND
C                     C(K) = -NUMIT IF THE SOLUTION FAILED TO CONVERGE.
C        (B) DIGITX GIVES AN ESTIMATE OF THE NUMBER OF CORRECT DIGITS
C            IN THE INITIAL SOLUTION OF THE COEFFICIENTS FOR THE K-TH
C            RIGHT-HAND SIDE.  ON EXIT, C(K+L) = DIGITX.
C IFAULT FAULT INDICATOR WHICH IS ZERO IF NO ERRORS WERE ENCOUNTERED
C        AND POSITIVE IF ERRORS WERE DETECTED OR IF EVIDENCE OF SEVERE
C        ILL-CONDITIONING WAS FOUND.  IF IFAULT IS SET TO 1, 2, 3, 4,
C        5, 6 OR 7, EXECUTION IS TERMINATED.  EXECUTION CONTINUES WHEN
C        IFAULT IS SET EQUAL TO 8, 9 OR 10 PROVIDED THAT A SOLUTION
C        WAS OBTAINED FOR AT LEAST ONE RIGHT-HAND SIDE.  THE VALUE OF
C        IFAULT IS USED TO INDICATE THE FOLLOWING --
C        0 = NO ERRORS ENCOUNTERED.
C        1 = BAD INPUT PARAMETER (M, N OR L).
C        2 = BAD INPUT PARAMETER (M1).
C        3 = BAD DIMENSION.  EITHER M.GT.MM, M.GT.MB, M.GT.MR,
C            N.GT.NN, OR M+N.GT.MMPNN.
C        4 = AT LEAST ONE WEIGHT IS NEGATIVE.
C        5 = EITHER MATRIX H OR MATRIX OF CONSTRAINTS EQUALS ZERO.
C        6 = CONSTRAINTS ARE LINEARLY DEPENDENT.
C        7 = ALL SOLUTIONS FAILED TO CONVERGE.
C        8 = SOLUTION FAILED TO CONVERGE FOR AT LEAST ONE RIGHT-HAND
C            SIDE.
C        9 = LARGE NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE.
C       10 = ESTIMATED NUMBER OF DIGITS IN INITIAL SOLUTION OF
C            COEFFICIENTS IS SMALL.
C       11 = DIAGONAL ELEMENT OF COVARIANCE MATRIX WAS COMPUTED TO BE
C            NEGATIVE OWING TO ROUNDING ERROR. NO SEVERE CONDITIONING
C            PROBLEMS WERE DETECTED.
C       12 = DIAGONAL ELEMENT OF COVARIANCE MATRIX WAS COMPUTED TO BE
C            NEGATIVE OWING TO ROUNDING ERROR. THE PROBLEM APPEARS TO
C            BE EXTREMELY ILL-CONDITIONED.
C
C ** SUBROUTINES REQUIRED **
C SUBROUTINE DECOM2
C        USES MODIFIED GRAM-SCHMIDT ALGORITHM WITH PIVOTING TO
C        OBTAIN AN ORTHOGONAL DECOMPOSITION OF THE INPUT MATRIX.
C SUBROUTINE SOLVE2
C        COMPUTES COEFFICIENTS AND RESIDUALS.  ITERATIVE REFINEMENT IS
C        USED TO IMPROVE THE ACCURACY OF THE INITIAL SOLUTION.
C SUBROUTINE SOLVE3
C        CALLED ONLY BY SUBROUTINE SOLVE2.
C SUBROUTINE COVAR
C        COMPUTES UNSCALED COVARIANCE MATRIX OF THE COEFFICIENTS.
C
C ** STORAGE REQUIREMENTS **
C THE STORAGE REQUIRED FOR THE DIMENSIONED ARRAYS IN L2SLV IS
C   M*(2*N + 2*L + 7) + N*(N + L + 7) + 2*L
C LOCATIONS.  ALL ARRAYS REQUIRED IN SUBROUTINES CALLED BY L2SLV ARE
C DECLARED HEREIN AND ARE TRANSMITTED ONLY THROUGH PARAMETER LISTS OF
C CALL-SEQUENCES.
C
C ** PRECISION OF ARITHMETIC CALCULATIONS **
C SINGLE PRECISION ARITHMETIC IS USED FOR ALL CALCULATIONS EXCEPT THE
C DOUBLE PRECISION ACCUMULATION OF INNER PRODUCTS.  (THE VARIABLE SUM
C IS DECLARED TO BE DOUBLE PRECISION IN SUBROUTINES DECOM2, SOLVE2,
C SOLVE3 AND COVAR.)  IT IS ESSENTIAL FOR THE SUCCESS OF THE ITERATIVE
C REFINEMENT PROCEDURE THAT INNER PRODUCTS BE ACCUMULATED IN DOUBLE
C PRECISION.
C
C ** CONVERSION OF THE PROGRAM TO DOUBLE PRECISION **
C *********************************************************************
C * ON COMPUTERS HAVING SHORT WORD LENGTH (AS THE IBM 360/370) IT MAY *
C * BE DESIRABLE TO PERFORM ALL CALCULATIONS IN DOUBLE PRECISION.  IN *
C * THIS CASE, THE ITERATIVE REFINEMENT PRESENTLY INCLUDED IN SOLVE2  *
C * SHOULD BE OMITTED.                                                *
C * TO CONVERT THE PROGRAM TO DOUBLE PRECISION, THE FOLLOWING         *
C * APPROACH IS SUGGESTED.                                            *
C *                                                                   *
C * 1. VARIABLES PRESENTLY DECLARED TO BE REAL SHOULD BE DECLARED     *
C *    DOUBLE PRECISION.  THOSE TYPED INTEGER, DOUBLE PRECISION AND   *
C *    LOGICAL SHOULD NOT BE CHANGED.                                 *
C * 2. THE USE OF FAIL, NUMIT AND DIGITX SHOULD BE OMITTED.           *
C * 3. DESCRIPTION OF VARIABLE C (AT L2B 690-790) SHOULD READ --      *
C *    C  VECTOR HAVING AT LEAST 6*(M+N) ELEMENTS USED ONLY FOR       *
C *       INTERNAL WORK SPACE.                                        *
C * 4. THE VALUE OF ETA (AT L2B 1960) SHOULD BE SET SO THAT IT IS THE *
C *    SMALLEST POSITIVE DOUBLE PRECISION NUMBER SUCH THAT 1.0 + ETA  *
C *    IS GREATER THAN 1.0 IN DOUBLE PRECISION ARITHMETIC.            *
C *    FOR IBM COMPUTER TYPE, ETA = 16.**(-13)                        *
C *    FOR UNIVAC COMPUTER TYPE, ETA = 2.**(-59)                      *
C * 5. THE FOLLOWING FORTRAN FUNCTIONS SHOULD BE CHANGED --           *
C *      SINGLE PRECISION NAME     DOUBLE PRECISION NAME              *
C *             DBLE(X)                X                              *
C *             FLOAT(N)               DBLE(FLOAT(N))                 *
C *             SQRT(X)                DSQRT(X)                       *
C *    DBLE(X) IS USED IN SUBROUTINES DECOM2, SOLVE2, SOLVE3 AND      *
C *       COVAR.                                                      *
C *    FLOAT(N) IS USED IN SUBROUTINE DECOM2.                         *
C *    SQRT(X) IS USED IN SUBROUTINE L2SLV.                           *
C * 6. REPLACE STATEMENT L2B 2500 BY A STATEMENT READING              *
C *          K3 = 1                                                   *
C * 7. FURTHER DETAILS ARE GIVEN IN SUBROUTINE SOLVE2 IN CONNECTION   *
C *    WITH THE OMISSION OF ITERATIVE REFINEMENT.                     *
C * 8. IN SUBROUTINE L2SLV, STATEMENTS L2B 950-1000, 1820-1830, 2020, *
C *    2350-2360, 2480-2490, 3070, 3280-3570 AND 3590-3620 SHOULD BE  *
C *    OMITTED.                                                       *
C *    STATEMENT NUMBERS GIVEN HERE REFER TO THOSE IN THE RIGHT-HAND  *
C *    MARGIN.  CERTAIN COMMENTS IN SUBROUTINE L2SLV DO NOT APPLY TO  *
C *    THE DOUBLE PRECISION VERSION.                                  *
C *                                                                   *
C *********************************************************************
C
      INTEGER IPIVOT(N)
      REAL A(MM,N), B(MB,L), C(1), ETA, QR(MMPNN,N),
     * RES(MR,L), TOL, W(M), X(NN,L), Z
      REAL DIGITX
      LOGICAL FAIL
      LOGICAL SING
C
C SET VALUE OF ETA, A MACHINE-DEPENDENT PARAMETER.
C ETA, THE RELATIVE MACHINE PRECISION, IS THE SMALLEST POSITIVE REAL
C NUMBER SUCH THAT 1.0 + ETA IS GREATER THAN 1.0 IN FLOATING-POINT
C ARITHMETIC.
C
      ETA = SPMPAR(1)
C
C DEFAULT VALUE FOR TOL IS ZERO.
C
      IF (TOL.LT.0.0) TOL = 0.0
      IFAULT = 0
      KSUM = 0
C
C PERFORM INITIAL CHECKING OF INPUT PARAMETERS, DIMENSIONS AND
C WEIGHTS FOR POSSIBLE ERRORS.
C
      IF (M.GT.0 .AND. N.GT.0 .AND. L.GT.0) GO TO 10
      IFAULT = 1
      RETURN
   10 IF (M1.LE.M .AND. M1.LE.N .AND. M1.GE.0) GO TO 20
      IFAULT = 2
      RETURN
   20 IF (M.LE.MM .AND. M.LE.MB .AND. M.LE.MR .AND. N.LE.NN
     *            .AND. M+N.LE.MMPNN) GO TO 30
      IFAULT = 3
      RETURN
   30 DO 40 I=1,M
        IF (M1.GT.0 .AND. I.LE.M1) W(I) = 1.0
        IF (W(I).GE.0.0) GO TO 40
        IFAULT = 4
        RETURN
   40 CONTINUE
C
C SET PARAMETERS WHICH ALLOCATE VECTOR C TO CONTAIN CERTAIN FINAL
C RESULTS AND ALSO TO BE USED AS WORK SPACE.
C
C K1 IS STARTING POINT FOR NUMIT AND FAIL, OF LENGTH L.
C K2 IS STARTING POINT FOR DIGITX, OF LENGTH L.
C K3 IS STARTING POINT FOR D, OF LENGTH N.
C K4 IS STARTING POINT FOR K-TH COLUMN OF B, OF LENGTH M.
C K5 IS STARTING POINT FOR K-TH COLUMN OF X, OF LENGTH N.
C K6 IS STARTING POINT FOR K-TH COLUMN OF RES, OF LENGTH M.
C K7 IS STARTING POINT FOR WORK SPACE OF LENGTH M.
C K8 IS STARTING POINT FOR WORK SPACE OF LENGTH M.
C K9 IS STARTING POINT FOR WORK SPACE OF LENGTH N.
C K10 IS STARTING POINT FOR WORK SPACE OF LENGTH N.
C K11 IS STARTING POINT FOR WORK SPACE OF LENGTH M + N.
C K12 IS STARTING POINT FOR WORK SPACE OF LENGTH M + N.
C
      K1 = 1
      K2 = K1 + L
      K3 = K2 + L
      K4 = K3 + N
      K5 = K4 + M
      K6 = K5 + N
      K7 = K6 + M
      K8 = K7 + M
      K9 = K8 + M
      K10 = K9 + N
      K11 = K10 + N
      K12 = K11 + M + N
      K = K12 + M + N - 1
C
C MULTIPLY EACH ROW OF MATRIX A (M BY N) BY ITS APPROPRIATE WEIGHT AND
C STORE THE RESULT IN THE FIRST M ROWS OF ARRAY QR.  SET ARRAY C AND
C THE LAST N ROWS OF ARRAY QR EQUAL TO ZERO.
C
      DO 60 I=1,K
        C(I) = 0.0
   60 CONTINUE
      MP1 = M + 1
      MPN = M + N
      DO 90 J=1,N
        DO 70 I=1,M
          QR(I,J) = A(I,J)*W(I)
   70   CONTINUE
        DO 80 I=MP1,MPN
          QR(I,J) = 0.0
   80   CONTINUE
   90 CONTINUE
C
C OBTAIN AN ORTHOGONAL DECOMPOSITION OF THE MATRIX STORED IN THE FIRST
C M ROWS OF ARRAY QR AND COMPUTE ITS RANK.
C
      CALL DECOM2(M, N, M1, ETA, TOL, QR, C(K3), N1, IPIVOT, SING,
     * MMPNN)
C
      IF (.NOT.SING) GO TO 110
      IF (N1.GT.0) GO TO 100
      IFAULT = 5
      RETURN
  100 IFAULT = 6
      RETURN
C
C SEEK A SOLUTION (COEFFICIENTS AND RESIDUALS) FOR EACH OF THE L LEAST
C SQUARES PROBLEMS WHOSE RIGHT-HAND SIDES ARE GIVEN IN THE ARRAY B.
C
  110 ITER = -ALOG10(ETA)
      DO 200 K=1,L
C K-TH RIGHT-HAND SIDE.
        K0 = K4 - 1
        DO 120 I=1,M
          K0 = K0 + 1
          C(K0) = B(I,K)
  120   CONTINUE
C
        CALL SOLVE2(M, N, M1, A, C(K4), W, N1, IPIVOT, QR, C(K3),
     *   ETA, FAIL, NUMIT, DIGITX,
     *   C(K5), C(K6), C(K7), C(K8), C(K9), C(K10), C(K11), C(K12),
     *   MM, MMPNN)
C
        K0 = K5 - 1
        DO 130 J=1,N
          K0 = K0 + 1
          X(J,K) = C(K0)
  130   CONTINUE
        IF (M1.EQ.0) GO TO 150
        DO 140 I=1,M1
          RES(I,K) = 0.0
  140   CONTINUE
  150   M1P1 = M1 + 1
        IF (M1P1.GT.M) GO TO 170
        K0 = K6 + M1 - 1
        DO 160 I=M1P1,M
          K0 = K0 + 1
          RES(I,K) = C(K0)
  160   CONTINUE
  170   CONTINUE
C
C FOR RIGHT-HAND SIDES WHERE CONVERGENCE OF A SOLUTION IS REPORTED,
C A CHECK IS MADE FOR EVIDENCE OF SEVERE ILL-CONDITIONING.  SUCH
C EVIDENCE IS FURNISHED BY LARGE VALUES OF NUMIT (NUMBER OF ITERATIONS
C BEFORE CONVERGENCE WAS OBTAINED) AND SMALL VALUES OF DIGITX
C (ESTIMATE OF THE NUMBER OF CORRECT DIGITS IN THE INITIAL SOLUTION
C OF THE COEFFICIENTS).  IF NUMIT EXCEEDS -ALOG10(ETA) THEN IFAULT
C IS SET TO 9.  IF DIGITX IS LESS THAN 0.5 (HALF A DECIMAL DIGIT)
C THEN IFAULT IS SET TO 10.
C
        C(K) = FLOAT(NUMIT)
        IF (FAIL) C(K) = -C(K)
        K0 = K2 + K - 1
        C(K0) = DIGITX
        IF (.NOT.FAIL) GO TO 180
C KSUM IS A TALLY OF SOLUTIONS WHICH FAILED TO CONVERGE.
        KSUM = KSUM + 1
        IFAULT = 8
        GO TO 200
  180   IF (NUMIT.LE.ITER) GO TO 190
        IFAULT = 9
  190   IF (DIGITX.GE.0.5) GO TO 200
        IFAULT = 10
  200 CONTINUE
      IF (KSUM.LT.L) GO TO 210
      IFAULT = 7
      RETURN
  210 IF (N1.LT.N) RETURN
      DO 230 I=1,N
        MPI = M + I
        DO 220 J=1,N
          QR(I,J) = QR(MPI,J)
  220   CONTINUE
        QR(I,I) = 0.0
  230 CONTINUE
C
C COMPUTE THE UNSCALED COVARIANCE MATRIX OF THE COEFFICIENTS.
C
      CALL COVAR(N, M1, N1, IPIVOT, QR, C(K3), C(K9), MMPNN)
      IF (N.EQ.1) GO TO 260
      DO 250 J=2,N
      JM1 = J - 1
        DO 240 I=1,JM1
  240   QR(I,J) = QR(J,I)
  250 CONTINUE
C
C IN CERTAIN PROBLEMS, SOME DIAGONAL TERMS OF THE UNSCALED COVARIANCE
C MATRIX ARE EQUAL TO ZERO OR TO SMALL POSITIVE NUMBERS.  BECAUSE OF
C ROUNDING ERRORS, COMPUTED VALUES FOR THESE TERMS MAY BE SMALL
C NEGATIVE NUMBERS.  IFAULT IS SET TO 11 IF THIS OCCURS.
C
  260 DO 270 J=1,N
        IF (QR(J,J).LT.0.0) GO TO 280
  270 CONTINUE
      RETURN
  280 IF (IFAULT.NE.0) GO TO 290
      IFAULT = 11
      RETURN
  290 IFAULT = 12
      RETURN
      END
