SUBROUTINE LSI(E,F,G,H,LE,ME,LG,MG,N,X,XNORM,W,JW,MODE)

C     FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF
C     INEQUALITY CONSTRAINED LINEAR LEAST SQUARES PROBLEM:

C                    MIN ||E*X-F||
C                     X

C                    S.T.  G*X >= H

C     THE ALGORITHM IS BASED ON QR DECOMPOSITION AS DESCRIBED IN
C     CHAPTER 23.5 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS

C     THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM
C     ARE NECESSARY
C     DIM(E) :   FORMAL (LE,N),    ACTUAL (ME,N)
C     DIM(F) :   FORMAL (LE  ),    ACTUAL (ME  )
C     DIM(G) :   FORMAL (LG,N),    ACTUAL (MG,N)
C     DIM(H) :   FORMAL (LG  ),    ACTUAL (MG  )
C     DIM(X) :   N
C     DIM(W) :   (N+1)*(MG+2) + 2*MG
C     DIM(JW):   LG
C     ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS E, F, G, AND H.
C     ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE.
C     X     STORES THE SOLUTION VECTOR
C     XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM
C     W     STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST
C           MG ELEMENTS
C     MODE  IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS:
C          MODE=1: SUCCESSFUL COMPUTATION
C               2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1)
C               3: ITERATION COUNT EXCEEDED BY NNLS
C               4: INEQUALITY CONSTRAINTS INCOMPATIBLE
C               5: MATRIX E IS NOT OF FULL RANK

C     03.01.1980, DIETER KRAFT: CODED
C     20.03.1987, DIETER KRAFT: REVISED TO FORTRAN 77

      INTEGER          I,J,LE,LG,ME,MG,MODE,N,JW(LG)
      DOUBLE PRECISION E(LE,N),F(LE),G(LG,N),H(LG),X(N),W(*),
     .                 DDOT,XNORM,DNRM2,EPMACH,T,ONE
      DATA             EPMACH/2.22D-16/,ONE/1.0D+00/

C  QR-FACTORS OF E AND APPLICATION TO F

      DO 10 I=1,N
      J=MIN(I+1,N)
      CALL H12(1,I,I+1,ME,E(1,I),1,T,E(1,J),1,LE,N-I)
   10 CALL H12(2,I,I+1,ME,E(1,I),1,T,F     ,1,1 ,1  )

C  TRANSFORM G AND H TO GET LEAST DISTANCE PROBLEM

      MODE=5
      DO 30 I=1,MG
          DO 20 J=1,N
              IF (ABS(E(J,J)).LT.EPMACH) GOTO 50
   20         G(I,J)=(G(I,J)-DDOT(J-1,G(I,1),LG,E(1,J),1))/E(J,J)
   30     H(I)=H(I)-DDOT(N,G(I,1),LG,F,1)

C  SOLVE LEAST DISTANCE PROBLEM

      CALL LDP(G,LG,MG,N,H,X,XNORM,W,JW,MODE)
      IF (MODE.NE.1)                     GOTO 50

C  SOLUTION OF ORIGINAL PROBLEM

      CALL DAXPY(N,ONE,F,1,X,1)
      DO 40 I=N,1,-1
          J=MIN(I+1,N)
   40     X(I)=(X(I)-DDOT(N-I,E(I,J),LE,X(J),1))/E(I,I)
      J=MIN(N+1,ME)
      T=DNRM2(ME-N,F(J),1)
      XNORM=SQRT(XNORM*XNORM+T*T)

C  END OF SUBROUTINE LSI

   50 END