      SUBROUTINE LSI(W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, IP)
C     THIS IS A COMPANION SUBPROGRAM TO LSEI( ).
C     THE DOCUMENTATION FOR LSEI( ) HAS MORE COMPLETE
C     USAGE INSTRUCTIONS.
C     WRITTEN BY R. J. HANSON, SLA.
C
C     SOLVE..
C              AX = B,  A  MA BY N  (LEAST SQUARES EQUATIONS)
C     SUBJECT TO..
C
C              GX.GE.H, G  MG BY N  (INEQUALITY CONSTRAINTS)
C
C     INPUT..
C
C      W(*,*) CONTAINS  (A B) IN ROWS 1,...,MA+MG, COLS 1,...,N+1.
C                       (G H)
C
C     MDW,MA,MG,N
C              CONTAIN (RESP) VAR. DIMENSION OF W(*,*),
C              AND MATRIX DIMENSIONS.
C
C     PRGOPT(*),
C              PROGRAM OPTION VECTOR.
C
C     OUTPUT..
C
C      X(*),RNORM
C
C              SOLUTION VECTOR(UNLESS MODE=2), LENGTH OF AX-B.
C
C      MODE
C              =0   INEQUALITY CONSTRAINTS ARE COMPATIBLE.
C              =2   INEQUALITY CONSTRAINTS CONTRADICTORY.
C
C      WS(*),
C              WORKING STORAGE OF DIMENSION K+N+(MG+2)*(N+7),
C              WHERE K=MAX(MA+MG,N).
C      IP(MG+2*N+1)
C              INTEGER WORKING STORAGE
C
C     SUBROUTINES CALLED
C
C     LPDP          THIS SUBPROGRAM MINIMIZES A SUM OF SQUARES
C                   OF UNKNOWNS SUBJECT TO LINEAR INEQUALITY
C                   CONSTRAINTS.  PART OF THIS PACKAGE.
C
C     SDOT,SSCAL    SUBROUTINES FROM THE BLAS PACKAGE.
C     SAXPY,SASUM,  THESE ARE FOUND IN THE SANDIA MATH. LIBRARY,
C     SCOPY,SSWAP   AND ARE ALSO DOCUMENTED IN THE SANDIA TECH. REPT.,
C                   SAND77-0898.
C
C     HFTI          SOLVES AN UNCONSTRAINED LINEAR LEAST SQUARES
C                   PROBLEM.
C
C     H12           SUBROUTINE TO CONSTRUCT AND APPLY A HOUSEHOLDER
C                   TRANSFORMATION.
C
C     SPMPAR        FUNCTION TO COMPUTE THE RELATIVE MACHINE
C                   PRECISION.
C
C     REVISED 1984 BY A.H. MORRIS (NSWC)
C
      DIMENSION W(MDW,1), PRGOPT(1), X(N), IP(N), WS(1), RNRM(1)
      LOGICAL COV
      DATA ZERO/0.0/, HALF/0.5/, ONE/1.0/
C
      EPS = SPMPAR(1)
   30 MODE = 0
      M = MA + MG
      KRANK = 0
      NP1 = N + 1
      IF (N.LE.0 .OR. M.LE.0) GO TO 70
      GO TO 480
C
C     PROCESS-OPTION-VECTOR
C
C     COMPUTE MATRIX NORM OF LEAST SQUARES EQUAS.
   40 ANORM = ZERO
      DO 50 J=1,N
         ANORM = AMAX1(ANORM,SASUM(MA,W(1,J),1))
   50 CONTINUE
C
C     SET TOL FOR HFTI( ) RANK TEST.
      TAU = TOL*ANORM
C
C     COMPUTE HOUSEHOLDER ORTHOGONAL DECOMP OF MATRIX.
      CALL SCOPY(MA, W(1,NP1), 1, WS, 1)
      K = MAX0(M,N)
      MINMAN = MIN0(MA,N)
      N1 = K + 1
      N2 = N1 + N
      CALL HFTI(W, MDW, MA, N, WS, 1, 1, TAU, KRANK, RNRM, WS(N2),
     *             WS(N1), IP)
      RNORM = RNRM(1)
      FAC = ONE
      IF (KRANK.LT.MA) FAC = RNORM**2/FLOAT(MA-KRANK)
      GO TO 80
C
C     REDUCE-TO-LPDP-AND-SOLVE
   60 CONTINUE
   70 IP(1) = KRANK
      IP(2) = N + MAX0(M,N) + (MG+2)*(N+7)
      RETURN
   80 CONTINUE
C
C     TO REDUCE-TO-LPDP-AND-SOLVE
      MAP1 = MA + 1
C
C     COMPUTE INEQ. RT-HAND SIDE FOR LPDP.
      IF (.NOT.(MA.LT.M)) GO TO 240
      DO 90 I=MAP1,M
         W(I,NP1) = W(I,NP1) - SDOT(N,W(I,1),MDW,WS,1)
   90 CONTINUE
      DO 100 I=1,MINMAN
         J = IP(I)
C
C     APPLY PERMUTATIONS TO COLS OF INEQ. CONSTRAINT MATRIX.
         CALL SSWAP(MG, W(MAP1,I), 1, W(MAP1,J), 1)
  100 CONTINUE
C
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO CONSTRAINT MATRIX.
      IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 120
      DO 110 II=1,KRANK
         I = KRANK + 1 - II
         L = N1 + I
         CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), W(MAP1,1),
     *    MDW, 1, MG)
  110 CONTINUE
C
C     COMPUTE PERMUTED INEQ. CONSTR. MATRIX TIMES R-INVERSE.
  120 DO 150 I=MAP1,M
         IF (.NOT.(0.LT.KRANK)) GO TO 140
         DO 130 J=1,KRANK
            W(I,J) = (W(I,J)-SDOT(J-1,W(1,J),1,W(I,1),MDW))/W(J,J)
  130    CONTINUE
  140    CONTINUE
  150 CONTINUE
C
C     SOLVE THE REDUCED PROBLEM WITH LPDP ALGORITHM,
C     THE LEAST PROJECTED DISTANCE PROBLEM.
      CALL LPDP(W(MAP1,1), MDW, MG, KRANK, N-KRANK, PRGOPT, X, XNORM,
     * MDLPDP, WS(N2), IP(N+1))
      IF (.NOT.(MDLPDP.EQ.1)) GO TO 220
      IF (.NOT.(KRANK.GT.0)) GO TO 170
C
C     COMPUTE SOLN IN ORIGINAL COORDINATES.
      DO 160 II=1,KRANK
         I = KRANK + 1 - II
         X(I) = (X(I)-SDOT(II-1,W(I,I+1),MDW,X(I+1),1))/W(I,I)
  160 CONTINUE
C
C     APPLY HOUSEHOLDER TRANS. TO SOLN VECTOR.
  170 IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 190
      DO 180 I=1,KRANK
         L = N1 + I
         CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), X, 1, 1, 1)
  180 CONTINUE
C
C     REPERMUTE VARIABLES TO THEIR INPUT ORDER.
  190 DO 200 II=1,MINMAN
         I = MINMAN + 1 - II
         J = IP(I)
         CALL SSWAP(1, X(I), 1, X(J), 1)
  200 CONTINUE
C
C     VARIABLES ARE NOW IN ORIG. COORDINATES.
C     ADD SOLN OF UNSCONSTRAINED PROB.
      DO 210 I=1,N
         X(I) = X(I) + WS(I)
  210 CONTINUE
C
C     COMPUTE THE RESIDUAL VECTOR NORM.
      RNORM = SQRT(RNORM**2+XNORM**2)
      MODE = 0
      GO TO 230
  220 MODE = 2
  230 GO TO 250
  240 CALL SCOPY(N, WS, 1, X, 1)
  250 IF (.NOT.(COV .AND. KRANK.GT.0)) GO TO 470
C
C     COMPUTE COVARIANCE MATRIX BASED ON THE ORTHOGONAL DECOMP.
C     FROM HFTI( ).
C
      KRM1 = KRANK - 1
      KRP1 = KRANK + 1
C
C     COPY DIAG. TERMS TO WORKING ARRAY.
      CALL SCOPY(KRANK, W, MDW+1, WS(N2), 1)
C
C     RECIPROCATE DIAG. TERMS.
      DO 260 J=1,KRANK
         W(J,J) = ONE/W(J,J)
  260 CONTINUE
      IF (.NOT.(KRANK.GT.1)) GO TO 290
C
C     INVERT THE UPPER TRIANGULAR QR FACTOR ON ITSELF.
      DO 280 I=1,KRM1
         IP1 = I + 1
         DO 270 J=IP1,KRANK
            W(I,J) = -SDOT(J-I,W(I,I),MDW,W(I,J),1)*W(J,J)
  270    CONTINUE
  280 CONTINUE
C
C     COMPUTE THE INVERTED FACTOR TIMES ITS TRANSPOSE.
  290 DO 310 I=1,KRANK
         DO 300 J=I,KRANK
            W(I,J) = SDOT(KRANK+1-J,W(I,J),MDW,W(J,J),MDW)
  300    CONTINUE
  310 CONTINUE
      IF (.NOT.(KRANK.LT.N)) GO TO 430
C
C     ZERO OUT LOWER TRAPEZOIDAL PART.
C     COPY UPPER TRI. TO LOWER TRI. PART.
      DO 320 J=1,KRANK
         CALL SCOPY(J, W(1,J), 1, W(J,1), MDW)
  320 CONTINUE
      DO 330 I=KRP1,N
         W(I,1) = ZERO
         CALL SCOPY(I, W(I,1), 0, W(I,1), MDW)
  330 CONTINUE
C
C     APPLY RIGHT SIDE TRANSFORMATIONS TO LOWER TRI.
      N3 = N2 + KRP1
      DO 410 I=1,KRANK
         L = N1 + I
         K = N2 + I
         RB = WS(L-1)*WS(K-1)
         IF (.NOT.(RB.LT.ZERO)) GO TO 400
C
C     IF RB.GE.ZERO, TRANSFORMATION CAN BE REGARDED AS ZERO.
         RB = ONE/RB
C
C     STORE UNSCALED RANK-ONE HOUSEHOLDER UPDATE IN WORK ARRAY.
         WS(N3) = ZERO
         CALL SCOPY(N, WS(N3), 0, WS(N3), 1)
         L = N1 + I
         K = N3 + I
         WS(K-1) = WS(L-1)
         DO 340 J=KRP1,N
            K = N3 + I
            WS(K-1) = W(I,J)
  340    CONTINUE
         DO 350 J=1,N
            L = N3 + I
            K = N3 + J
            WS(J) = SDOT(J-I,W(J,I),MDW,WS(L-1),1) +
     *       SDOT(N-J+1,W(J,J),1,WS(K-1),1)
            WS(J) = WS(J)*RB
  350    CONTINUE
         L = N3 + I
         GAM = SDOT(N-I+1,WS(L-1),1,WS(I),1)*RB
         GAM = GAM*HALF
         CALL SAXPY(N-I+1, GAM, WS(L-1), 1, WS(I), 1)
         DO 390 J=I,N
            IF (.NOT.(I.GT.1)) GO TO 370
            IM1 = I - 1
            K = N3 + J
            DO 360 L=1,IM1
               W(J,L) = W(J,L) + WS(K-1)*WS(L)
  360       CONTINUE
  370       K = N3 + J
            DO 380 L=I,J
               IL = N3 + L
               W(J,L) = W(J,L) + WS(J)*WS(IL-1) + WS(L)*WS(K-1)
  380       CONTINUE
  390    CONTINUE
  400    CONTINUE
  410 CONTINUE
C
C     COPY LOWER TRI. TO UPPER TRI. TO SYMMETRIZE THE COVARIANCE MATRIX.
      DO 420 I=1,N
         CALL SCOPY(I, W(I,1), MDW, W(1,I), 1)
  420 CONTINUE
C
C     REPERMUTE ROWS AND COLS.
  430 DO 450 II=1,MINMAN
         I = MINMAN + 1 - II
         K = IP(I)
         IF (.NOT.(I.NE.K)) GO TO 440
         CALL SSWAP(1, W(I,I), 1, W(K,K), 1)
         CALL SSWAP(I-1, W(1,I), 1, W(1,K), 1)
         CALL SSWAP(K-I-1, W(I,I+1), MDW, W(I+1,K), 1)
         CALL SSWAP(N-K, W(I,K+1), MDW, W(K,K+1), MDW)
  440    CONTINUE
  450 CONTINUE
C
C     PUT IN NORMALIZED RESIDUAL SUM OF SQUARES SCALE FACTOR
C     AND SYMMETRIZE THE RESULTING COVARIANCE MARIX.
      DO 460 J=1,N
         CALL SSCAL(J, FAC, W(1,J), 1)
         CALL SCOPY(J, W(1,J), 1, W(J,1), MDW)
  460 CONTINUE
  470 GO TO 60
  480 CONTINUE
C
C     TO PROCESS-OPTION-VECTOR
C
C     THE NOMINAL TOLERANCE USED IN THE CODE,
      TOL = SQRT(EPS)
      COV = .FALSE.
      LAST = 1
      LINK = PRGOPT(1)
  490 IF (.NOT.(LINK.GT.1)) GO TO 40
      KEY = PRGOPT(LAST+1)
      IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO
      IF (KEY.EQ.5) TOL = AMAX1(EPS,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      LAST = LINK
      LINK = NEXT
      GO TO 490
      END
