      SUBROUTINE LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML,
     * MODE, WS, IP)
C     DIMENSION W(MDW,N+1),PRGOPT(*),X(N),
C     WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
C     ABOVE, K=MAX(MA+MG,N).
C
C     WRITTEN BY R. J. HANSON AND K. H. HASKELL.  FOR FURTHER MATH.
C     AND ALGORITHMIC DETAILS SEE SANDIA LABORATORIES TECH. REPTS.
C     SAND77-0552, (1978), AND SAND78-1290, (1979).
C
C     ABSTRACT
C
C     THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES
C     PROBLEM WITH BOTH EQUALITY AND INEQUALITY CONSTRAINTS, AND, IF THE
C     USER REQUESTS, OBTAINS A COVARIANCE MATRIX OF THE SOLUTION
C     PARAMETERS.
C
C     SUPPOSE THERE ARE GIVEN MATRICES E, A AND G OF RESPECTIVE
C     DIMENSIONS ME BY N, MA BY N AND MG BY N, AND VECTORS F, B AND H OF
C     RESPECTIVE LENGTHS ME, MA AND MG.  THIS SUBROUTINE SOLVES THE
C     LINEARLY CONSTRAINED LEAST SQUARES PROBLEM
C
C                   EX = F, (E ME BY N) (EQUATIONS TO BE EXACTLY
C                                       SATISFIED)
C                   AX = B, (A MA BY N) (EQUATIONS TO BE
C                                       APPROXIMATELY SATISFIED,
C                                       LEAST SQUARES SENSE)
C                   GX.GE.H,(G MG BY N) (INEQUALITY CONSTRAINTS)
C
C     THE INEQUALITIES GX.GE.H MEAN THAT EVERY COMPONENT OF THE PRODUCT
C     GX MUST BE .GE. THE CORRESPONDING COMPONENT OF H.
C
C     IN CASE THE EQUALITY CONSTRAINTS CANNOT BE SATISFIED, A
C     GENERALIZED INVERSE SOLUTION RESIDUAL VECTOR LENGTH IS OBTAINED
C     FOR F-EX. THIS IS THE MINIMAL LENGTH POSSIBLE FOR F-EX.
C
C
C     ANY VALUES ME.GE.0, MA.GE.0, OR MG.GE.0 ARE PERMITTED.  THE
C     RANK OF THE MATRIX E IS ESTIMATED DURING THE COMPUTATION. WE CALL
C     THIS VALUE KRANKE. IT IS AN OUTPUT PARAMETER IN IP(1) DEFINED
C     BELOW. USING A GENERALIZED INVERSE SOLUTION OF EX=F, A REDUCED
C     LEAST SQUARES PROBLEM WITH INEQUALITY CONSTRAINTS IS OBTAINED.
C     THE TOLERANCES USED IN THESE TESTS FOR DETERMINING THE RANK
C     OF E AND THE RANK OF THE REDUCED LEAST SQUARES PROBLEM ARE
C     GIVEN IN SANDIA TECH. REPT. SAND78-1290. THEY CAN BE
C     MODIFIED BY THE USER IF NEW VALUES ARE PROVIDED IN
C     THE OPTION LIST OF THE ARRAY PRGOPT(*).
C
C     THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST..
C     W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
C     WHERE K=MAX(MA+MG,N).  THIS ALLOWS FOR A SOLUTION OF A RANGE OF
C     PROBLEMS IN THE GIVEN WORKING SPACE.  THE DIMENSION OF WS(*)
C     GIVEN IS A NECESSARY OVERESTIMATE.  ONCE A PARTICULAR PROBLEM
C     HAS BEEN RUN, THE OUTPUT PARAMETER IP(3) GIVES THE ACTUAL
C     DIMENSION REQUIRED FOR THAT PROBLEM.
C
C     THE PARAMETERS FOR LSEI( ) ARE
C
C     INPUT..
C
C     W(*,*),MDW,   THE ARRAY W(*,*) IS DOUBLY SUBSCRIPTED WITH
C     ME,MA,MG,N    FIRST DIMENSIONING PARAMETER EQUAL TO MDW.
C                   FOR THIS DISCUSSION LET US CALL M = ME+MA+MG.  THEN
C                   MDW MUST SATISFY MDW.GE.M.  THE CONDITION
C                   MDW.LT.M IS AN ERROR.
C
C                   THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS
C
C                                  (E  F)
C                                  (A  B)
C                                  (G  H)
C
C                   IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1
C                   RESPECTIVELY.
C
C                   THE INTEGERS ME, MA, AND MG ARE THE
C                   RESPECTIVE MATRIX ROW DIMENSIONS
C                   OF E, A AND G. EACH MATRIX HAS N COLUMNS.
C
C     PRGOPT(*)    THIS REAL-VALUED ARRAY IS THE OPTION VECTOR.
C                  IF THE USER IS SATISFIED WITH THE NOMINAL
C                  SUBPROGRAM FEATURES SET
C
C                  PRGOPT(1)=1 (OR PRGOPT(1)=1.0)
C
C                  OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF
C                  GROUPS OF DATA OF THE FOLLOWING FORM
C
C                  LINK
C                  KEY
C                  DATA SET
C
C                  THE PARAMETERS LINK AND KEY ARE EACH ONE WORD.
C                  THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS.
C                  THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY.
C                  THE VALUE OF LINK POINTS TO THE FIRST
C                  ENTRY OF THE NEXT GROUP OF DATA WITHIN
C                  PRGOPT(*).  THE EXCEPTION IS WHEN THERE ARE
C                  NO MORE OPTIONS TO CHANGE.  IN THAT
C                  CASE LINK=1 AND THE VALUES KEY AND DATA SET
C                  ARE NOT REFERENCED. THE GENERAL LAYOUT OF
C                  PRGOPT(*) IS AS FOLLOWS.
C
C               ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP)
C               .  PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE)
C               .       .
C               .       .
C               .       .
C               ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF
C               .                       NEXT GROUP)
C               .  PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(LINK1+2)=DATA VALUE
C               ...     .
C               .       .
C               .       .
C               ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS.
C                  A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR.
C                  THIS HELPS PREVENT USING INVALID BUT POSITIVE
C                  VALUES OF LINK THAT WILL PROBABLY EXTEND
C                  BEYOND THE PROGRAM LIMITS OF PRGOPT(*).
C                  UNRECOGNIZED VALUES OF KEY ARE IGNORED.  THE
C                  ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER
C                  OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING
C                  RESTRICTION.  TO PREVENT CYCLING IN THE
C                  PROCESSING OF THE OPTION ARRAY A COUNT OF THE
C                  NUMBER OF OPTIONS CHANGED IS MAINTAINED.
C                  WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR
C                  MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS.
C
C                  OPTIONS..
C
C                  KEY=1
C                         COMPUTE IN W(*,*) THE N BY N
C                  COVARIANCE MATRIX OF THE SOLUTION VARIABLES
C                  AS AN OUTPUT PARAMTER.  NOMINALLY THE
C                  COVARIANCE MATRIX WILL NOT BE COMPUTED.
C                  (THIS REQUIRES NO USER INPUT.)
C                  THE DATA SET FOR THIS OPTION IS A SINGLE VALUE.
C                  IT MUST BE NONZERO WHEN THE COVARIANCE MATRIX
C                  IS DESIRED.  IF IT IS ZERO, THE COVARIANCE
C                  MATRIX IS NOT COMPUTED.  WHEN THE COVARIANCE MATRIX
C                  IS COMPUTED, THE FIRST DIMENSIONING PARAMETER
C                  OF THE ARRAY W(*,*) MUST SATISFY MDW.GE.MAX0(M,N).
C
C                  KEY=2
C                         SCALE THE NONZERO COLUMNS OF THE
C                         ENTIRE DATA MATRIX.
C                  (E)
C                  (A)
C                  (G)
C
C                  TO HAVE LENGTH ONE.  THE DATA SET FOR THIS
C                  OPTION IS A SINGLE VALUE.  IT MUST BE
C                  NONZERO IF UNIT LENGTH COLUMN SCALING
C                  IS DESIRED.
C
C                  KEY=3
C                         SCALE COLUMNS OF THE ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  (G)
C
C                  WITH A USER-PROVIDED DIAGONAL MATRIX.
C                  THE DATA SET FOR THIS OPTION CONSISTS
C                  OF THE N DIAGONAL SCALING FACTORS, ONE FOR
C                  EACH MATRIX COLUMN.
C
C                  KEY=4
C                         CHANGE THE RANK DETERMINATION TOLERANCE FOR
C                  THE EQUALITY CONSTRAINT EQUATIONS FROM
C                  THE NOMINAL VALUE OF SQRT(EPS).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN EPS, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY EPS IS THE
C                  LARGEST POSITIVE NUMBER SUCH THAT T=1.+EPS
C                  SATISFIES T.EQ.1.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST EPS.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  KEY=5
C                         CHANGE THE RANK DETERMINATION TOLERANCE FOR
C                  THE REDUCED LEAST SQUARES EQUATIONS FROM
C                  THE NOMINAL VALUE OF SQRT(EPS).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN EPS, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST EPS.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  FOR EXAMPLE, SUPPOSE WE WANT TO CHANGE
C                  THE TOLERANCE FOR THE REDUCED LEAST SQUARES
C                  PROBLEM, COMPUTE THE COVARIANCE MATRIX OF
C                  THE SOLUTION PARAMETERS, AND PROVIDE
C                  COLUMN SCALING FOR THE DATA MATRIX.  FOR
C                  THESE OPTIONS THE DIMENSION OF PRGOPT(*)
C                  MUST BE AT LEAST N+9.  THE FORTRAN STATEMENTS
C                  DEFINING THESE OPTIONS WOULD BE AS FOLLOWS.
C
C                  PRGOPT(1)=4 (LINK TO ENTRY 4 IN PRGOPT(*))
C                  PRGOPT(2)=1 (COVARIANCE MATRIX KEY)
C                  PRGOPT(3)=1 (COVARIANCE MATRIX WANTED)
C
C                  PRGOPT(4)=7 (LINK TO ENTRY 7 IN PRGOPT(*))
C                  PRGOPT(5)=5 (LEAST SQUARES EQUAS. TOLERANCE KEY)
C                  PRGOPT(6)=... (NEW VALUE OF THE TOLERANCE)
C
C                  PRGOPT(7)=N+9 (LINK TO ENTRY N+9 IN PRGOPT(*))
C                  PRGOPT(8)=3 (USER-PROVIDED COLUMN SCALING KEY)
C
C                  CALL SCOPY(N,D,1,PRGOPT(9),1)  (COPY THE N
C                    SCALING FACTORS FROM THE USER ARRAY D(*)
C                    TO PRGOPT(9)-PRGOPT(N+8))
C
C                  PRGOPT(N+9)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  THE CONTENTS OF PRGOPT(*) ARE NOT MODIFIED
C                  BY THE SUBPROGRAM.
C                  THE OPTIONS FOR WNNLS( ) CAN ALSO BE INCLUDED
C                  IN THIS ARRAY.  THE VALUES OF KEY RECOGNIZED
C                  BY WNNLS( ) ARE 6, 7 AND 8.  THEIR FUNCTIONS
C                  ARE DOCUMENTED IN THE USAGE INSTRUCTIONS FOR
C                  SUBPROGRAM WNNLS( ).  NORMALLY THESE OPTIONS
C                  DO NOT NEED TO BE MODIFIED WHEN USING LSEI( ).
C
C     OUTPUT..
C
C     X(*),RNORME,  THE ARRAY X(*) CONTAINS THE SOLUTION PARAMETERS
C     RNORML        IF THE INTEGER OUTPUT FLAG MODE = 0 OR 1.
C                   THE DEFINITION OF MODE IS GIVEN DIRECTLY BELOW.
C                   WHEN MODE = 0 OR 1, RNORME AND RNORML
C                   RESPECTIVELY CONTAIN THE RESIDUAL VECTOR
C                   EUCLIDEAN LENGTHS OF F - EX AND B - AX.  WHEN
C                   MODE=1 THE EQUALITY CONSTRAINT EQUATIONS EX=F
C                   ARE CONTRADICTORY, SO RNORME.NE.0. THE RESIDUAL
C                   VECTOR F-EX HAS MINIMAL EUCLIDEAN LENGTH. FOR
C                   MODE.GE.2, NONE OF THESE PARAMETERS ARE
C                   DEFINED.
C
C     MODE          INTEGER FLAG THAT INDICATES THE SUBPROGRAM
C                   STATUS AFTER COMPLETION.  IF MODE.GE.2, NO
C                   SOLUTION HAS BEEN COMPUTED.
C
C                   MODE =
C
C                   0  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
C                      ARE COMPATIBLE AND HAVE BEEN SATISFIED.
C
C                   1  EQUALITY CONSTRAINTS ARE CONTRADICTORY.
C                      A GENERALIZED INVERSE SOLUTION OF EX=F WAS USED
C                      TO MINIMIZE THE RESIDUAL VECTOR LENGTH F-EX.
C                      IN THIS SENSE, THE SOLUTION IS STILL MEANINGFUL.
C
C                   2  INEQUALITY CONSTRAINTS ARE CONTRADICTORY.
C
C                   3  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
C                      ARE CONTRADICTORY
C
C                   4  USAGE ERROR OCCURRED.  THE VALUE
C                      OF MDW IS .LT. ME+MA+MG, MDW IS
C                      .LT. N AND A COVARIANCE MATRIX IS
C                      REQUESTED, OR THE OPTION VECTOR
C                      PRGOPT(*) IS NOT PROPERLY DEFINED.
C
C     W(*,*)        THE ARRAY W(*,*) CONTAINS THE N BY N SYMMETRIC
C                   COVARIANCE MATRIX OF THE SOLUTION PARAMETERS,
C                   PROVIDED THIS WAS REQUESTED ON INPUT WITH
C                   THE OPTION VECTOR PRGOPT(*) AND THE OUTPUT
C                   FLAG IS RETURNED WITH MODE = 0 OR 1.
C
C     IP(*)         THE INTEGER WORKING ARRAY HAS THREE ENTRIES
C                   THAT PROVIDE RANK AND WORKING ARRAY LENGTH
C                   INFORMATION AFTER COMPLETION.
C
C                      IP(1) = RANK OF EQUALITY CONSTRAINT
C                              MATRIX.  DEFINE THIS QUANTITY
C                              AS KRANKE.
C
C                      IP(2) = RANK OF REDUCED LEAST SQUARES
C                              PROBLEM.
C
C                      IP(3) = THE AMOUNT OF STORAGE IN THE
C                              WORKING ARRAY WS(*) THAT WAS
C                              ACTUALLY USED BY THE SUBPROGRAM.
C                              THE FORMULA GIVEN ABOVE FOR THE LENGTH
C                              OF WS(*) IS A NECESSARY OVERESTIMATE.
C                              IF EXACTLY THE SAME PROBLEM MATRICES
C                              ARE USED IN SUBSEQUENT EXECUTIONS,
C                              THE DECLARED DIMENSION OF WS(*) CAN
C                              BE REDUCED TO THIS OUTPUT VALUE.
C     USER DESIGNATED
C     WORKING ARRAYS..
C
C     WS(*),IP(*)              THESE ARE RESPECTIVELY TYPE REAL
C                              AND TYPE INTEGER WORKING ARRAYS.
C                              THEIR REQUIRED MINIMAL LENGTHS ARE
C                              GIVEN ABOVE.
C
C
C     SUBROUTINES CALLED
C
C     LSI           PART OF THIS PACKAGE.  SOLVES A
C                   CONSTRAINED LEAST SQUARES PROBLEM WITH
C                   INEQUALITY CONSTRAINTS.
C
C     SDOT,SSCAL,   SUBROUTINES FROM THE BLAS PACKAGE.
C     SAXPY,SASUM,  THESE ARE FOUND IN THE SANDIA MATH. LIBRARY.
C     SCOPY,SNRM2,  THESE CODES ARE ALSO DOCUMENTED IN THE SANDIA
C     SSWAP,ISAMAX  TECH. REPT., SAND77-0898.
C
C     H12           SUBROUTINE TO CONSTRUCT AND APPLY A
C                   HOUSEHOLDER TRANSFORMATION.
C
C     SPMPAR        FUNCTION TO COMPUTE THE RELATIVE MACHINE
C                   PRECISION.
C
C     REVISED FEBRUARY 26, 1980.
      DIMENSION W(MDW,1), PRGOPT(1), X(N), WS(1), IP(3)
      LOGICAL COV
      DATA ZERO/0.0/, HALF/0.5/, ONE/1.0/
C
C     COMPUTE MACHINE PRECISION EPS
      EPS = SPMPAR(1)
C
C     COMPUTE NUMBER OF POSSIBLE RIGHT MULTIPLYING HOUSEHOLDER
C     TRANSFORMATIONS.
   30 M = ME + MA + MG
      MODE = 0
      IF (N.LE.0 .OR. M.LE.0) RETURN
      IF (.NOT.(MDW.LT.M)) GO TO 40
      MODE = 4
      RETURN
   40 NP1 = N + 1
      KRANKE = MIN0(ME,N)
      N1 = 2*KRANKE + 1
      N2 = N1 + N
C
C     PROCESS-OPTION-VECTOR
      GO TO 410
   50 IF (.NOT.(COV .AND. MDW.LT.N)) GO TO 60
      MODE = 4
      RETURN
   60 L = KRANKE
C
C     COMPUTE NORM OF EQUALITY CONSTRAINT MATRIX AND RT SIDE.
      ENORM = ZERO
      DO 70 J=1,N
         ENORM = AMAX1(ENORM,SASUM(ME,W(1,J),1))
   70 CONTINUE
      FNORM = SASUM(ME,W(1,NP1),1)
      IF (.NOT.(L.GT.0)) GO TO 120
      SNMAX = ZERO
      RNMAX = ZERO
      DO 110 I=1,L
C
C     COMPUTE MAXIMUM RATIO OF VECTOR LENGTHS. PARTITION
C     IS AT COLUMN I.
      DO 80 K=I,ME
      SN = SDOT(N-I+1, W(K,I), MDW, W(K,I), MDW)
      RN = SDOT(I-1, W(K,1), MDW, W(K,1), MDW)
      IF (.NOT.(RN.EQ.ZERO.AND.SN.GT.SNMAX)) GO TO 85
      SNMAX = SN
      IMAX = K
      GO TO 80
   85 IF (.NOT.(K.EQ.I.OR.(SN*RNMAX.GT.RN*SNMAX))) GO TO 95
      SNMAX = SN
      RNMAX = RN
      IMAX = K
   95 CONTINUE
   80 CONTINUE
C
C     INTERCHANGE ROWS IF REQUIRED.
         IF (I.NE.IMAX) CALL SSWAP(NP1, W(I,1), MDW, W(IMAX,1), MDW)
         IF (.NOT.(SNMAX.GT.TAU**2*RNMAX)) GO TO 90
C
C     ELIMINATE ELEMS I+1,...,N IN ROW I.
         CALL H12(1, I, I+1, N, W(I,1), MDW, WS(I), W(I+1,1), MDW, 1,
     *    M-I)
         GO TO 100
   90    KRANKE = I - 1
         GO TO 130
  100    CONTINUE
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
C
C     SAVE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL SCOPY(KRANKE, W, MDW+1, WS(KRANKE+1), 1)
C
C     USE HOUSEHOLDER TRANS FROM LEFT TO ACHIEVE KRANKE BY KRANKE UPPER
C     TRIANGULAR FORM.
      IF (.NOT.(KRANKE.GT.0 .AND. KRANKE.LT.ME)) GO TO 150
      DO 140 KK=1,KRANKE
         K = KRANKE + 1 - KK
C
C     APPLY TRANFORMATION TO MATRIX COLS. 1,...,K-1.
         CALL H12(1, K, KRANKE+1, ME, W(1,K), 1, UP, W, 1, MDW, K-1)
C
C     APPLY TO RT SIDE VECTOR.
         CALL H12(2, K, KRANKE+1, ME, W(1,K), 1, UP, W(1,NP1), 1, 1, 1)
  140 CONTINUE
  150 IF (.NOT.(KRANKE.GT.0)) GO TO 170
C
C     SOLVE FOR VARIABLES 1,...,KRANKE IN NEW COORDINATES.
      CALL SCOPY(KRANKE, W(1,NP1), 1, X, 1)
      DO 160 I=1,KRANKE
         X(I) = (X(I)-SDOT(I-1,W(I,1),MDW,X,1))/W(I,I)
  160 CONTINUE
C
C     COMPUTE RESIDUALS FOR REDUCED PROBLEM.
  170 MEP1 = ME + 1
      RNORML = ZERO
      IF (.NOT.(ME.LT.M)) GO TO 190
      DO 180 I=MEP1,M
         W(I,NP1) = W(I,NP1) - SDOT(KRANKE,W(I,1),MDW,X,1)
  180 CONTINUE
C
C     COMPUTE EQUAL. CONSTRAINT EQUAS. RESIDUAL LENGTH.
  190 RNORME = SNRM2(ME-KRANKE,W(KRANKE+1,NP1),1)
C
C     MOVE REDUCED PROBLEM DATA UPWARD IF KRANKE.LT.ME.
      IF (.NOT.(KRANKE.LT.ME)) GO TO 210
      DO 200 J=1,NP1
         CALL SCOPY(M-KRANKE, W(ME+1,J), 1, W(KRANKE+1,J), 1)
  200 CONTINUE
C
C     COMPUTE SOLN OF REDUCED PROBLEM.
  210 CALL LSI(W(KRANKE+1,KRANKE+1), MDW, MA, MG, N-KRANKE, PRGOPT,
     * X(KRANKE+1), RNORML, MODE, WS(N2), IP(2))
      IF (.NOT.(ME.GT.0)) GO TO 250
C
C     TEST FOR CONSISTENCY OF EQUALITY CONSTRAINTS.
      MDEQC = 0
      XNRME = SASUM(KRANKE,W(1,NP1),1)
      IF (RNORME.GT.TAU*(ENORM*XNRME+FNORM)) MDEQC = 1
      MODE = MODE + MDEQC
C
C     CHECK IF SOLN TO EQUAL. CONSTRAINTS SATISFIES INEQUAL.
C     CONSTRAINTS WHEN THERE ARE NO DEGREES OF FREEDOM LEFT.
      IF (.NOT.(KRANKE.EQ.N .AND. MG.GT.0)) GO TO 240
      XNORM = SASUM(N,X,1)
      MAPME1 = MA + ME + 1
      DO 230 I=MAPME1,M
         T = SDOT(N,W(I,1),MDW,X,1) - W(I,NP1)
         SIZE = SASUM(N,W(I,1),MDW)*XNORM + ABS(W(I,NP1))
         IF (.NOT.(T.LT.(-TAU*SIZE))) GO TO 220
         MODE = MODE + 2
         GO TO 380
  220    CONTINUE
  230 CONTINUE
  240 CONTINUE
  250 IF (.NOT.(KRANKE.GT.0)) GO TO 340
C
C     REPLACE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL SCOPY(KRANKE, WS(KRANKE+1), 1, W, MDW+1)
C
C     REAPPLY TRANS TO PUT SOLN IN ORIGINAL COORDINATES.
      DO 260 II=1,KRANKE
         I = KRANKE + 1 - II
         CALL H12(2, I, I+1, N, W(I,1), MDW, WS(I), X, 1, 1, 1)
  260 CONTINUE
C
C     COMPUTE COV MATRIX OF EQUAL. CONSTRAINED PROBLEM.
      IF (.NOT.(COV)) GO TO 330
      DO 320 JJ=1,KRANKE
         J = KRANKE + 1 - JJ
         RB = WS(J)*W(J,J)
         IF (RB.NE.ZERO) RB = ONE/RB
         IF (.NOT.(J.LT.N)) GO TO 310
         JP1 = J + 1
         DO 270 I=JP1,N
            W(I,J) = SDOT(N-J,W(I,JP1),MDW,W(J,JP1),MDW)*RB
  270    CONTINUE
         GAM = SDOT(N-J,W(JP1,J),1,W(J,JP1),MDW)*RB
         GAM = HALF*GAM
         CALL SAXPY(N-J, GAM, W(J,JP1), MDW, W(JP1,J), 1)
         DO 290 I=JP1,N
            DO 280 K=I,N
               W(I,K) = W(I,K) + W(J,I)*W(K,J) + W(I,J)*W(J,K)
               W(K,I) = W(I,K)
  280       CONTINUE
  290    CONTINUE
         UJ = WS(J)
         VJ = GAM*UJ
         W(J,J) = UJ*VJ + UJ*VJ
         DO 300 I=JP1,N
            W(J,I) = UJ*W(I,J) + VJ*W(J,I)
  300    CONTINUE
         CALL SCOPY(N-J, W(J,JP1), MDW, W(JP1,J), 1)
  310    CONTINUE
  320 CONTINUE
  330 CONTINUE
C
C     APPLY THE SCALING TO THE COVARIANCE MATRIX.
  340 IF (.NOT.(COV)) GO TO 360
      DO 350 I=1,N
         L = N1 + I
         CALL SSCAL(N, WS(L-1), W(I,1), MDW)
         CALL SSCAL(N, WS(L-1), W(1,I), 1)
  350 CONTINUE
  360 CONTINUE
  380 CONTINUE
C
C     RESCALE SOLN. VECTOR.
      IF (.NOT.(MODE.LE.1)) GO TO 400
      DO 390 J=1,N
         L = N1 + J
         X(J) = X(J)*WS(L-1)
  390 CONTINUE
  400 IP(1) = KRANKE
      IP(3) = IP(3) + 2*KRANKE + N
      RETURN
  410 CONTINUE
C     TO PROCESS-OPTION-VECTOR
C
C     THE NOMINAL TOLERANCE USED IN THE CODE
C     FOR THE EQUALITY CONSTRAINT EQUATIONS.
      TAU = SQRT(EPS)
C
C     THE NOMINAL COLUMN SCALING USED IN THE CODE IS
C     THE IDENTITY SCALING.
      WS(N1) = ONE
      CALL SCOPY(N, WS(N1), 0, WS(N1), 1)
C
C     NO COVARIANCE MATRIX IS NOMINALLY COMPUTED.
      COV = .FALSE.
C
C     DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE.
      NOPT = 1000
      NTIMES = 0
C
C     DEFINE BOUND FOR POSITIVE VALUES OF LINK.
      NLINK = 100000
      LAST = 1
      LINK = PRGOPT(1)
      IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 420
      MODE = 4
      RETURN
  420 IF (.NOT.(LINK.GT.1)) GO TO 470
      NTIMES = NTIMES + 1
      IF (.NOT.(NTIMES.GT.NOPT)) GO TO 430
      MODE = 4
      RETURN
  430 KEY = PRGOPT(LAST+1)
      IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO
      IF (.NOT.(KEY.EQ.2 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 450
      DO 440 J=1,N
         T = SNRM2(M,W(1,J),1)
         IF (T.NE.ZERO) T = ONE/T
         L = N1 + J
         WS(L-1) = T
  440 CONTINUE
  450 IF (KEY.EQ.3) CALL SCOPY(N, PRGOPT(LAST+2), 1, WS(N1), 1)
      IF (KEY.EQ.4) TAU = AMAX1(EPS,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 460
      MODE = 4
      RETURN
  460 LAST = LINK
      LINK = NEXT
      GO TO 420
  470 DO 480 J=1,N
         L = N1 + J
         CALL SSCAL(M, WS(L-1), W(1,J), 1)
  480 CONTINUE
      GO TO 50
      END
