	SUBROUTINE LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML,
     * MODE, WS, IP)
	IMPLICIT INTEGER*4 (I-N)
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                  SUBROUTINE 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                   THE FOLLOWING INTERPRETATION OF
C                   MODE=1,2 OR 3 MUST BE MADE.  THE
C                   SETS CONSISTING OF ALL SOLUTIONS
C                   OF THE EQUALITY CONSTRAINTS EX=F
C                   AND ALL VECTORS SATISFYING GX.GE.H
C                   HAVE NO POINTS ON COMMON.  (IN
C                   PARTICULAR THIS DOES NOT SAY THAT
C                   EACH INDIVIDUAL SET HAS NO POINTS
C                   AT ALL, ALTHOUGH THIS COULD BE THE
C                   CASE.)
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.  LOCATED IN SANDIA
C                   MATH. LIBRARY.
C
C     XERROR        FROM SLATEC ERROR PROCESSING PACKAGE.
C                   THIS IS DOCUMENTED IN SANDIA TECH. REPT.,
C                   SAND78-1189.
C
C     SUBROUTINE LSEI(W,MDW,ME,MA,MG,N,PRGOPT,X,RNORME,RNORML,MODE,WS,
C    1 IP)
C
C     REVISED FEBRUARY 26, 1980
	DIMENSION W(MDW,1), PRGOPT(1), X(N), WS(1), IP(3)
	LOGICAL COV
	DATA ZERO /0.E0/, EPS /0.E0/, HALF /.5E0/, ONE /1.E0/
C
C     COMPUTE MACHINE PRECISION=EPS ONLY WHEN NECESSARY.
	IF (.NOT.(EPS.EQ.ZERO)) GO TO 30
	EPS = ONE
10	IF (ONE+EPS.EQ.ONE) GO TO 20
	EPS = EPS*HALF
	GO TO 10
20	EPS = EPS + EPS
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
	NERR = 1
	IOPT = 1
	CALL XERROR(36HLSEI( ), MDW.LT.ME+MA+MG IS AN ERROR, 36, NERR,
     * IOPT)
	MODE = 4
	RETURN
40	NP1 = N + 1
	KRANKE = MIN0(ME,N)
	N1 = 2*KRANKE + 1
	N2 = N1 + N
C
C     PROCESS-OPTION-VECTOR
	ASSIGN 50 TO I99994
	GO TO 400
50	IF (.NOT.(COV .AND. MDW.LT.N)) GO TO 60
	NERR = 2
	IOPT = 1
	CALL XERROR(
     * 54HLSEI( ), MDW.LT.N, WHEN COV MATRIX NEEDED, IS AN ERROR, 54,
     * NERR, IOPT)
	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 COL. 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 370
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
370	CONTINUE
C
C     RESCALE SOLN. VECTOR.
	IF (.NOT.(MODE.LE.1)) GO TO 390
	DO 380 J=1,N
	   L = N1 + J
	   X(J) = X(J)*WS(L-1)
380	CONTINUE
390	IP(1) = KRANKE
	IP(3) = IP(3) + 2*KRANKE + N
	RETURN
400	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 410
	NERR = 3
	IOPT = 1
	CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
	MODE = 4
	RETURN
410	IF (.NOT.(LINK.GT.1)) GO TO 460
	NTIMES = NTIMES + 1
	IF (.NOT.(NTIMES.GT.NOPT)) GO TO 420
	NERR = 3
	IOPT = 1
	CALL XERROR(
     * 52HLSEI( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 52,
     * NERR, IOPT)
	MODE = 4
	RETURN
420	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 440
	DO 430 J=1,N
	   T = SNRM2(M,W(1,J),1)
	   IF (T.NE.ZERO) T = ONE/T
	   L = N1 + J
	   WS(L-1) = T
430	CONTINUE
440	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 450
	NERR = 3
	IOPT = 1
	CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
	MODE = 4
	RETURN
450	LAST = LINK
	LINK = NEXT
	GO TO 410
460	DO 470 J=1,N
	   L = N1 + J
	   CALL SSCAL(M, WS(L-1), W(1,J), 1)
470	CONTINUE
	GO TO 480
480	GO TO I99994, (50)
	END

	REAL FUNCTION SNRM2 ( N, SX, INCX)
	IMPLICIT INTEGER*4 (I-N)
	INTEGER*4          NEXT
	REAL   SX(1),  CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
	DATA   ZERO, ONE /0.0E0, 1.0E0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
	DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
C
	IF(N .GT. 0) GO TO 10
	   SNRM2  = ZERO
	   GO TO 300
C
10	ASSIGN 30 TO NEXT
	SUM = ZERO
	NN = N * INCX
C                                                 BEGIN MAIN LOOP
	I = 1
20	   GO TO NEXT,(30, 50, 70, 110)
30	IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
	ASSIGN 50 TO NEXT
	XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
50	IF( SX(I) .EQ. ZERO) GO TO 200
	IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
	ASSIGN 70 TO NEXT
	GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
100	I = J
	ASSIGN 110 TO NEXT
	SUM = (SUM / SX(I)) / SX(I)
105	XMAX = ABS(SX(I))
	GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70	IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
110	IF( ABS(SX(I)) .LE. XMAX ) GO TO 115
	   SUM = ONE + SUM * (XMAX / SX(I))**2
	   XMAX = ABS(SX(I))
	   GO TO 200
C
115	SUM = SUM + (SX(I)/XMAX)**2
	GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
75	SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
85	HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
	DO 95 J =I,NN,INCX
	IF(ABS(SX(J)) .GE. HITEST) GO TO 100
95	   SUM = SUM + SX(J)**2
	SNRM2 = SQRT( SUM )
	GO TO 300
C
200	CONTINUE
	I = I + INCX
	IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
	SNRM2 = XMAX * SQRT(SUM)
300	CONTINUE
	RETURN
	END

	SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)
	IMPLICIT INTEGER*4 (I-N)
C
C     COPY SINGLE PRECISION SX TO SINGLE PRECISION SY.
C     FOR I = 0 TO N-1, COPY  SX(LX+I*INCX) TO SY(LY+I*INCY),
C     WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C     DEFINED IN A SIMILAR WAY USING INCY.
C
	REAL SX(1),SY(1)
	IF(N.LE.0)RETURN
	IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
5	CONTINUE
C
C        CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
	IX = 1
	IY = 1
	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
	DO 10 I = 1,N
	  SY(IY) = SX(IX)
	  IX = IX + INCX
	  IY = IY + INCY
10	CONTINUE
	RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7.
C
20	M = MOD(N,7)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  SY(I) = SX(I)
30	CONTINUE
	IF( N .LT. 7 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,7
	  SY(I) = SX(I)
	  SY(I + 1) = SX(I + 1)
	  SY(I + 2) = SX(I + 2)
	  SY(I + 3) = SX(I + 3)
	  SY(I + 4) = SX(I + 4)
	  SY(I + 5) = SX(I + 5)
	  SY(I + 6) = SX(I + 6)
50	CONTINUE
	RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
60	CONTINUE
	NS = N*INCX
	    DO 70 I=1,NS,INCX
	    SY(I) = SX(I)
70	    CONTINUE
	RETURN
	END

	SUBROUTINE SSCAL(N,SA,SX,INCX)
	IMPLICIT INTEGER*4 (I-N)
C
C     REPLACE SINGLE PRECISION SX BY SINGLE PRECISION SA*SX.
C     FOR I = 0 TO N-1, REPLACE SX(1+I*INCX) WITH  SA * SX(1+I*INCX)
C
	REAL SA,SX(1)
	IF(N.LE.0)RETURN
	IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
	NS = N*INCX
	    DO 10 I = 1,NS,INCX
	    SX(I) = SA*SX(I)
10	    CONTINUE
	RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
20	M = MOD(N,5)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  SX(I) = SA*SX(I)
30	CONTINUE
	IF( N .LT. 5 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,5
	  SX(I) = SA*SX(I)
	  SX(I + 1) = SA*SX(I + 1)
	  SX(I + 2) = SA*SX(I + 2)
	  SX(I + 3) = SA*SX(I + 3)
	  SX(I + 4) = SA*SX(I + 4)
50	CONTINUE
	RETURN
	END

	SUBROUTINE LSI(W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, IP)
	IMPLICIT INTEGER*4 (I-N)
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     REVISED SEP. 6, 1979.
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.  PART OF THIS PACKAGE.
C
C     H12           SUBROUTINE TO CONSTRUCT AND APPLY A HOUSEHOLDER
C                   TRANSFORMATION.  LOCATED IN SANDIA MATH. LIB.
C
C     SUBROUTINE LSI(W,MDW,MA,MG,N,PRGOPT,X,RNORM,MODE,WS,IP)
C
	DIMENSION W(MDW,1), PRGOPT(1), X(N), IP(N), WS(1)
	LOGICAL COV
C
	DATA ZERO /0.E0/, EPS /0.E0/, ONE /1.E0/, HALF /.5E0/
C
C     COMPUTE MACHINE PRECISION=EPS ONLY WHEN NECESSARY.
	IF (.NOT.(EPS.EQ.ZERO)) GO TO 30
	EPS = ONE
10	IF (ONE+EPS.EQ.ONE) GO TO 20
	EPS = EPS*HALF
	GO TO 10
20	EPS = EPS + EPS
30	MODE = 0
	M = MA + MG
	NP1 = N + 1
	KRANK = 0
	IF (N.LE.0 .OR. M.LE.0) GO TO 70
	ASSIGN 40 TO I99994
	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, RNORM, WS(N2),
     * WS(N1), IP)
	FAC = ONE
	IF (KRANK.LT.MA) FAC = RNORM**2/FLOAT(MA-KRANK)
	ASSIGN 60 TO I99990
	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 520
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 500
	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
500	GO TO 510
510	GO TO I99994, (40)
520	GO TO I99990, (60)
	END

	SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP)
	IMPLICIT INTEGER*4 (I-N)
C     DIMENSION A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N)
C
C     WRITTEN BY C. L. LAWSON AND R. J. HANSON.  FROM THE BOOK SOLVING
C     LEAST SQUARES PROBLEMS, PRENTICE-HALL, INC. (1974). FOR FURTHER
C     ALGORITHMIC DETAILS SEE ALGORITHM HFTI IN CHAPTER 14.
C
C     ABSTRACT
C
C     THIS SUBROUTINE SOLVES A LINEAR LEAST SQUARES PROBLEM OR A SET OF
C     LINEAR LEAST SQUARES PROBLEMS HAVING THE SAME MATRIX BUT DIFFERENT
C     RIGHT-SIDE VECTORS.  THE PROBLEM DATA CONSISTS OF AN M BY N MATRIX
C     A, AN M BY NB MATRIX B, AND AN ABSOLUTE TOLERANCE PARAMETER TAU
C     WHOSE USAGE IS DESCRIBED BELOW.  THE NB COLUMN VECTORS OF B
C     REPRESENT RIGHT-SIDE VECTORS FOR NB DISTINCT LINEAR LEAST SQUARES
C     PROBLEMS.
C
C     THIS SET OF PROBLEMS CAN ALSO BE WRITTEN AS THE MATRIX LEAST
C     SQUARES PROBLEM
C
C                       AX = B,
C
C     WHERE X IS THE N BY NB SOLUTION MATRIX.
C
C     NOTE THAT IF B IS THE M BY M IDENTITY MATRIX, THEN X WILL BE THE
C     PSEUDO-INVERSE OF A.
C
C     THIS SUBROUTINE FIRST TRANSFORMS THE AUGMENTED MATRIX (A B) TO A
C     MATRIX (R C) USING PREMULTIPLYING HOUSEHOLDER TRANSFORMATIONS WITH
C     COLUMN INTERCHANGES.  ALL SUBDIAGONAL ELEMENTS IN THE MATRIX R ARE
C     ZERO AND ITS DIAGONAL ELEMENTS SATISFY
C
C                       ABS(R(I,I)).GE.ABS(R(I+1,I+1)),
C
C                       I = 1,...,L-1, WHERE
C
C                       L = MIN(M,N).
C
C     THE SUBROUTINE WILL COMPUTE AN INTEGER, KRANK, EQUAL TO THE NUMBER
C     OF DIAGONAL TERMS OF R THAT EXCEED TAU IN MAGNITUDE.  THEN A
C     SOLUTION OF MINIMUM EUCLIDEAN LENGTH IS COMPUTED USING THE FIRST
C     KRANK ROWS OF (R C).
C
C     TO BE SPECIFIC WE SUGGEST THAT THE USER CONSIDER AN EASILY
C     COMPUTABLE MATRIX NORM, SUCH AS, THE MAXIMUM OF ALL COLUMN SUMS OF
C     MAGNITUDES.
C
C     NOW IF THE RELATIVE UNCERTAINTY OF B IS EPS, (NORM OF UNCERTAINTY/
C     NORM OF B), IT IS SUGGESTED THAT TAU BE SET APPROXIMATELY EQUAL TO
C     EPS*(NORM OF A).
C
C     THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST..
C     A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N).  THIS
C     PERMITS THE SOLUTION OF A RANGE OF PROBLEMS IN THE SAME ARRAY
C     SPACE.
C
C     THE ENTIRE SET OF PARAMETERS FOR HFTI ARE
C
C     INPUT..
C
C     A(*,*),MDA,M,N    THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY N
C                       MATRIX A OF THE LEAST SQUARES PROBLEM AX = B.
C                       THE FIRST DIMENSIONING PARAMETER OF THE ARRAY
C                       A(*,*) IS MDA, WHICH MUST SATISFY MDA.GE.M
C                       EITHER M.GE.N OR M.LT.N IS PERMITTED.  THERE
C                       IS NO RESTRICTION ON THE RANK OF A.  THE
C                       CONDITION MDA.LT.M IS CONSIDERED AN ERROR.
C
C     B(*),MDB,NB       IF NB = 0 THE SUBROUTINE WILL PERFORM THE
C                       ORTHOGONAL DECOMPOSITION BUT WILL MAKE NO
C                       REFERENCES TO THE ARRAY B(*).  IF NB.GT.0
C                       THE ARRAY B(*) MUST INITIALLY CONTAIN THE M BY
C                       NB MATRIX B OF THE LEAST SQUARES PROBLEM AX =
C                       B.  IF NB.GE.2 THE ARRAY B(*) MUST BE DOUBLY
C                       SUBSCRIPTED WITH FIRST DIMENSIONING PARAMETER
C                       MDB.GE.MAX(M,N).  IF NB = 1 THE ARRAY B(*) MAY
C                       BE EITHER DOUBLY OR SINGLY SUBSCRIPTED.  IN
C                       THE LATTER CASE THE VALUE OF MDB IS ARBITRARY
C                       BUT IT SHOULD BE SET TO SOME VALID INTEGER
C                       VALUE SUCH AS MDB = M.
C
C                       THE CONDITION OF NB.GT.1.AND.MDB.LT. MAX(M,N)
C                       IS CONSIDERED AN ERROR.
C
C     TAU               ABSOLUTE TOLERANCE PARAMETER PROVIDED BY USER
C                       FOR PSEUDORANK DETERMINATION.
C
C     H(*),G(*),IP(*)   ARRAYS OF WORKING SPACE USED BY HFTI.
C
C     OUTPUT..
C
C     A(*,*)            THE CONTENTS OF THE ARRAY A(*,*) WILL BE
C                       MODIFIED BY THE SUBROUTINE.  THESE CONTENTS
C                       ARE NOT GENERALLY REQUIRED BY THE USER.
C
C     B(*)              ON RETURN THE ARRAY B(*) WILL CONTAIN THE N BY
C                       NB SOLUTION MATRIX X.
C
C     KRANK             SET BY THE SUBROUTINE TO INDICATE THE
C                       PSEUDORANK OF A.
C
C     RNORM(*)          ON RETURN, RNORM(J) WILL CONTAIN THE EUCLIDEAN
C                       NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM
C                       DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY
C                       B(*,*) FOR J = 1,...,NB.
C
C     H(*),G(*)         ON RETURN THESE ARRAYS RESPECTIVELY CONTAIN
C                       ELEMENTS OF THE PRE- AND POST-MULTIPLYING
C                       HOUSEHOLDER TRANSFORMATIONS USED TO COMPUTE
C                       THE MINIMUM EUCLIDEAN LENGTH SOLUTION.
C
C     IP(*)             ARRAY IN WHICH THE SUBROUTINE RECORDS INDICES
C                       DESCRIBING THE PERMUTATION OF COLUMN VECTORS.
C                       THE CONTENTS OF ARRAYS H(*),G(*) AND IP(*)
C                       ARE NOT GENERALLY REQUIRED BY THE USER.
C
	DIMENSION A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB)
	INTEGER IP(N)
	DOUBLE PRECISION SM,DZERO
	SZERO=0.
	DZERO=0.D0
	FACTOR=0.001
C
	K=0
	LDIAG=MIN0(M,N)
	IF (LDIAG.LE.0) GO TO 270
	IF (.NOT.MDA.LT.M) GO TO 5
	NERR=1
	IOPT=2
	CALL XERROR(31HHFTI MDA.LT.M.. PROBABLE ERROR.,31,NERR,IOPT)
	RETURN
5	CONTINUE
C
	IF (.NOT.(NB.GT.1.AND.MAX0(M,N).GT.MDB)) GO TO 6
	NERR=2
	IOPT=2
	CALL XERROR(49HHFTI MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR.,
     1 49,NERR,IOPT)
	RETURN
6	CONTINUE
C
	    DO 80 J=1,LDIAG
	    IF (J.EQ.1) GO TO 20
C
C     UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
C    ..
	    LMAX=J
	        DO 10 L=J,N
	        H(L)=H(L)-A(J-1,L)**2
	        IF (H(L).GT.H(LMAX)) LMAX=L
10	        CONTINUE
	    IF(HMAX+FACTOR*H(LMAX) .GT. HMAX) GO TO 50
C
C     COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
C    ..
20	    LMAX=J
	        DO 40 L=J,N
	        H(L)=0.
	            DO 30 I=J,M
30	            H(L)=H(L)+A(I,L)**2
	        IF (H(L).GT.H(LMAX)) LMAX=L
40	        CONTINUE
	    HMAX=H(LMAX)
C    ..
C     LMAX HAS BEEN DETERMINED
C
C     DO COLUMN INTERCHANGES IF NEEDED.
C    ..
50	    CONTINUE
	    IP(J)=LMAX
	    IF (IP(J).EQ.J) GO TO 70
	        DO 60 I=1,M
	        TMP=A(I,J)
	        A(I,J)=A(I,LMAX)
60	        A(I,LMAX)=TMP
	    H(LMAX)=H(J)
C
C     COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.
C    ..
70	    CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)
80	    CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)
C
C     DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
C    ..
	    DO 90 J=1,LDIAG
	    IF (ABS(A(J,J)).LE.TAU) GO TO 100
90	    CONTINUE
	K=LDIAG
	GO TO 110
100	K=J-1
110	KP1=K+1
C
C     COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
C
	IF (NB.LE.0) GO TO 140
	    DO 130 JB=1,NB
	    TMP=SZERO
	    IF (KP1.GT.M) GO TO 130
	        DO 120 I=KP1,M
120	        TMP=TMP+B(I,JB)**2
130	    RNORM(JB)=SQRT(TMP)
140	CONTINUE
C                                           SPECIAL FOR PSEUDORANK = 0
	IF (K.GT.0) GO TO 160
	IF (NB.LE.0) GO TO 270
	    DO 150 JB=1,NB
	        DO 150 I=1,N
150	        B(I,JB)=SZERO
	GO TO 270
C
C     IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
C     DECOMPOSITION OF FIRST K ROWS.
C    ..
160	IF (K.EQ.N) GO TO 180
	    DO 170 II=1,K
	    I=KP1-II
170	    CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)
180	CONTINUE
C
C
	IF (NB.LE.0) GO TO 270
	    DO 260 JB=1,NB
C
C     SOLVE THE K BY K TRIANGULAR SYSTEM.
C    ..
	        DO 210 L=1,K
	        SM=DZERO
	        I=KP1-L
	        IF (I.EQ.K) GO TO 200
	        IP1=I+1
	            DO 190 J=IP1,K
190	            SM=SM+A(I,J)*DBLE(B(J,JB))
200	        SM1=SM
210	        B(I,JB)=(B(I,JB)-SM1)/A(I,I)
C
C     COMPLETE COMPUTATION OF SOLUTION VECTOR.
C    ..
	    IF (K.EQ.N) GO TO 240
	        DO 220 J=KP1,N
220	        B(J,JB)=SZERO
	        DO 230 I=1,K
230	        CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1)
C
C      RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
C      COLUMN INTERCHANGES.
C    ..
240	        DO 250 JJ=1,LDIAG
	        J=LDIAG+1-JJ
	        IF (IP(J).EQ.J) GO TO 250
	        L=IP(J)
	        TMP=B(L,JB)
	        B(L,JB)=B(J,JB)
	        B(J,JB)=TMP
250	        CONTINUE
260	    CONTINUE
C    ..
C     THE SOLUTION VECTORS, X, ARE NOW
C     IN THE FIRST  N  ROWS OF THE ARRAY B(,).
C
270	KRANK=K
	RETURN
	END

	SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
	IMPLICIT INTEGER*4 (I-N)
	DIMENSION U(IUE,M), C(1)
	DOUBLE PRECISION SM,B
	ONE=1.
C
	IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
	CL=ABS(U(1,LPIVOT))
	IF (MODE.EQ.2) GO TO 60
C                            ****** CONSTRUCT THE TRANSFORMATION. ******
	    DO 10 J=L1,M
10	    CL=AMAX1(ABS(U(1,J)),CL)
	IF (CL) 130,130,20
20	CLINV=ONE/CL
	SM=(DBLE(U(1,LPIVOT))*CLINV)**2
	    DO 30 J=L1,M
30	    SM=SM+(DBLE(U(1,J))*CLINV)**2
C                              CONVERT DBLE. PREC. SM TO SNGL. PREC. SM1
	SM1=SM
	CL=CL*SQRT(SM1)
	IF (U(1,LPIVOT)) 50,50,40
40	CL=-CL
50	UP=U(1,LPIVOT)-CL
	U(1,LPIVOT)=CL
	GO TO 70
C            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******
C
60	IF (CL) 130,130,70
70	IF (NCV.LE.0) RETURN
	B=DBLE(UP)*U(1,LPIVOT)
C                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.
C
	IF (B) 80,130,130
80	B=ONE/B
	I2=1-ICV+ICE*(LPIVOT-1)
	INCR=ICE*(L1-LPIVOT)
	    DO 120 J=1,NCV
	    I2=I2+ICV
	    I3=I2+INCR
	    I4=I3
	    SM=C(I2)*DBLE(UP)
	        DO 90 I=L1,M
	        SM=SM+C(I3)*DBLE(U(1,I))
90	        I3=I3+ICE
	    IF (SM) 100,120,100
100	    SM=SM*B
	    C(I2)=C(I2)+SM*DBLE(UP)
	        DO 110 I=L1,M
	        C(I4)=C(I4)+SM*DBLE(U(1,I))
110	        I4=I4+ICE
120	    CONTINUE
130	RETURN
	END

	SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
	IMPLICIT INTEGER*4 (I-N)
C
C     OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY.
C     FOR I = 0 TO N-1, REPLACE  SY(LY+I*INCY) WITH SA*SX(LX+I*INCX) +
C       SY(LY+I*INCY), WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N,
C       AND LY IS DEFINED IN A SIMILAR WAY USING INCY.
C
	REAL SX(1),SY(1),SA
	IF(N.LE.0.OR.SA.EQ.0.E0) RETURN
	IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
5	CONTINUE
C
C        CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS.
C
	IX = 1
	IY = 1
	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
	DO 10 I = 1,N
	  SY(IY) = SY(IY) + SA*SX(IX)
	  IX = IX + INCX
	  IY = IY + INCY
10	CONTINUE
	RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4.
C
20	M = MOD(N,4)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  SY(I) = SY(I) + SA*SX(I)
30	CONTINUE
	IF( N .LT. 4 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,4
	  SY(I) = SY(I) + SA*SX(I)
	  SY(I + 1) = SY(I + 1) + SA*SX(I + 1)
	  SY(I + 2) = SY(I + 2) + SA*SX(I + 2)
	  SY(I + 3) = SY(I + 3) + SA*SX(I + 3)
50	CONTINUE
	RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
60	CONTINUE
	NS = N*INCX
	    DO 70 I=1,NS,INCX
	    SY(I) = SA*SX(I) + SY(I)
70	    CONTINUE
	RETURN
	END

	SUBROUTINE SSWAP (N,SX,INCX,SY,INCY)
	IMPLICIT INTEGER*4 (I-N)
C
C     INTERCHANGE SINGLE PRECISION SX AND SINGLE PRECISION SY.
C     FOR I = 0 TO N-1, INTERCHANGE  SX(LX+I*INCX) AND SY(LY+I*INCY),
C     WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C     DEFINED IN A SIMILAR WAY USING INCY.
C
	REAL SX(1),SY(1),STEMP1,STEMP2,STEMP3
	IF(N.LE.0)RETURN
	IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
5	CONTINUE
C
C       CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
	IX = 1
	IY = 1
	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
	DO 10 I = 1,N
	  STEMP1 = SX(IX)
	  SX(IX) = SY(IY)
	  SY(IY) = STEMP1
	  IX = IX + INCX
	  IY = IY + INCY
10	CONTINUE
	RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C       CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3.
C
20	M = MOD(N,3)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  STEMP1 = SX(I)
	  SX(I) = SY(I)
	  SY(I) = STEMP1
30	CONTINUE
	IF( N .LT. 3 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,3
	  STEMP1 = SX(I)
	  STEMP2 = SX(I+1)
	  STEMP3 = SX(I+2)
	  SX(I) = SY(I)
	  SX(I+1) = SY(I+1)
	  SX(I+2) = SY(I+2)
	  SY(I) = STEMP1
	  SY(I+1) = STEMP2
	  SY(I+2) = STEMP3
50	CONTINUE
	RETURN
60	CONTINUE
C
C     CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
	NS = N*INCX
	  DO 70 I=1,NS,INCX
	  STEMP1 = SX(I)
	  SX(I) = SY(I)
	  SY(I) = STEMP1
70	  CONTINUE
	RETURN
	END

	SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        XERROR PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER
C        DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE
C        OF THE LIBRARY ERROR CONTROL FLAG, KONTRL.
C        (SEE SUBROUTINE XSETF FOR DETAILS.)
C
C     DESCRIPTION OF PARAMETERS
C      --INPUT--
C        MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED, CONTAINING
C                NO MORE THAN 72 CHARACTERS.
C        NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG.
C        NERR  - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE.
C                NERR MUST NOT BE ZERO.
C        LEVEL - ERROR CATEGORY.
C                =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR.
C                =1 MEANS THIS IS A RECOVERABLE ERROR.  (I.E., IT IS
C                   NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.)
C                =0 MEANS THIS IS A WARNING MESSAGE ONLY.
C                =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE
C                   PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY
C                   TIMES THIS CALL IS EXECUTED.
C
C     EXAMPLES
C        CALL XERROR(23HSMOOTH -- NUM WAS ZERO.,23,1,2)
C        CALL XERROR(43HINTEG  -- LESS THAN FULL ACCURACY ACHIEVED.,
C                    43,2,1)
C        CALL XERROR(65HROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL
C    1 FULLY COLLAPSED.,65,3,0)
C        CALL XERROR(39HEXP    -- UNDERFLOWS BEING SET TO ZERO.,39,1,-1)
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     LATEST REVISION ---  7 FEB 1979
C
	DIMENSION MESSG(NMESSG)
	CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.)
	RETURN
	END

	REAL FUNCTION SASUM(N,SX,INCX)
	IMPLICIT INTEGER*4 (I-N)
C
C     RETURNS SUM OF MAGNITUDES OF SINGLE PRECISION SX.
C     SASUM = SUM FROM 0 TO N-1 OF  ABS(SX(1+I*INCX))
C
	REAL SX(1)
	SASUM = 0.0E0
	IF(N.LE.0)RETURN
	IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
	NS = N*INCX
	    DO 10 I=1,NS,INCX
	    SASUM = SASUM + ABS(SX(I))
10	    CONTINUE
	RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6.
C
20	M = MOD(N,6)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  SASUM = SASUM + ABS(SX(I))
30	CONTINUE
	IF( N .LT. 6 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,6
	  SASUM = SASUM + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2))
     1  + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5))
50	CONTINUE
	RETURN
	END

	REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
	IMPLICIT INTEGER*4 (I-N)
C
C     RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY.
C     SDOT = SUM FOR I = 0 TO N-1 OF  SX(LX+I*INCX) * SY(LY+I*INCY),
C     WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C     DEFINED IN A SIMILAR WAY USING INCY.
C
	REAL SX(1),SY(1)
	SDOT = 0.0E0
	IF(N.LE.0)RETURN
	IF(INCX.EQ.INCY) IF(INCX-1)5,20,60
5	CONTINUE
C
C        CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS.
C
	IX = 1
	IY = 1
	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
	DO 10 I = 1,N
	  SDOT = SDOT + SX(IX)*SY(IY)
	  IX = IX + INCX
	  IY = IY + INCY
10	CONTINUE
	RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
20	M = MOD(N,5)
	IF( M .EQ. 0 ) GO TO 40
	DO 30 I = 1,M
	  SDOT = SDOT + SX(I)*SY(I)
30	CONTINUE
	IF( N .LT. 5 ) RETURN
40	MP1 = M + 1
	DO 50 I = MP1,N,5
	  SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     1   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
50	CONTINUE
	RETURN
C
C        CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
60	CONTINUE
	NS=N*INCX
	DO 70 I=1,NS,INCX
	  SDOT = SDOT + SX(I)*SY(I)
70	  CONTINUE
	RETURN
	END

	SUBROUTINE LPDP(A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, IS)
	IMPLICIT INTEGER*4 (I-N)
C     DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1),
C     WHERE N=N1+N2.  THIS IS A SLIGHT OVERESTIMATE FOR WS(*).
C
C     WRITTEN BY R. J. HANSON AND K. H. HASKELL, SANDIA LABS
C     REVISED SEP. 6, 1979.
C
C     DETERMINE AN N1-VECTOR W, AND
C               AN N2-VECTOR Z
C     WHICH MINIMIZES THE EUCLIDEAN LENGTH OF W
C     SUBJECT TO G*W+H*Z .GE. Y.
C     THIS IS THE LEAST PROJECTED DISTANCE PROBLEM, LPDP.
C     THE MATRICES G AND H ARE OF RESPECTIVE
C     DIMENSIONS M BY N1 AND M BY N2.
C
C     CALLED BY SUBPROGRAM LSI( ).
C
C     THE MATRIX
C                (G H Y)
C
C     OCCUPIES ROWS 1,...,M AND COLS 1,...,N1+N2+1 OF A(*,*).
C
C     THE SOLUTION (W) IS RETURNED IN X(*).
C                  (Z)
C
C     THE VALUE OF MODE INDICATES THE STATUS OF
C     THE COMPUTATION AFTER RETURNING TO THE USER.
C
C          MODE=1  THE SOLUTION WAS SUCCESSFULLY OBTAINED.
C
C          MODE=2  THE INEQUALITIES ARE INCONSISTENT.
C
C     SUBROUTINES CALLED
C
C     WNNLS         SOLVES A NONNEGATIVELY CONSTRAINED LINEAR LEAST
C                   SQUARES PROBLEM WITH LINEAR EQUALITY CONSTRAINTS.
C                   PART OF THIS PACKAGE.
C
C     SDOT,         SUBROUTINES FROM THE BLAS PACKAGE.
C     SSCAL,SNRM2,  THESE ARE FOUND IN THE SANDIA MATH. LIBRARY,
C     SCOPY         AND ARE ALSO DOCUMENTED IN THE SANDIA TECH.
C                   REPT., SAND77-0898.
C
	DIMENSION A(MDA,1), PRGOPT(1), X(1), WS(1), IS(1)
	DATA ZERO, ONE /0.E0,1.E0/, FAC /0.1E0/
	N = N1 + N2
	MODE = 1
	IF (.NOT.(M.LE.0)) GO TO 20
	IF (.NOT.(N.GT.0)) GO TO 10
	X(1) = ZERO
	CALL SCOPY(N, X, 0, X, 1)
10	WNORM = ZERO
	RETURN
20	NP1 = N + 1
C
C     SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE.
	DO 40 I=1,M
	   SC = SNRM2(N,A(I,1),MDA)
	   IF (.NOT.(SC.NE.ZERO)) GO TO 30
	   SC = ONE/SC
	   CALL SSCAL(NP1, SC, A(I,1), MDA)
30	   CONTINUE
40	CONTINUE
C
C     SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO).
	YNORM = SNRM2(M,A(1,NP1),1)
	IF (.NOT.(YNORM.NE.ZERO)) GO TO 50
	SC = ONE/YNORM
	CALL SSCAL(M, SC, A(1,NP1), 1)
C
C     SCALE COLS OF MATRIX H.
50	J = N1 + 1
60	IF (.NOT.(J.LE.N)) GO TO 70
	SC = SNRM2(M,A(1,J),1)
	IF (SC.NE.ZERO) SC = ONE/SC
	CALL SSCAL(M, SC, A(1,J), 1)
	X(J) = SC
	J = J + 1
	GO TO 60
70	IF (.NOT.(N1.GT.0)) GO TO 130
C
C     COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*).
	IW = 0
	DO 80 I=1,M
C
C     MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY.
	   CALL SCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1)
	   IW = IW + N2
C
C     MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY.
	   CALL SCOPY(N1, A(I,1), MDA, WS(IW+1), 1)
	   IW = IW + N1
C
C     MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY.
	   WS(IW+1) = A(I,NP1)
	   IW = IW + 1
80	CONTINUE
	WS(IW+1) = ZERO
	CALL SCOPY(N, WS(IW+1), 0, WS(IW+1), 1)
	IW = IW + N
	WS(IW+1) = ONE
	IW = IW + 1
C
C     SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0.  THE
C     MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR
C     F = TRANSPOSE OF (0,...,0,1).
	IX = IW + 1
	IW = IW + M
	CALL WNNLS(WS, NP1, N2, NP1-N2, M, 0, PRGOPT, WS(IX), RNORM,
     * MODEW, IS, WS(IW+1))
C
C     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W.
	SC = ONE - SDOT(M,A(1,NP1),1,WS(IX),1)
	IF (.NOT.(ONE+FAC*ABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 110
	SC = ONE/SC
	DO 90 J=1,N1
	   X(J) = SC*SDOT(M,A(1,J),1,WS(IX),1)
90	CONTINUE
C
C     COMPUTE THE VECTOR Q=Y-GW.  OVERWRITE Y WITH THIS VECTOR.
	DO 100 I=1,M
	   A(I,NP1) = A(I,NP1) - SDOT(N1,A(I,1),MDA,X,1)
100	CONTINUE
	GO TO 120
110	MODE = 2
	RETURN
120	CONTINUE
130	IF (.NOT.(N2.GT.0)) GO TO 180
C
C     COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*).
	IW = 0
	DO 140 I=1,M
	   CALL SCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1)
	   IW = IW + N2
	   WS(IW+1) = A(I,NP1)
	   IW = IW + 1
140	CONTINUE
	WS(IW+1) = ZERO
	CALL SCOPY(N2, WS(IW+1), 0, WS(IW+1), 1)
	IW = IW + N2
	WS(IW+1) = ONE
	IW = IW + 1
	IX = IW + 1
	IW = IW + M
C
C     SOLVE RV=S SUBJECT TO V.GE.0.  THE MATRIX R =(TRANSPOSE
C     OF (H Q)), WHERE Q=Y-GW.  THE (N2+1)-VECTOR S =(TRANSPOSE
C     OF (0,...,0,1)).
	CALL WNNLS(WS, N2+1, 0, N2+1, M, 0, PRGOPT, WS(IX), RNORM, MODEW,
     * IS, WS(IW+1))
C
C     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z.
	SC = ONE - SDOT(M,A(1,NP1),1,WS(IX),1)
	IF (.NOT.(ONE+FAC*ABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 160
	SC = ONE/SC
	DO 150 J=1,N2
	   L = N1 + J
	   X(L) = SC*SDOT(M,A(1,L),1,WS(IX),1)*X(L)
150	CONTINUE
	GO TO 170
160	MODE = 2
	RETURN
170	CONTINUE
C
C     ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION.
180	CALL SSCAL(N, YNORM, X, 1)
	WNORM = SNRM2(N1,X,1)
	RETURN
	END

	SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        XERRWV PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER
C        DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE
C        OF THE LIBRARY ERROR CONTROL FLAG, KONTRL.
C        (SEE SUBROUTINE XSETF FOR DETAILS.)
C        IN ADDITION, UP TO TWO INTEGER VALUES AND TWO REAL
C        VALUES MAY BE PRINTED ALONG WITH THE MESSAGE.
C
C     DESCRIPTION OF PARAMETERS
C      --INPUT--
C        MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED.
C        NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG.
C        NERR  - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE.
C                NERR MUST NOT BE ZERO.
C        LEVEL - ERROR CATEGORY.
C                =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR.
C                =1 MEANS THIS IS A RECOVERABLE ERROR.  (I.E., IT IS
C                   NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.)
C                =0 MEANS THIS IS A WARNING MESSAGE ONLY.
C                =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE
C                   PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY
C                   TIMES THIS CALL IS EXECUTED.
C        NI    - NUMBER OF INTEGER VALUES TO BE PRINTED. (O TO 2)
C        I1    - FIRST INTEGER VALUE.
C        I2    - SECOND INTEGER VALUE.
C        NR    - NUMBER OF REAL VALUES TO BE PRINTED. (0 TO 2)
C        R1    - FIRST REAL VALUE.
C        R2    - SECOND REAL VALUE.
C
C     EXAMPLES
C        CALL XERROR(29HSMOOTH -- NUM (=I1) WAS ZERO.,29,1,2,
C    1   1,NUM,0,0,0.,0.)
C        CALL XERRWV(54HQUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM
C    1 (R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN)
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     LATEST REVISION ---  19 MAR 1980
C
	DIMENSION MESSG(NMESSG),LUN(5)
C     GET FLAGS
	LKNTRL = J4SAVE(2,0,.FALSE.)
	MAXMES = J4SAVE(4,0,.FALSE.)
C     CHECK FOR VALID INPUT
	IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND.
     1    (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10
	   IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17)
	   CALL XERPRT(23HXERROR -- INVALID INPUT,23)
	   IF (LKNTRL.GT.0) CALL FDUMP
	   IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR.,
     1   29)
	   IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY)
	   CALL XERABT(23HXERROR -- INVALID INPUT,23)
	   RETURN
10	CONTINUE
C     RECORD MESSAGE
	JUNK = J4SAVE(1,NERR,.TRUE.)
	CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT)
C     LET USER OVERRIDE
	LFIRST = MESSG(1)
	LMESSG = NMESSG
	LERR = NERR
	LLEVEL = LEVEL
	CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL)
C     RESET TO ORIGINAL VALUES
	LMESSG = NMESSG
	LERR = NERR
	LLEVEL = LEVEL
	LKNTRL = MAX0(-2,MIN0(2,LKNTRL))
	MKNTRL = IABS(LKNTRL)
C     DECIDE WHETHER TO PRINT MESSAGE
	IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100
	IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES)))
     1.OR.((LLEVEL.EQ.0)   .AND.(KOUNT.GT.MAXMES))
     2.OR.((LLEVEL.EQ.1)   .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1))
     3.OR.((LLEVEL.EQ.2)   .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100
	   IF (LKNTRL.LE.0) GO TO 20
	      CALL XERPRT(1H ,1)
C           INTRODUCTION
	      IF (LLEVEL.EQ.(-1)) CALL XERPRT
     1(57HWARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.,57)
	      IF (LLEVEL.EQ.0) CALL XERPRT(13HWARNING IN...,13)
	      IF (LLEVEL.EQ.1) CALL XERPRT
     1      (23HRECOVERABLE ERROR IN...,23)
	      IF (LLEVEL.EQ.2) CALL XERPRT(17HFATAL ERROR IN...,17)
20	   CONTINUE
C        MESSAGE
	   CALL XERPRT(MESSG,LMESSG)
	   CALL XGETUA(LUN,NUNIT)
	   DO 50 KUNIT=1,NUNIT
	      IUNIT = LUN(KUNIT)
	      IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
	      IF (NI.GE.1) WRITE (IUNIT,22) I1
	      IF (NI.GE.2) WRITE (IUNIT,23) I2
	      IF (NR.GE.1) WRITE (IUNIT,24) R1
	      IF (NR.GE.2) WRITE (IUNIT,25) R2
22	      FORMAT (11X,21HIn above message, I1=,I10)
23	      FORMAT (11X,21HIn above message, I2=,I10)
24	      FORMAT (11X,21HIn above message, R1=,E20.10)
25	      FORMAT (11X,21HIn above message, R2=,E20.10)
	      IF (LKNTRL.LE.0) GO TO 40
C              ERROR NUMBER
	         WRITE (IUNIT,30) LERR
30	         FORMAT (15H Error number =,I10)
40	      CONTINUE
50	   CONTINUE
C        TRACE-BACK
	   IF (LKNTRL.GT.0) CALL FDUMP
100	CONTINUE
	IFATAL = 0
	IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2)))
     1IFATAL = 1
C     QUIT HERE IF MESSAGE IS NOT FATAL
	IF (IFATAL.LE.0) RETURN
	IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120
C        PRINT REASON FOR ABORT
	   IF (LLEVEL.EQ.1) CALL XERPRT
     1   (35HJOB ABORT DUE TO UNRECOVERED ERROR.,35)
	   IF (LLEVEL.EQ.2) CALL XERPRT
     1   (29HJOB ABORT DUE TO FATAL ERROR.,29)
C        PRINT ERROR SUMMARY
	   CALL XERSAV(1H ,-1,0,0,KDUMMY)
120	CONTINUE
C     ABORT
	IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0
	CALL XERABT(MESSG,LMESSG)
	RETURN
	END

	SUBROUTINE XERPRT(MESSG,NMESSG)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        PRINT THE HOLLERITH MESSAGE IN MESSG, OF LENGTH MESSG,
C        ON EACH FILE INDICATED BY XGETUA.
C        THIS VERSION PRINTS EXACTLY THE RIGHT NUMBER OF CHARACTERS,
C        NOT A NUMBER OF WORDS, AND THUS SHOULD WORK ON MACHINES
C        WHICH DO NOT BLANK FILL THE LAST WORD OF THE HOLLERITH.
C
C     RON JONES, JUNE 1980
C
	INTEGER F(10),G(14),LUN(5)
	DIMENSION MESSG(NMESSG)
	DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10)
     1   / 1H( ,1H1 ,1HX ,1H, ,1H  ,1H  ,1HA ,1H  ,1H  ,1H) /
	DATA G(1),G(2),G(3),G(4),G(5),G(6),G(7),G(8),G(9),G(10)
     1   / 1H( ,1H1 ,1HX ,1H  ,1H  ,1H  ,1H  ,1H  ,1H  ,1H  /
	DATA G(11),G(12),G(13),G(14)
     1   / 1H   ,1H   ,1H   ,1H)  /
	DATA LA/1HA/,LCOM/1H,/,LBLANK/1H /
C     PREPARE FORMAT FOR WHOLE LINES
	NCHAR = I1MACH(6)
	NFIELD = 72/NCHAR
	CALL S88FMT(2,NFIELD,F(5))
	CALL S88FMT(2,NCHAR,F(8))
C     PREPARE FORMAT FOR LAST, PARTIAL LINE, IF NEEDED
	NCHARL = NFIELD*NCHAR
	NLINES = NMESSG/NCHARL
	NWORD  = NLINES*NFIELD
	NCHREM = NMESSG - NLINES*NCHARL
	IF (NCHREM.LE.0) GO TO 40
	   DO 10 I=4,13
10	      G(I) = LBLANK
	   NFIELD = NCHREM/NCHAR
	   IF (NFIELD.LE.0) GO TO 20
C        PREPARE WHOLE WORD FIELDS
	      G(4) = LCOM
	      CALL S88FMT(2,NFIELD,G(5))
	      G(7) = LA
	      CALL S88FMT(2,NCHAR,G(8))
20	   CONTINUE
	   NCHLST = MOD(NCHREM,NCHAR)
	   IF (NCHLST.LE.0) GO TO 30
C        PREPARE PARTIAL WORD FIELD
	      G(10) = LCOM
	      G(11) = LA
	      CALL S88FMT(2,NCHLST,G(12))
30	   CONTINUE
40	CONTINUE
C     PRINT THE MESSAGE
	NWORD1 = NWORD+1
	NWORD2 = (NMESSG+NCHAR-1)/NCHAR
	CALL XGETUA(LUN,NUNIT)
	DO 50 KUNIT = 1,NUNIT
	   IUNIT = LUN(KUNIT)
	   IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
	   IF (NWORD.GT.0) WRITE (IUNIT,F) (MESSG(I),I=1,NWORD)
	   IF (NCHREM.GT.0) WRITE (IUNIT,G) (MESSG(I),I=NWORD1,NWORD2)
50	CONTINUE
	RETURN
	END

	SUBROUTINE FDUMP
	IMPLICIT INTEGER*4 (I-N)
C
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   *                 ISSUED BY SANDIA LABORATORIES,
C   *                   A PRIME CONTRACTOR TO THE
C   *                UNITED STATES DEPARTMENT OF ENERGY
C   * * * * * * * * * * * * * *  NOTICE   * * * * * * * * * * * * * * *
C   * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
C   * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
C   * UNITED STATES DEPARTMENT OF ENERGY
C   * NOR ANY OF THEIR EMPLOYEES, NOR ANY OF THEIR CONTRACTORS,
C   * SUBCONTRACTORS, OR THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS
C   * OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY
C   * FOR THE ACCURACY, COMPLETENESS OR USEFULNESS OF ANY INFORMATION,
C   * APPARATUS, PRODUCT OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS
C   * USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     ABSTRACT
C        ***NOTE*** MACHINE DEPENDENT ROUTINE
C        FDUMP IS INTENDED TO BE REPLACED BY A LOCALLY WRITTEN
C        VERSION WHICH PRODUCES A SYMBOLIC DUMP.  FAILING THIS,
C        IT SHOULD BE REPLACED BY A VERSION WHICH PRINTS THE
C        SUBPROGRAM NESTING LIST.  NOTE THAT THIS DUMP MUST BE
C        PRINTED ON EACH OF UP TO FIVE FILES, AS INDICATED BY THE
C        XGETUA ROUTINE.  SEE XSETUA AND XGETUA FOR DETAILS.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     LATEST REVISION ---  23 MAY 1979
	RETURN
	END

	SUBROUTINE S88FMT(N,IVALUE,IFMT)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        S88FMT REPLACES IFMT(1), ... ,IFMT(N) WITH THE
C        CHARACTERS CORRESPONDING TO THE N LEAST SIGNIFICANT
C        DIGITS OF IVALUE.
C
C     TAKEN FROM THE BELL LABORATORIES PORT LIBRARY ERROR HANDLER
C     LATEST REVISION ---  7 JUNE 1978
C
	DIMENSION IFMT(N),IDIGIT(10)
	DATA IDIGIT(1),IDIGIT(2),IDIGIT(3),IDIGIT(4),IDIGIT(5),
     1     IDIGIT(6),IDIGIT(7),IDIGIT(8),IDIGIT(9),IDIGIT(10)
     2     /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
	NT = N
	IT = IVALUE
10	   IF (NT .EQ. 0) RETURN
	   INDEX = MOD(IT,10)
	   IFMT(NT) = IDIGIT(INDEX+1)
	   IT = IT/10
	   NT = NT - 1
	   GO TO 10
	END

	SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        RECORD THAT THIS ERROR OCCURRED.
C
C     DESCRIPTION OF PARAMETERS
C     --INPUT--
C       MESSG, NMESSG, NERR, LEVEL ARE AS IN XERROR,
C       EXCEPT THAT WHEN NMESSG=0 THE TABLES WILL BE
C       DUMPED AND CLEARED, AND WHEN NMESSG IS LESS THAN ZERO THE
C       TABLES WILL BE DUMPED AND NOT CLEARED.
C     --OUTPUT--
C       ICOUNT WILL BE THE NUMBER OF TIMES THIS MESSAGE HAS
C       BEEN SEEN, OR ZERO IF THE TABLE HAS OVERFLOWED AND
C       DOES NOT CONTAIN THIS MESSAGE SPECIFICALLY.
C       WHEN NMESSG=0, ICOUNT WILL NOT BE ALTERED.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     LATEST REVISION ---  19 MAR 1980
C
	INTEGER F(17),LUN(5)
	DIMENSION MESSG(NMESSG)
	DIMENSION MESTAB(10),NERTAB(10),LEVTAB(10),KOUNT(10)
C     NEXT THREE DATA STATEMENTS ARE NEEDED MERELY TO SATISFY
C     CERTAIN CONVENTIONS FOR COMPILERS WHICH DYNAMICALLY
C     ALLOCATE STORAGE.
	DATA MESTAB(1),MESTAB(2),MESTAB(3),MESTAB(4),MESTAB(5),
     1     MESTAB(6),MESTAB(7),MESTAB(8),MESTAB(9),MESTAB(10)
     2     /0,0,0,0,0,0,0,0,0,0/
	DATA NERTAB(1),NERTAB(2),NERTAB(3),NERTAB(4),NERTAB(5),
     1     NERTAB(6),NERTAB(7),NERTAB(8),NERTAB(9),NERTAB(10)
     2     /0,0,0,0,0,0,0,0,0,0/
	DATA LEVTAB(1),LEVTAB(2),LEVTAB(3),LEVTAB(4),LEVTAB(5),
     1     LEVTAB(6),LEVTAB(7),LEVTAB(8),LEVTAB(9),LEVTAB(10)
     2     /0,0,0,0,0,0,0,0,0,0/
C     NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK
C     ERROR TABLE INITIALLY
	DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5),
     1     KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10)
     2     /0,0,0,0,0,0,0,0,0,0/
	DATA KOUNTX/0/
C     NEXT DATA STATEMENT SETS UP OUTPUT FORMAT
	DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10),
     1     F(11),F(12),F(13),F(14),F(15),F(16),F(17)
     2     /1H( ,1H1 ,1HX ,1H, ,1HA ,1H  ,1H  ,1H, ,1HI ,1H  ,
     3      1H  ,1H, ,1H2 ,1HI ,1H1 ,1H0 ,1H) /
	IF (NMESSG.GT.0) GO TO 80
C     DUMP THE TABLE
	   IF (KOUNT(1).EQ.0) RETURN
C        PREPARE FORMAT
	   NCHAR = I1MACH(6)
	   CALL S88FMT(2,NCHAR,F(6))
	   NCOL = 20 - NCHAR
	   CALL S88FMT(2,NCOL,F(10))
C        PRINT TO EACH UNIT
	   CALL XGETUA(LUN,NUNIT)
	   DO 60 KUNIT=1,NUNIT
	      IUNIT = LUN(KUNIT)
	      IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
C           PRINT TABLE HEADER
	      WRITE (IUNIT,10)
10	      FORMAT (32H0          Error Message Summary/
     1              41H First Word      Nerr     Level     Count)
C           PRINT BODY OF TABLE
	      DO 20 I=1,10
	         IF (KOUNT(I).EQ.0) GO TO 30
	         WRITE (IUNIT,F) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I)
20	      CONTINUE
30	      CONTINUE
C           PRINT NUMBER OF OTHER ERRORS
	      IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX
40	      FORMAT (41H0Other errors not individually tabulated=,I10)
	      WRITE (IUNIT,50)
50	      FORMAT (1X)
60	   CONTINUE
	   IF (NMESSG.LT.0) RETURN
C        CLEAR THE ERROR TABLES
	   DO 70 I=1,10
70	      KOUNT(I) = 0
	   KOUNTX = 0
	   RETURN
80	CONTINUE
C     PROCESS A MESSAGE...
C     SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
C     OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
	DO 90 I=1,10
	   II = I
	   IF (KOUNT(I).EQ.0) GO TO 110
	   IF (MESSG(1).NE.MESTAB(I)) GO TO 90
	   IF (NERR.NE.NERTAB(I)) GO TO 90
	   IF (LEVEL.NE.LEVTAB(I)) GO TO 90
	   GO TO 100
90	CONTINUE
C     THREE POSSIBLE CASES...
C     TABLE IS FULL
	   KOUNTX = KOUNTX+1
	   ICOUNT = 1
	   RETURN
C     MESSAGE FOUND IN TABLE
100	   KOUNT(II) = KOUNT(II) + 1
	   ICOUNT = KOUNT(II)
	   RETURN
C     EMPTY SLOT FOUND FOR NEW MESSAGE
110	   MESTAB(II) = MESSG(1)
	   NERTAB(II) = NERR
	   LEVTAB(II) = LEVEL
	   KOUNT(II)  = 1
	   ICOUNT = 1
	   RETURN
	END

	SUBROUTINE XERABT(MESSG,NMESSG)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        ***NOTE*** MACHINE DEPENDENT ROUTINE
C        XERABT ABORTS THE EXECUTION OF THE PROGRAM.
C        THE ERROR MESSAGE CAUSING THE ABORT IS GIVEN IN THE CALLING
C        SEQUENCE IN CASE ONE NEEDS IT FOR PRINTING ON A DAYFILE,
C        FOR EXAMPLE.
C
C     DESCRIPTION OF PARAMETERS
C        MESSG AND NMESSG ARE AS IN XERROR, EXCEPT THAT NMESSG MAY
C        BE ZERO, IN WHICH CASE NO MESSAGE IS BEING SUPPLIED.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     LATEST REVISION ---  7 JUNE 1978
C
	STOP
	END

	FUNCTION J4SAVE(IWHICH,IVALUE,ISET)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        J4SAVE SAVES AND RECALLS SEVERAL GLOBAL VARIABLES NEEDED
C        BY THE LIBRARY ERROR HANDLING ROUTINES.
C
C     DESCRIPTION OF PARAMETERS
C      --INPUT--
C        IWHICH - INDEX OF ITEM DESIRED.
C                 = 1 REFERS TO CURRENT ERROR NUMBER.
C                 = 2 REFERS TO CURRENT ERROR CONTROL FLAG.
C                 = 3 REFERS TO CURRENT UNIT NUMBER TO WHICH ERROR
C                     MESSAGES ARE TO BE SENT.  (0 MEANS USE STANDARD.)
C                 = 4 REFERS TO THE MAXIMUM NUMBER OF TIMES ANY
C                     MESSAGE IS TO BE PRINTED (AS SET BY XERMAX).
C                 = 5 REFERS TO THE TOTAL NUMBER OF UNITS TO WHICH
C                     EACH ERROR MESSAGE IS TO BE WRITTEN.
C                 = 6 REFERS TO THE 2ND UNIT FOR ERROR MESSAGES
C                 = 7 REFERS TO THE 3RD UNIT FOR ERROR MESSAGES
C                 = 8 REFERS TO THE 4TH UNIT FOR ERROR MESSAGES
C                 = 9 REFERS TO THE 5TH UNIT FOR ERROR MESSAGES
C        IVALUE - THE VALUE TO BE SET FOR THE IWHICH-TH PARAMETER,
C                 IF ISET IS .TRUE. .
C        ISET   - IF ISET=.TRUE., THE IWHICH-TH PARAMETER WILL BE
C                 GIVEN THE VALUE, IVALUE.  IF ISET=.FALSE., THE
C                 IWHICH-TH PARAMETER WILL BE UNCHANGED, AND IVALUE
C                 IS A DUMMY PARAMETER.
C      --OUTPUT--
C        THE (OLD) VALUE OF THE IWHICH-TH PARAMETER WILL BE RETURNED
C        IN THE FUNCTION VALUE, J4SAVE.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C     ADAPTED FROM BELL LABORATORIES PORT LIBRARY ERROR HANDLER
C     LATEST REVISION ---  23 MAY 1979
C
	LOGICAL ISET
	INTEGER IPARAM(9)
	DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,1,0,10/
	DATA IPARAM(5)/1/
	DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
	J4SAVE = IPARAM(IWHICH)
	IF (ISET) IPARAM(IWHICH) = IVALUE
	RETURN
	END

	SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        ALLOWS USER CONTROL OVER HANDLING OF INDIVIDUAL ERRORS.
C        JUST AFTER EACH MESSAGE IS RECORDED, BUT BEFORE IT IS
C        PROCESSED ANY FURTHER (I.E., BEFORE IT IS PRINTED OR
C        A DECISION TO ABORT IS MADE) A CALL IS MADE TO XERCTL.
C        IF THE USER HAS PROVIDED HIS OWN VERSION OF XERCTL, HE
C        CAN THEN OVERRIDE THE VALUE OF KONTROL USED IN PROCESSING
C        THIS MESSAGE BY REDEFINING ITS VALUE.
C        KONTRL MAY BE SET TO ANY VALUE FROM -2 TO 2.
C        THE MEANINGS FOR KONTRL ARE THE SAME AS IN XSETF, EXCEPT
C        THAT THE VALUE OF KONTRL CHANGES ONLY FOR THIS MESSAGE.
C        IF KONTRL IS SET TO A VALUE OUTSIDE THE RANGE FROM -2 TO 2,
C        IT WILL BE MOVED BACK INTO THAT RANGE.
C
C     DESCRIPTION OF PARAMETERS
C
C      --INPUT--
C        MESSG1 - THE FIRST WORD (ONLY) OF THE ERROR MESSAGE.
C        NMESSG - SAME AS IN THE CALL TO XERROR OR XERRWV.
C        NERR   - SAME AS IN THE CALL TO XERROR OR XERRWV.
C        LEVEL  - SAME AS IN THE CALL TO XERROR OR XERRWV.
C        KONTRL - THE CURRENT VALUE OF THE CONTROL FLAG AS SET
C                 BY A CALL TO XSETF.
C
C      --OUTPUT--
C        KONTRL - THE NEW VALUE OF KONTRL.  IF KONTRL IS NOT
C                 DEFINED, IT WILL REMAIN AT ITS ORIGINAL VALUE.
C                 THIS CHANGED VALUE OF CONTROL AFFECTS ONLY
C                 THE CURRENT OCCURRENCE OF THE CURRENT MESSAGE.
C
	RETURN
	END

	SUBROUTINE XGETUA(IUNIT,N)
	IMPLICIT INTEGER*4 (I-N)
C
C     ABSTRACT
C        XGETUA MAY BE CALLED TO DETERMINE THE UNIT NUMBER OR NUMBERS
C        TO WHICH ERROR MESSAGES ARE BEING SENT.
C        THESE UNIT NUMBERS MAY HAVE BEEN SET BY A CALL TO XSETUN,
C        OR A CALL TO XSETUA, OR MAY BE A DEFAULT VALUE.
C
C     DESCRIPTION OF PARAMETERS
C      --OUTPUT--
C        IUNIT - AN ARRAY OF ONE TO FIVE UNIT NUMBERS, DEPENDING
C                ON THE VALUE OF N.  A VALUE OF ZERO REFERS TO THE
C                DEFAULT UNIT, AS DEFINED BY THE I1MACH MACHINE
C                CONSTANT ROUTINE.  ONLY IUNIT(1),...,IUNIT(N) ARE
C                DEFINED BY XGETUA.  THE VALUES OF IUNIT(N+1),...,
C                IUNIT(5) ARE NOT DEFINED (FOR N.LT.5) OR ALTERED
C                IN ANY WAY BY XGETUA.
C        N     - THE NUMBER OF UNITS TO WHICH COPIES OF THE
C                ERROR MESSAGES ARE BEING SENT.  N WILL BE IN THE
C                RANGE FROM 1 TO 5.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C
	DIMENSION IUNIT(5)
	N = J4SAVE(5,0,.FALSE.)
	DO 30 I=1,N
	   INDEX = I+4
	   IF (I.EQ.1) INDEX = 3
	   IUNIT(I) = J4SAVE(INDEX,0,.FALSE.)
30	CONTINUE
	RETURN
	END

	INTEGER*4 FUNCTION I1MACH(I)
	IMPLICIT INTEGER*4 (I-N)
C
C  I/O UNIT NUMBERS.
C
C    I1MACH( 1) = THE STANDARD INPUT UNIT.
C
C    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
C
C    I1MACH( 3) = THE STANDARD PUNCH UNIT.
C
C    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
C
C  WORDS.
C
C    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
C
C    I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
C
C  INTEGERS.
C
C    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C
C               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C
C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
C
C    I1MACH( 7) = A, THE BASE.
C
C    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
C
C    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
C
C  FLOATING-POINT NUMBERS.
C
C    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C    BASE-B FORM
C
C               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
C
C    I1MACH(10) = B, THE BASE.
C
C  SINGLE-PRECISION
C
C    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C  DOUBLE-PRECISION
C
C    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C  WITH THE LOCAL OPERATING SYSTEM.
C
	INTEGER IMACH(16),OUTPUT
C
	EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA IMACH( 1) /    7 /
C     DATA IMACH( 2) /    2 /
C     DATA IMACH( 3) /    2 /
C     DATA IMACH( 4) /    2 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    4 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   33 /
C     DATA IMACH( 9) / Z1FFFFFFFF /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   24 /
C     DATA IMACH(12) / -256 /
C     DATA IMACH(13) /  255 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) / -256 /
C     DATA IMACH(16) /  255 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA IMACH( 1) /   5 /
C     DATA IMACH( 2) /   6 /
C     DATA IMACH( 3) /   7 /
C     DATA IMACH( 4) /   6 /
C     DATA IMACH( 5) /  48 /
C     DATA IMACH( 6) /   6 /
C     DATA IMACH( 7) /   2 /
C     DATA IMACH( 8) /  39 /
C     DATA IMACH( 9) / O0007777777777777 /
C     DATA IMACH(10) /   8 /
C     DATA IMACH(11) /  13 /
C     DATA IMACH(12) / -50 /
C     DATA IMACH(13) /  76 /
C     DATA IMACH(14) /  26 /
C     DATA IMACH(15) / -50 /
C     DATA IMACH(16) /  76 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA IMACH( 1) /   5 /
C     DATA IMACH( 2) /   6 /
C     DATA IMACH( 3) /   7 /
C     DATA IMACH( 4) /   6 /
C     DATA IMACH( 5) /  48 /
C     DATA IMACH( 6) /   6 /
C     DATA IMACH( 7) /   2 /
C     DATA IMACH( 8) /  39 /
C     DATA IMACH( 9) / O0007777777777777 /
C     DATA IMACH(10) /   8 /
C     DATA IMACH(11) /  13 /
C     DATA IMACH(12) / -50 /
C     DATA IMACH(13) /  76 /
C     DATA IMACH(14) /  26 /
C     DATA IMACH(15) / -32754 /
C     DATA IMACH(16) /  32780 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    7 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   60 /
C     DATA IMACH( 6) /   10 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   48 /
C     DATA IMACH( 9) / 00007777777777777777B /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   48 /
C     DATA IMACH(12) / -974 /
C     DATA IMACH(13) / 1070 /
C     DATA IMACH(14) /   96 /
C     DATA IMACH(15) / -927 /
C     DATA IMACH(16) / 1070 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA IMACH( 1) /   100 /
C     DATA IMACH( 2) /   101 /
C     DATA IMACH( 3) /   102 /
C     DATA IMACH( 4) /   101 /
C     DATA IMACH( 5) /    64 /
C     DATA IMACH( 6) /     8 /
C     DATA IMACH( 7) /     2 /
C     DATA IMACH( 8) /    63 /
C     DATA IMACH( 9) /  777777777777777777777B /
C     DATA IMACH(10) /     2 /
C     DATA IMACH(11) /    48 /
C     DATA IMACH(12) / -8192 /
C     DATA IMACH(13) /  8191 /
C     DATA IMACH(14) /    96 /
C     DATA IMACH(15) / -8192 /
C     DATA IMACH(16) /  8191 /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     DATA IMACH( 1) /   11 /
C     DATA IMACH( 2) /   12 /
C     DATA IMACH( 3) /    8 /
C     DATA IMACH( 4) /   10 /
C     DATA IMACH( 5) /   16 /
C     DATA IMACH( 6) /    2 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   15 /
C     DATA IMACH( 9) /32767 /
C     DATA IMACH(10) /   16 /
C     DATA IMACH(11) /    6 /
C     DATA IMACH(12) /  -64 /
C     DATA IMACH(13) /   63 /
C     DATA IMACH(14) /   14 /
C     DATA IMACH(15) /  -64 /
C     DATA IMACH(16) /   63 /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA IMACH( 1) /       5 /
C     DATA IMACH( 2) /       6 /
C     DATA IMACH( 3) /       0 /
C     DATA IMACH( 4) /       6 /
C     DATA IMACH( 5) /      24 /
C     DATA IMACH( 6) /       3 /
C     DATA IMACH( 7) /       2 /
C     DATA IMACH( 8) /      23 /
C     DATA IMACH( 9) / 8388607 /
C     DATA IMACH(10) /       2 /
C     DATA IMACH(11) /      23 /
C     DATA IMACH(12) /    -127 /
C     DATA IMACH(13) /     127 /
C     DATA IMACH(14) /      38 /
C     DATA IMACH(15) /    -127 /
C     DATA IMACH(16) /     127 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /   43 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    6 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -127 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   63 /
C     DATA IMACH(15) / -127 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA I1MACH( 1) /    5 /
C     DATA I1MACH( 2) /    6 /
C     DATA I1MACH( 3) /    4 /
C     DATA I1MACH( 4) /    1 /
C     DATA I1MACH( 5) /   16 /
C     DATA I1MACH( 6) /    2 /
C     DATA I1MACH( 7) /    2 /
C     DATA I1MACH( 8) /   15 /
C     DATA I1MACH( 9) / 32767 /
C     DATA I1MACH(10) /    2 /
C     DATA I1MACH(11) /   23 /
C     DATA I1MACH(12) / -128 /
C     DATA I1MACH(13) /  127 /
C     DATA I1MACH(14) /   39 /
C     DATA I1MACH(15) / -128 /
C     DATA I1MACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA I1MACH( 1) /    5 /
C     DATA I1MACH( 2) /    6 /
C     DATA I1MACH( 3) /    4 /
C     DATA I1MACH( 4) /    1 /
C     DATA I1MACH( 5) /   16 /
C     DATA I1MACH( 6) /    2 /
C     DATA I1MACH( 7) /    2 /
C     DATA I1MACH( 8) /   15 /
C     DATA I1MACH( 9) / 32767 /
C     DATA I1MACH(10) /    2 /
C     DATA I1MACH(11) /   23 /
C     DATA I1MACH(12) / -128 /
C     DATA I1MACH(13) /  127 /
C     DATA I1MACH(14) /   55 /
C     DATA I1MACH(15) / -128 /
C     DATA I1MACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA IMACH( 1) /   5 /
C     DATA IMACH( 2) /   6 /
C     DATA IMACH( 3) /   7 /
C     DATA IMACH( 4) /   6 /
C     DATA IMACH( 5) /  32 /
C     DATA IMACH( 6) /   4 /
C     DATA IMACH( 7) /   2 /
C     DATA IMACH( 8) /  31 /
C     DATA IMACH( 9) / Z7FFFFFFF /
C     DATA IMACH(10) /  16 /
C     DATA IMACH(11) /   6 /
C     DATA IMACH(12) / -64 /
C     DATA IMACH(13) /  63 /
C     DATA IMACH(14) /  14 /
C     DATA IMACH(15) / -64 /
C     DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    5 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / "377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   54 /
C     DATA IMACH(15) / -101 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    5 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / "377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   62 /
C     DATA IMACH(15) / -128 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN'S SUPPORTING
C     32-BIT INTEGER ARITHMETIC.
C
      DATA IMACH( 1) /    5 /
      DATA IMACH( 2) /    6 /
      DATA IMACH( 3) /    5 /
      DATA IMACH( 4) /    6 /
      DATA IMACH( 5) /   32 /
      DATA IMACH( 6) /    4 /
      DATA IMACH( 7) /    2 /
      DATA IMACH( 8) /   31 /
      DATA IMACH( 9) / 2147483647 /
      DATA IMACH(10) /    2 /
      DATA IMACH(11) /   24 /
      DATA IMACH(12) / -127 /
      DATA IMACH(13) /  127 /
      DATA IMACH(14) /   56 /
      DATA IMACH(15) / -127 /
      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN'S SUPPORTING
C     16-BIT INTEGER ARITHMETIC.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   16 /
C     DATA IMACH( 6) /    2 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   15 /
C     DATA IMACH( 9) / 32767 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   24 /
C     DATA IMACH(12) / -127 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   56 /
C     DATA IMACH(15) / -127 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    1 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    4 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) /-1024 /
C     DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
C
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    1 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    6 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) /-1024 /
C     DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR THE VAX-11/780
C
C	DATA IMACH( 1) /    5 /
C	DATA IMACH( 2) /    6 /
C	DATA IMACH( 3) /    5 /
C	DATA IMACH( 4) /    6 /
C	DATA IMACH( 5) /   32 /
C	DATA IMACH( 6) /    4 /
C	DATA IMACH( 7) /    2 /
C	DATA IMACH( 8) /   31 /
C	DATA IMACH( 9) / 2147483647 /
C	DATA IMACH(10) /    2 /
C	DATA IMACH(11) /   24 /
C	DATA IMACH(12) / -127 /
C	DATA IMACH(13) /  127 /
C	DATA IMACH(14) /   56 /
C	DATA IMACH(15) / -127 /
C	DATA IMACH(16) /  127 /
	IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
C
	I1MACH=IMACH(I)
	RETURN
C
10	WRITE(OUTPUT,9000)
9000	FORMAT(39H1Error    1 in I1MACH - I out of bounds)
C
C     CALL FDUMP
C
	STOP
C
	END

	SUBROUTINE WNNLS(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,
     * IWORK, WORK)
	IMPLICIT INTEGER*4 (I-N)
C     DIMENSION W(MDW,N+1),PRGOPT(*),X(N),IWORK(M+N),WORK(M+5*N)
C
C     WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES,
C     AND R.J. HANSON, SANDIA LABORATORIES.
C
C     ABSTRACT
C
C     THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES
C     PROBLEM.  SUPPOSE THERE ARE GIVEN MATRICES E AND A OF
C     RESPECTIVE DIMENSIONS ME BY N AND MA BY N, AND VECTORS F
C     AND B OF RESPECTIVE LENGTHS ME AND MA.  THIS SUBROUTINE
C     SOLVES THE PROBLEM
C
C               EX = F, (EQUATIONS TO BE EXACTLY SATISFIED)
C
C               AX = B, (EQUATIONS TO BE APPROXIMATELY SATISFIED,
C                        IN THE LEAST SQUARES SENSE)
C
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE
C
C     ANY VALUES ME.GE.0, MA.GE.0 AND 0.LE. L .LE.N ARE PERMITTED.
C
C     THE PROBLEM IS REPOSED AS PROBLEM WNNLS
C
C               (WT*E)X = (WT*F)
C               (   A)    (   B), (LEAST SQUARES)
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE.
C
C     THE SUBPROGRAM CHOOSES THE HEAVY WEIGHT (OR PENALTY PARAMETER) WT.
C
C     THE PARAMETERS FOR WNNLS ARE
C
C     INPUT..
C
C     W(*,*),MDW,  THE ARRAY W(*,*) IS DOUBLE SUBSCRIPTED WITH FIRST
C     ME,MA,N,L    DIMENSIONING PARAMETER EQUAL TO MDW.  FOR THIS
C                  DISCUSSION LET US CALL M = ME + MA.  THEN MDW
C                  MUST SATISFY MDW.GE.M.  THE CONDITION MDW.LT.M
C                  IS AN ERROR.
C
C                  THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS
C
C                       (E  F)
C                       (A  B)
C
C                  IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1
C                  RESPECTIVELY.  COLUMNS 1,...,L CORRESPOND TO
C                  UNCONSTRAINED VARIABLES X(1),...,X(L).  THE
C                  REMAINING VARIABLES ARE CONSTRAINED TO BE
C                  NONNEGATIVE.  THE CONDITION L.LT.0 .OR. L.GT.N IS
C                  AN ERROR.
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=6
C                         SCALE THE NONZERO COLUMNS OF THE
C                         ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  TO HAVE LENGTH ONE.  THE DATA SET FOR
C                  THIS OPTION IS A SINGLE VALUE.  IT MUST
C                  BE NONZERO IF UNIT LENGTH COLUMN SCALING IS
C                  DESIRED.
C
C                  KEY=7
C                         SCALE COLUMNS OF THE ENTIRE DATA MATRIX
C                  (E)
C                  (A)
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=8
C                         CHANGE THE RANK DETERMINATION TOLERANCE 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 PROVIDE
C                  A DIAGONAL MATRIX TO SCALE THE PROBLEM
C                  MATRIX AND CHANGE THE TOLERANCE USED FOR
C                  DETERMINING LINEAR DEPENDENCE OF DROPPED COL
C                  VECTORS.  FOR THESE OPTIONS THE DIMENSIONS OF
C                  PRGOPT(*) MUST BE AT LEAST N+6.  THE FORTRAN
C                  STATEMENTS DEFINING THESE OPTIONS WOULD
C                  BE AS FOLLOWS.
C
C                  PRGOPT(1)=N+3 (LINK TO ENTRY N+3 IN PRGOPT(*))
C                  PRGOPT(2)=7 (USER-PROVIDED SCALING KEY)
C
C                  CALL SCOPY(N,D,1,PRGOPT(3),1) (COPY THE N
C                  SCALING FACTORS FROM A USER ARRAY CALLED D(*)
C                  INTO PRGOPT(3)-PRGOPT(N+2))
C
C                  PRGOPT(N+3)=N+6 (LINK TO ENTRY N+6 OF PRGOPT(*))
C                  PRGOPT(N+4)=8 (LINEAR DEPENDENCE TOLERANCE KEY)
C                  PRGOPT(N+5)=... (NEW VALUE OF THE TOLERANCE)
C
C                  PRGOPT(N+6)=1 (NO MORE OPTIONS TO CHANGE)
C
C
C     OUTPUT..
C
C     X(*)         AN ARRAY DIMENSIONED AT LEAST N, WHICH WILL
C                  CONTAIN THE N COMPONENTS OF THE SOLUTION VECTOR
C                  ON OUTPUT.
C
C     RNORM        THE RESIDUAL NORM OF THE SOLUTION.  THE VALUE OF
C                  RNORM CONTAINS THE RESIDUAL VECTOR LENGTH OF THE
C                  EQUALITY CONSTRAINTS AND LEAST SQUARES EQUATIONS.
C
C     MODE         THE VALUE OF MODE INDICATES THE SUCCESS OR FAILURE
C                  OF THE SUBPROGRAM.
C
C                  MODE = 0  SUBPROGRAM COMPLETED SUCCESSFULLY.
C
C                       = 1  MAX. NUMBER OF ITERATIONS (EQUAL TO
C                            3*(N-L)) EXCEEDED. NEARLY ALL PROBLEMS
C                            SHOULD COMPLETE IN FEWER THAN THIS
C                            NUMBER OF ITERATIONS. AN APPROXIMATE
C                            SOLUTION AND ITS CORRESPONDING RESIDUAL
C                            VECTOR LENGTH ARE IN X(*) AND RNORM.
C
C                       = 2  USAGE ERROR OCCURRED.  THE OFFENDING
C                            CONDITION IS NOTED WITH THE ERROR
C                            PROCESSING SUBPROGRAM, XERROR( ).
C
C     USER-DESIGNATED
C     WORKING ARRAYS..
C
C     WORK(*)      A REAL-VALUED WORKING ARRAY OF LENGTH AT LEAST
C                  M + 5*N.
C
C     IWORK(*)     AN INTEGER-VALUED WORKING ARRAY OF LENGTH AT LEAST
C                  M+N.
C
C     SUBROUTINES CALLED
C     BY WNNLS( ) AND ITS
C     SUBPROGRAMS WNLSM( ) AND WNLIT( )
C
C     SROTMG,SROTM, SUBROUTINES FROM THE BLAS PACKAGE.
C     SAXPY,SASUM,  THESE ARE FOUND IN THE SANDIA MATH. LIBRARY,
C     SCOPY,SNRM2,  AND ARE ALSO DOCUMENTED IN THE SANDIA
C     SSWAP,ISAMAX, TECH. REPT., SAND77-0898.
C     SSCAL
C
C     H12           SUBROUTINE TO CONSTRUCT AND APPLY A
C                   HOUSEHOLDER TRANSFORMATION.  LOCATED IN SANDIA
C                   MATH. LIBRARY.
C
C     XERROR        FROM SLATEC ERROR PROCESSING PACKAGE.
C                   THIS IS DOCUMENTED IN SANDIA TECH. REPT.,
C                   SAND78-1189.
C
C     REFERENCES
C
C     1. SOLVING LEAST SQUARES PROBLEMS, BY C.L. LAWSON
C        AND R.J. HANSON.  PRENTICE-HALL, INC. (1974).
C
C     2. BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, BY
C        C.L. LAWSON, R.J. HANSON, D.R. KINCAID, AND F.T. KROGH.
C        SUBMITTED TO TOMS, JULY 1977.  ALSO AVAILABLE AS
C        SANDIA TECHNICAL REPORT NO. SAND77-0898.
C
C     3. AN ALGORITHM FOR LINEAR LEAST SQUARES WITH EQUALITY
C        AND NONNEGATIVITY CONSTRAINTS, BY K.H. HASKELL AND
C        R.J. HANSON.  AVAILABLE AS SANDIA TECHNICAL REPORT NO.
C        SAND77-0552.
C
C     4. SLATEC COMMON MATH. LIBRARY ERROR HANDLING
C        PACKAGE.  BY R. E. JONES.  AVAILABLE AS SANDIA
C        TECHNICAL REPORT SAND78-1189.
C
	DIMENSION W(MDW,1), PRGOPT(1), X(1), IWORK(1), WORK(1)
C     REVISED SEP. 6, 1979.
C
	MODE = 0
	IF (MA+ME.LE.0 .OR. N.LE.0) RETURN
	IF (.NOT.(MDW.LT.ME+MA)) GO TO 10
	NERR = 1
	IOPT = 1
C
C     THE OUTPUT UNIT FOR ERROR MESSAGES WITH XERROR( )
C     IS UNIT 6.
	CALL XERROR(44HWNNLS( ), THE VALUE MDW.LT.ME+MA IS AN ERROR, 44,
     * NERR, IOPT)
	MODE = 2
	RETURN
10	IF (0.LE.L .AND. L.LE.N) GO TO 20
	NERR = 2
	IOPT = 1
	CALL XERROR(39HWNNLS( ), L.LE.0.AND.L.LE.N IS REQUIRED, 39, NERR,
     * IOPT)
	MODE = 2
	RETURN
C
C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS
C     REQUIRED BY THE MAIN SUBROUTINE WNLSM( ).
C
20	L1 = N + 1
	L2 = L1 + N
	L3 = L2 + ME + MA
	L4 = L3 + N
	L5 = L4 + N
C
	CALL WNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK,
     * IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3), WORK(L4),
     * WORK(L5))
	RETURN
	END

	SUBROUTINE WNLSM(W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE,
     * IPIVOT, ITYPE, WD, H, SCALE, Z, TEMP, D)
	IMPLICIT INTEGER*4 (I-N)
C
C     THIS IS A COMPANION SUBPROGRAM TO WNNLS( ).
C     THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE
C     USAGE INSTRUCTIONS.
C
C     WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES,
C     WITH THE HELP OF R.J. HANSON, SANDIA LABORATORIES,
C     DECEMBER 1976 - JANUARY 1978.
C     REVISED SEP. 6, 1979.
C
C     IN ADDITION TO THE PARAMETERS DISCUSSED IN THE PROLOGUE TO
C     SUBROUTINE WNNLS, THE FOLLOWING WORK ARRAYS ARE USED IN
C     SUBROUTINE WNLSM  (THEY ARE PASSED THROUGH THE CALLING
C     SEQUENCE FROM WNNLS FOR PURPOSES OF VARIABLE DIMENSIONING).
C     THEIR CONTENTS WILL IN GENERAL BE OF NO INTEREST TO THE USER.
C
C         IPIVOT(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         PIVOTING INFORMATION FOR THE COLS OF W(*,*).
C
C         ITYPE(*)
C            AN ARRAY OF LENGTH M WHICH IS USED TO KEEP TRACK
C         OF THE CLASSIFICATION OF THE EQUATIONS.  ITYPE(I)=0
C         DENOTES EQUATION I AS AN EQUALITY CONSTRAINT.
C         ITYPE(I)=1 DENOTES EQUATION I AS A LEAST SQUARES
C         EQUATION.
C
C         WD(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         DUAL SOLUTION VECTOR.
C
C         H(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         PIVOT SCALARS OF THE HOUSEHOLDER TRANSFORMATIONS PERFORMED
C         IN THE CASE KRANK.LT.L.
C
C         SCALE(*)
C            AN ARRAY OF LENGTH M WHICH IS USED BY THE SUBROUTINE
C         TO STORE THE DIAGONAL MATRIX OF WEIGHTS.
C         THESE ARE USED TO APPLY THE MODIFIED GIVENS
C         TRANSFORMATIONS.
C
C         Z(*),TEMP(*)
C            WORKING ARRAYS OF LENGTH N.
C
C         D(*)
C            AN ARRAY OF LENGTH N THAT CONTAINS THE
C         COLUMN SCALING FOR THE MATRIX (E).
C                                       (A)
C
C     SUBROUTINE WNLSM (W,MDW,MME,MA,N,L,PRGOPT,X,RNORM,MODE,
C    1                  IPIVOT,ITYPE,WD,H,SCALE,Z,TEMP,D)
	DIMENSION W(MDW,1), X(1)
	DIMENSION WD(1), H(1), IPIVOT(1), ITYPE(1), SCALE(1)
	DIMENSION Z(1), TEMP(1), PRGOPT(1), D(1)
	DIMENSION SPARAM(5), IDOPE(8), DOPE(4)
	LOGICAL HITCON, FEASBL, DONE, POS
	DATA ZERO /0.E0/, ONE /1.E0/, TWO /2.E0/, EPS /0.E0/
C
C
C
C     INITIALIZE-VARIABLES
	ASSIGN 10 TO I99998
	GO TO 180
C
C     PERFORM INITIAL TRIANGULARIZATION IN THE SUBMATRIX
C     CORRESPONDING TO THE UNCONSTRAINED VARIABLES USING
C     THE PROCEDURE INITIALLY-TRIANGULARIZE.
10	ASSIGN 20 TO I99995
	GO TO 280
C
C     PERFORM WNNLS ALGORITHM USING THE FOLLOWING STEPS.
C
C     UNTIL(DONE)
C
C        COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
C        WHEN (HITCON) ADD-CONSTRAINTS
C
C        ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
C
C        FIN
C
C     COMPUTE-FINAL-SOLUTION
C
20	IF (DONE) GO TO 80
C
	ASSIGN 30 TO I99991
	GO TO 300
C
C     COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
30	IF (.NOT.(HITCON)) GO TO 50
	ASSIGN 40 TO I99986
	GO TO 370
40	GO TO 70
C
C     WHEN (HITCON) ADD-CONSTRAINTS
C
50	ASSIGN 60 TO I99983
	GO TO 640
60	CONTINUE
C
C     ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
C
70	GO TO 20
C
80	ASSIGN 90 TO I99980
	GO TO 1000
C
C     COMPUTE-FINAL-SOLUTION
C
90	RETURN
100	CONTINUE
C
C     TO PROCESS-OPTION-VECTOR
	FAC = 1.E-4
C
C     THE NOMINAL TOLERANCE USED IN THE CODE,
	TAU = SQRT(EPS)
C
C     THE NOMINAL COLUMN SCALING USED IN THE CODE IS
C     THE IDENTITY SCALING.
	D(1) = ONE
	CALL SCOPY(N, D, 0, D, 1)
C
C     DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE.
	NOPT = 1000
C
C     DEFINE BOUND FOR POSITIVE VALUE OF LINK.
	NLINK = 100000
	NTIMES = 0
	LAST = 1
	LINK = PRGOPT(1)
	IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 110
	NERR = 3
	IOPT = 1
	CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR,
     * IOPT)
	MODE = 2
	RETURN
110	IF (.NOT.(LINK.GT.1)) GO TO 160
	NTIMES = NTIMES + 1
	IF (.NOT.(NTIMES.GT.NOPT)) GO TO 120
	NERR = 3
	IOPT = 1
	CALL XERROR(
     * 53HWNNLS( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 53,
     * NERR, IOPT)
	MODE = 2
	RETURN
120	KEY = PRGOPT(LAST+1)
	IF (.NOT.(KEY.EQ.6 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 140
	DO 130 J=1,N
	   T = SNRM2(M,W(1,J),1)
	   IF (T.NE.ZERO) T = ONE/T
	   D(J) = T
130	CONTINUE
140	IF (KEY.EQ.7) CALL SCOPY(N, PRGOPT(LAST+2), 1, D, 1)
	IF (KEY.EQ.8) TAU = AMAX1(EPS,PRGOPT(LAST+2))
	NEXT = PRGOPT(LINK)
	IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 150
	NERR = 3
	IOPT = 1
	CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR,
     * IOPT)
	MODE = 2
	RETURN
150	LAST = LINK
	LINK = NEXT
	GO TO 110
160	DO 170 J=1,N
	   CALL SSCAL(M, D(J), W(1,J), 1)
170	CONTINUE
	GO TO 1260
180	CONTINUE
C
C     TO INITIALIZE-VARIABLES
C
C     EPS IS THE PRECISION FOR THE PARTICULAR MACHINE
C     BEING USED.  THIS LOGIC AVOIDS RECOMPUTING IT EVERY ENTRY.
	IF (.NOT.(EPS.EQ.ZERO)) GO TO 210
	EPS = ONE
190	IF (ONE+EPS.EQ.ONE) GO TO 200
	EPS = EPS/TWO
	GO TO 190
200	EPS = EPS*TWO
210	M = MA + MME
	ME = MME
	MEP1 = ME + 1
	ASSIGN 220 TO I99977
	GO TO 100
C
C     PROCESS-OPTION-VECTOR
220	DONE = .FALSE.
	ITER = 0
	ITMAX = 3*(N-L)
	MODE = 0
	LP1 = L + 1
	NSOLN = L
	NSP1 = NSOLN + 1
	NP1 = N + 1
	NM1 = N - 1
	L1 = MIN0(M,L)
C
C     COMPUTE SCALE FACTOR TO APPLY TO EQUAL. CONSTRAINT EQUAS.
	DO 230 J=1,N
	   WD(J) = SASUM(M,W(1,J),1)
230	CONTINUE
	IMAX = ISAMAX(N,WD,1)
	EANORM = WD(IMAX)
	BNORM = SASUM(M,W(1,NP1),1)
	RATNRM = BNORM/EANORM
	RATNSQ = RATNRM**2
	ALAMDA = EANORM/(EPS*FAC)
	EPENSQ = EPS*EANORM**2
C
C     DEFINE SCALING DIAG MATRIX FOR MOD GIVENS USAGE AND
C     CLASSIFY EQUATION TYPES.
	ALSQ = ALAMDA**2
	DO 260 I=1,M
C
C     WHEN EQU I IS HEAVILY WEIGHTED ITYPE(I)=0, ELSE ITYPE(I)=1.
	   IF (.NOT.(I.LE.ME)) GO TO 240
	   T = ALSQ
	   ITEMP = 0
	   GO TO 250
240	   T = ONE
	   ITEMP = 1
250	   SCALE(I) = T
	   ITYPE(I) = ITEMP
260	CONTINUE
C
C     SET THE SOLN VECTOR X(*) TO ZERO AND THE COL INTERCHANGE
C     MATRIX TO THE IDENTITY.
	X(1) = ZERO
	CALL SCOPY(N, X, 0, X, 1)
	DO 270 I=1,N
	   IPIVOT(I) = I
270	CONTINUE
	GO TO 1230
280	CONTINUE
C
C     TO INITIALLY-TRIANGULARIZE
C
C     SET FIRST L COMPS. OF DUAL VECTOR TO ZERO BECAUSE
C     THESE CORRESPOND TO THE UNCONSTRAINED VARIABLES.
	IF (.NOT.(L.GT.0)) GO TO 290
	WD(1) = ZERO
	CALL SCOPY(L, WD, 0, WD, 1)
C
C     THE ARRAYS IDOPE(*) AND DOPE(*) ARE USED TO PASS
C     INFORMATION TO WNLIT().  THIS WAS DONE TO AVOID
C     A LONG CALLING SEQUENCE OR THE USE OF COMMON.
290	IDOPE(1) = ME
	IDOPE(2) = MEP1
C     IDOPE(3)=KRANK
C     IDOPE(4)=KRP1
	IDOPE(5) = NSOLN
	IDOPE(6) = NIV
	IDOPE(7) = NIV1
	IDOPE(8) = L1
C
	DOPE(1) = ALSQ
	DOPE(2) = EPENSQ
	DOPE(3) = FAC
	DOPE(4) = TAU
	CALL WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM,
     * IDOPE, DOPE, DONE)
	ME = IDOPE(1)
	MEP1 = IDOPE(2)
	KRANK = IDOPE(3)
	KRP1 = IDOPE(4)
	NSOLN = IDOPE(5)
	NIV = IDOPE(6)
	NIV1 = IDOPE(7)
	L1 = IDOPE(8)
	GO TO 1240
300	CONTINUE
C
C     TO COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
C     SOLVE THE TRIANGULAR SYSTEM OF CURRENTLY NON-ACTIVE
C     VARIABLES AND STORE THE SOLUTION IN Z(*).
C
C     SOLVE-SYSTEM
	ASSIGN 310 TO I99958
	GO TO 1110
C
C     INCREMENT ITERATION COUNTER AND CHECK AGAINST MAX. NUMBER
C     OF ITERATIONS.
310	ITER = ITER + 1
	IF (.NOT.(ITER.GT.ITMAX)) GO TO 320
	MODE = 1
	DONE = .TRUE.
C
C     CHECK TO SEE IF ANY CONSTRAINTS HAVE BECOME ACTIVE.
C     IF SO, CALCULATE AN INTERPOLATION FACTOR SO THAT ALL
C     ACTIVE CONSTRAINTS ARE REMOVED FROM THE BASIS.
320	ALPHA = TWO
	HITCON = .FALSE.
	IF (.NOT.(L.LT.NSOLN)) GO TO 360
	DO 350 J=LP1,NSOLN
	   ZZ = Z(J)
	   IF (.NOT.(ZZ.LE.ZERO)) GO TO 340
	   T = X(J)/(X(J)-ZZ)
	   IF (.NOT.(T.LT.ALPHA)) GO TO 330
	   ALPHA = T
	   JCON = J
330	   HITCON = .TRUE.
340	   CONTINUE
350	CONTINUE
360	GO TO 1220
370	CONTINUE
C
C     TO ADD-CONSTRAINTS
C
C     USE COMPUTED ALPHA TO INTERPOLATE BETWEEN LAST
C     FEASIBLE SOLUTION X(*) AND CURRENT UNCONSTRAINED
C     (AND INFEASIBLE) SOLUTION Z(*).
	IF (.NOT.(LP1.LE.NSOLN)) GO TO 390
	DO 380 J=LP1,NSOLN
	   X(J) = X(J) + ALPHA*(Z(J)-X(J))
380	CONTINUE
390	FEASBL = .FALSE.
	GO TO 410
400	IF (FEASBL) GO TO 610
C
C     REMOVE COL JCON AND SHIFT COLS JCON+1 THROUGH N TO THE
C     LEFT. SWAP COL JCON INTO THE N-TH POSITION.  THIS ACHIEVES
C     UPPER HESSENBERG FORM FOR THE NONACTIVE CONSTRAINTS AND
C     LEAVES AN UPPER HESSENBERG MATRIX TO RETRIANGULARIZE.
410	DO 420 I=1,M
	   T = W(I,JCON)
	   CALL SCOPY(N-JCON, W(I,JCON+1), MDW, W(I,JCON), MDW)
	   W(I,N) = T
420	CONTINUE
C
C     UPDATE PERMUTED INDEX VECTOR TO REFLECT THIS SHIFT AND SWAP.
	ITEMP = IPIVOT(JCON)
	IF (.NOT.(JCON.LT.N)) GO TO 440
	DO 430 I=JCON,NM1
	   IPIVOT(I) = IPIVOT(I+1)
430	CONTINUE
440	IPIVOT(N) = ITEMP
C
C     SIMILARLY REPERMUTE X(*) VECTOR.
	CALL SCOPY(N-JCON, X(JCON+1), 1, X(JCON), 1)
	X(N) = ZERO
	NSP1 = NSOLN
	NSOLN = NSOLN - 1
	NIV1 = NIV
	NIV = NIV - 1
C
C     RETRIANGULARIZE UPPER HESSENBERG MATRIX AFTER ADDING CONSTRAINTS.
	J = JCON
	I = KRANK + JCON - L
450	IF (.NOT.(J.LE.NSOLN)) GO TO 570
	IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0)) GO TO 470
	ASSIGN 460 TO I99938
	GO TO 620
C
C     (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0) ZERO-IP1-TO-I-IN-COL-J
460	GO TO 560
470	IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1)) GO TO 490
	ASSIGN 480 TO I99938
	GO TO 620
C
C     (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1) ZERO-IP1-TO-I-IN-COL-J
480	GO TO 560
490	IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.0)) GO TO 510
	CALL SSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW)
	CALL SSWAP(1, SCALE(I), 1, SCALE(I+1), 1)
	ITEMP = ITYPE(I+1)
	ITYPE(I+1) = ITYPE(I)
	ITYPE(I) = ITEMP
C
C     SWAPPED ROW WAS FORMERLY A PIVOT ELT., SO IT WILL
C     BE LARGE ENOUGH TO PERFORM ELIM.
	ASSIGN 500 TO I99938
	GO TO 620
C
C     ZERO-IP1-TO-I-IN-COL-J
500	GO TO 560
510	IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.1)) GO TO 550
	T = SCALE(I)*W(I,J)**2/ALSQ
	IF (.NOT.(T*FAC+EPENSQ.NE.EPENSQ)) GO TO 530
	ASSIGN 520 TO I99938
	GO TO 620
520	GO TO 540
C
C     WHEN (T*FAC+EPENSQ.NE.EPENSQ) ZERO-IP1-TO-I-IN-COL-J
530	CALL SSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW)
	CALL SSWAP(1, SCALE(I), 1, SCALE(I+1), 1)
	ITEMP = ITYPE(I+1)
	ITYPE(I+1) = ITYPE(I)
	ITYPE(I) = ITEMP
	W(I+1,J) = ZERO
540	CONTINUE
550	CONTINUE
560	I = I + 1
	J = J + 1
	GO TO 450
C
C     SEE IF THE REMAINING COEFFS IN THE SOLN SET ARE FEASIBLE.  THEY
C     SHOULD BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.  IF ANY ARE
C     INFEASIBLE IT IS DUE TO ROUNDOFF ERROR.  ANY THAT ARE NON-
C     POSITIVE WILL BE SET TO ZERO AND REMOVED FROM THE SOLN SET.
570	IF (.NOT.(LP1.LE.NSOLN)) GO TO 590
	DO 580 JCON=LP1,NSOLN
	   IF (X(JCON).LE.ZERO) GO TO 600
580	CONTINUE
590	FEASBL = .TRUE.
600	CONTINUE
	GO TO 400
610	GO TO 1200
620	CONTINUE
C
C     TO ZERO-IP1-TO-I-IN-COL-J
	IF (.NOT.(W(I+1,J).NE.ZERO)) GO TO 630
	CALL SROTMG(SCALE(I), SCALE(I+1), W(I,J), W(I+1,J), SPARAM)
	W(I+1,J) = ZERO
	CALL SROTM(NP1-J, W(I,J+1), MDW, W(I+1,J+1), MDW, SPARAM)
630	GO TO 1290
640	CONTINUE
C
C     TO PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
	CALL SCOPY(NSOLN, Z, 1, X, 1)
	IF (.NOT.(NSOLN.LT.N)) GO TO 650
	X(NSP1) = ZERO
	CALL SCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1)
650	I = NIV1
660	IF (.NOT.(I.LE.ME)) GO TO 690
C
C     RECLASSIFY LEAST SQUARES EQATIONS AS EQUALITIES AS
C     NECESSARY.
	IF (.NOT.(ITYPE(I).EQ.0)) GO TO 670
	I = I + 1
	GO TO 680
670	CALL SSWAP(NP1, W(I,1), MDW, W(ME,1), MDW)
	CALL SSWAP(1, SCALE(I), 1, SCALE(ME), 1)
	ITEMP = ITYPE(I)
	ITYPE(I) = ITYPE(ME)
	ITYPE(ME) = ITEMP
	MEP1 = ME
	ME = ME - 1
680	GO TO 660
C
C     FORM INNER PRODUCT VECTOR WD(*) OF DUAL COEFFS.
690	IF (.NOT.(NSP1.LE.N)) GO TO 730
	DO 720 J=NSP1,N
	   SM = ZERO
	   IF (.NOT.(NSOLN.LT.M)) GO TO 710
	   DO 700 I=NSP1,M
	      SM = SM + SCALE(I)*W(I,J)*W(I,NP1)
700	   CONTINUE
710	   WD(J) = SM
720	CONTINUE
730	GO TO 750
740	IF (POS .OR. DONE) GO TO 970
C
C     FIND J SUCH THAT WD(J)=WMAX IS MAXIMUM.  THIS DETERMINES
C     THAT THE INCOMING COL J WILL REDUCE THE RESIDUAL VECTOR
C     AND BE POSITIVE.
750	WMAX = ZERO
	IWMAX = NSP1
	IF (.NOT.(NSP1.LE.N)) GO TO 780
	DO 770 J=NSP1,N
	   IF (.NOT.(WD(J).GT.WMAX)) GO TO 760
	   WMAX = WD(J)
	   IWMAX = J
760	   CONTINUE
770	CONTINUE
780	IF (.NOT.(WMAX.LE.ZERO)) GO TO 790
	DONE = .TRUE.
	GO TO 960
C
C     SET DUAL COEFF TO ZERO FOR INCOMING COL.
790	WD(IWMAX) = ZERO
C
C     WMAX .GT. ZERO, SO OKAY TO MOVE COL IWMAX TO SOLN SET.
C     PERFORM TRANSFORMATION TO RETRIANGULARIZE, AND TEST
C     FOR NEAR LINEAR DEPENDENCE.
C     SWAP COL IWMAX INTO NSOLN-TH POSITION TO MAINTAIN UPPER
C     HESSENBERG FORM OF ADJACENT COLS, AND ADD NEW COL TO
C     TRIANGULAR DECOMPOSITION.
	NSOLN = NSP1
	NSP1 = NSOLN + 1
	NIV = NIV1
	NIV1 = NIV + 1
	IF (.NOT.(NSOLN.NE.IWMAX)) GO TO 800
	CALL SSWAP(M, W(1,NSOLN), 1, W(1,IWMAX), 1)
	WD(IWMAX) = WD(NSOLN)
	WD(NSOLN) = ZERO
	ITEMP = IPIVOT(NSOLN)
	IPIVOT(NSOLN) = IPIVOT(IWMAX)
	IPIVOT(IWMAX) = ITEMP
C
C     REDUCE COL NOLN SO THAT THE MATRIX OF NONACTIVE
C     CONSTRAINTS VARIABLES IS TRIANGULAR.
800	J = M
810	IF (.NOT.(J.GT.NIV)) GO TO 870
	JM1 = J - 1
	JP = JM1
C
C     WHEN OPERATING NEAR THE ME LINE, TEST TO SEE IF THE PIVOT ELT.
C     IS NEAR ZERO.  IF SO, USE THE LARGEST ELT. ABOVE IT AS THE PIVOT.
C     THIS IS TO MAINTAIN THE SHARP INTERFACE BETWEEN WEIGHTED AND
C     NON-WEIGHTED ROWS IN ALL CASES.
	IF (.NOT.(J.EQ.MEP1)) GO TO 850
	IMAX = ME
	AMAX = SCALE(ME)*W(ME,NSOLN)**2
820	IF (.NOT.(JP.GE.NIV)) GO TO 840
	T = SCALE(JP)*W(JP,NSOLN)**2
	IF (.NOT.(T.GT.AMAX)) GO TO 830
	IMAX = JP
	AMAX = T
830	JP = JP - 1
	GO TO 820
840	JP = IMAX
850	IF (.NOT.(W(J,NSOLN).NE.ZERO)) GO TO 860
	CALL SROTMG(SCALE(JP), SCALE(J), W(JP,NSOLN), W(J,NSOLN), SPARAM)
	W(J,NSOLN) = ZERO
	CALL SROTM(NP1-NSOLN, W(JP,NSP1), MDW, W(J,NSP1), MDW, SPARAM)
860	J = JM1
	GO TO 810
C
C     SOLVE FOR Z(NSOLN)=PROPOSED NEW VALUE FOR X(NSOLN).
C     TEST IF THIS IS NONPOSITIVE OR TOO LARGE.
C     IF THIS WAS TRUE OR IF THE PIVOT TERM WAS ZERO REJECT
C     THE COL AS DEPENDENT.
870	IF (.NOT.(W(NIV,NSOLN).NE.ZERO)) GO TO 890
	ISOL = NIV
	ASSIGN 880 TO I99897
	GO TO 980
C
C     TEST-PROPOSED-NEW-COMPONENT
880	GO TO 940
890	IF (.NOT.(NIV.LE.ME .AND. W(MEP1,NSOLN).NE.ZERO)) GO TO 920
C
C     TRY TO ADD ROW MEP1 AS AN ADDITIONAL EQUALITY CONSTRAINT.
C     CHECK SIZE OF PROPOSED NEW SOLN COMPONENT.
C     REJECT IT IF IT IS TOO LARGE.
	ISOL = MEP1
	ASSIGN 900 TO I99897
	GO TO 980
C
C     TEST-PROPOSED-NEW-COMPONENT
900	IF (.NOT.(POS)) GO TO 910
C
C     SWAP ROWS MEP1 AND NIV, AND SCALE FACTORS FOR THESE ROWS.
	CALL SSWAP(NP1, W(MEP1,1), MDW, W(NIV,1), MDW)
	CALL SSWAP(1, SCALE(MEP1), 1, SCALE(NIV), 1)
	ITEMP = ITYPE(MEP1)
	ITYPE(MEP1) = ITYPE(NIV)
	ITYPE(NIV) = ITEMP
	ME = MEP1
	MEP1 = ME + 1
910	GO TO 930
920	POS = .FALSE.
930	CONTINUE
940	IF (POS) GO TO 950
	NSP1 = NSOLN
	NSOLN = NSOLN - 1
	NIV1 = NIV
	NIV = NIV - 1
950	CONTINUE
960	GO TO 740
970	GO TO 1250
980	CONTINUE
C
C     TO TEST-PROPOSED-NEW-COMPONENT
	Z2 = W(ISOL,NP1)/W(ISOL,NSOLN)
	Z(NSOLN) = Z2
	POS = Z2.GT.ZERO
	IF (.NOT.(Z2.GE.RATNSQ .AND. POS)) GO TO 990
	Z2 = Z2**2
	IF (Z2+RATNSQ.EQ.Z2) POS = .FALSE.
990	GO TO 1280
1000	CONTINUE
C     TO COMPUTE-FINAL-SOLUTION
C
C     SOLVE SYSTEM, STORE RESULTS IN X(*).
C
	ASSIGN 1010 TO I99958
	GO TO 1110
C     SOLVE-SYSTEM
1010	CALL SCOPY(NSOLN, Z, 1, X, 1)
C
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO X(*) IF KRANK.LT.L
	IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.L)) GO TO 1030
	DO 1020 I=1,KRANK
	   CALL H12(2, I, KRP1, L, W(I,1), MDW, H(I), X, 1, 1, 1)
1020	CONTINUE
C
C     FILL IN TRAILING ZEROES FOR CONSTRAINED VARIABLES NOT IN SOLN.
1030	IF (.NOT.(NSOLN.LT.N)) GO TO 1040
	X(NSP1) = ZERO
	CALL SCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1)
C
C     REPERMUTE SOLN VECTOR TO NATURAL ORDER.
1040	DO 1070 I=1,N
	   J = I
1050	   IF (IPIVOT(J).EQ.I) GO TO 1060
	   J = J + 1
	   GO TO 1050
1060	   IPIVOT(J) = IPIVOT(I)
	   IPIVOT(I) = J
	   CALL SSWAP(1, X(J), 1, X(I), 1)
1070	CONTINUE
C
C     RESCALE THE SOLN USING THE COL SCALING.
	DO 1080 J=1,N
	   X(J) = X(J)*D(J)
1080	CONTINUE
	IF (.NOT.(NIV.LT.M)) GO TO 1100
	DO 1090 I=NIV1,M
	   T = SCALE(I)*W(I,NP1)**2
	   IF (I.LE.ME) T = T/ALSQ
	   RNORM = RNORM + T
1090	CONTINUE
1100	RNORM = SQRT(RNORM)
	GO TO 1210
1110	CONTINUE
C
C     TO SOLVE-SYSTEM
C
	IF (.NOT.(DONE)) GO TO 1120
	ISOL = 1
	GO TO 1130
1120	ISOL = LP1
1130	IF (.NOT.(NSOLN.GE.ISOL)) GO TO 1190
C
C     COPY RT. HAND SIDE INTO TEMP VECTOR TO USE OVERWRITING METHOD.
	CALL SCOPY(NIV, W(1,NP1), 1, TEMP, 1)
	DO 1180 JJ=ISOL,NSOLN
	   J = NSOLN - JJ + ISOL
	   IF (.NOT.(J.GT.KRANK)) GO TO 1140
	   I = NIV - JJ + ISOL
	   GO TO 1150
1140	   I = J
1150	   IF (.NOT.(J.GT.KRANK .AND. J.LE.L)) GO TO 1160
	   Z(J) = ZERO
	   GO TO 1170
1160	   Z(J) = TEMP(I)/W(I,J)
	   CALL SAXPY(I-1, -Z(J), W(1,J), 1, TEMP, 1)
1170	   CONTINUE
1180	CONTINUE
1190	GO TO 1270
1200	GO TO I99986, (40)
1210	GO TO I99980, (90)
1220	GO TO I99991, (30)
1230	GO TO I99998, (10)
1240	GO TO I99995, (20)
1250	GO TO I99983, (60)
1260	GO TO I99977, (220)
1270	GO TO I99958, (310, 1010)
1280	GO TO I99897, (880, 900)
1290	GO TO I99938, (460, 480, 500, 520)
	END

	INTEGER*4 FUNCTION ISAMAX(N,SX,INCX)
	IMPLICIT INTEGER*4 (I-N)
C
C     FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF SINGLE PRECISION SX.
C     ISAMAX =  FIRST I, I = 1 TO N, TO MINIMIZE  ABS(SX(1-INCX+I*INCX))
C
	REAL SX(1),SMAX,XMAG
	ISAMAX = 0
	IF(N.LE.0) RETURN
	ISAMAX = 1
	IF(N.LE.1)RETURN
	IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
	SMAX = ABS(SX(1))
	NS = N*INCX
	II = 1
	    DO 10 I=1,NS,INCX
	    XMAG = ABS(SX(I))
	    IF(XMAG.LE.SMAX) GO TO 5
	    ISAMAX = II
	    SMAX = XMAG
5	    II = II + 1
10	    CONTINUE
	RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
20	SMAX = ABS(SX(1))
	DO 30 I = 2,N
	   XMAG = ABS(SX(I))
	   IF(XMAG.LE.SMAX) GO TO 30
	   ISAMAX = I
	   SMAX = XMAG
30	CONTINUE
	RETURN
	END

	SUBROUTINE SROTM (N,SX,INCX,SY,INCY,SPARAM)
	IMPLICIT INTEGER*4 (I-N)
C
C     APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
C
C     (SX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF SX ARE IN
C     (DX**T)
C
C     SX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
C     LX = (-INCX)*N, AND SIMILARLY FOR SY USING USING LY AND INCY.
C     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
C
C       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
C     H=(          )    (          )    (          )    (          )
C       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
C     SEE  SROTMG FOR A DESCRIPTION OF DATA STORAGE IN SPARAM.
C
	DIMENSION SX(1),SY(1),SPARAM(5)
	DATA ZERO,TWO/0.E0,2.E0/
C
	SFLAG=SPARAM(1)
	IF(N .LE. 0 .OR.(SFLAG+TWO.EQ.ZERO)) GO TO 140
	    IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70
C
	         NSTEPS=N*INCX
	         IF(SFLAG) 50,10,30
10	         CONTINUE
	         SH12=SPARAM(4)
	         SH21=SPARAM(3)
	              DO 20 I=1,NSTEPS,INCX
	              W=SX(I)
	              Z=SY(I)
	              SX(I)=W+Z*SH12
	              SY(I)=W*SH21+Z
20	              CONTINUE
	         GO TO 140
30	         CONTINUE
	         SH11=SPARAM(2)
	         SH22=SPARAM(5)
	              DO 40 I=1,NSTEPS,INCX
	              W=SX(I)
	              Z=SY(I)
	              SX(I)=W*SH11+Z
	              SY(I)=-W+SH22*Z
40	              CONTINUE
	         GO TO 140
50	         CONTINUE
	         SH11=SPARAM(2)
	         SH12=SPARAM(4)
	         SH21=SPARAM(3)
	         SH22=SPARAM(5)
	              DO 60 I=1,NSTEPS,INCX
	              W=SX(I)
	              Z=SY(I)
	              SX(I)=W*SH11+Z*SH12
	              SY(I)=W*SH21+Z*SH22
60	              CONTINUE
	         GO TO 140
70	    CONTINUE
	    KX=1
	    KY=1
	    IF(INCX .LT. 0) KX=1+(1-N)*INCX
	    IF(INCY .LT. 0) KY=1+(1-N)*INCY
C
	    IF(SFLAG)120,80,100
80	    CONTINUE
	    SH12=SPARAM(4)
	    SH21=SPARAM(3)
	         DO 90 I=1,N
	         W=SX(KX)
	         Z=SY(KY)
	         SX(KX)=W+Z*SH12
	         SY(KY)=W*SH21+Z
	         KX=KX+INCX
	         KY=KY+INCY
90	         CONTINUE
	    GO TO 140
100	    CONTINUE
	    SH11=SPARAM(2)
	    SH22=SPARAM(5)
	         DO 110 I=1,N
	         W=SX(KX)
	         Z=SY(KY)
	         SX(KX)=W*SH11+Z
	         SY(KY)=-W+SH22*Z
	         KX=KX+INCX
	         KY=KY+INCY
110	         CONTINUE
	    GO TO 140
120	    CONTINUE
	    SH11=SPARAM(2)
	    SH12=SPARAM(4)
	    SH21=SPARAM(3)
	    SH22=SPARAM(5)
	         DO 130 I=1,N
	         W=SX(KX)
	         Z=SY(KY)
	         SX(KX)=W*SH11+Z*SH12
	         SY(KY)=W*SH21+Z*SH22
	         KX=KX+INCX
	         KY=KY+INCY
130	         CONTINUE
140	    CONTINUE
	    RETURN
	    END

	SUBROUTINE WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM,
     * IDOPE, DOPE, DONE)
	IMPLICIT INTEGER*4 (I-N)
C     THIS IS A COMPANION SUBPROGRAM TO WNNLS( ).
C     THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE
C     USAGE INSTRUCTIONS.
C
C     NOTE  THE M BY (N+1) MATRIX W( , ) CONTAINS THE RT. HAND SIDE
C           B AS THE (N+1)ST COL.
C
C
C     TRIANGULARIZE L1 BY L1 SUBSYSTEM, WHERE L1=MIN(M,L), WITH
C     COL INTERCHANGES.
C     REVISED MAY 17, 1979.
C
	DIMENSION W(MDW,1), IPIVOT(1), ITYPE(1), H(1), SCALE(1), SPARAM(5)
	DIMENSION IDOPE(8), DOPE(4)
	LOGICAL INDEP, DONE, RECALC
	DATA TENM3 /1.E-3/, ZERO /0.E0/, ONE /1.E0/
C
	ME = IDOPE(1)
	MEP1 = IDOPE(2)
	KRANK = IDOPE(3)
	KRP1 = IDOPE(4)
	NSOLN = IDOPE(5)
	NIV = IDOPE(6)
	NIV1 = IDOPE(7)
	L1 = IDOPE(8)
C
	ALSQ = DOPE(1)
	EPENSQ = DOPE(2)
	FAC = DOPE(3)
	TAU = DOPE(4)
	NP1 = N + 1
	LB = MIN0(M-1,L)
	RECALC = .TRUE.
	RNORM = ZERO
	KRANK = 0
C     WE SET FACTOR=1.E0 SO THAT THE HEAVY WEIGHT ALAMDA WILL BE
C     INCLUDED IN THE TEST FOR COL INDEPENDENCE.
	FACTOR = 1.E0
	I = 1
	IP1 = 2
10	IF (.NOT.(I.LE.LB)) GO TO 130
C
C     SET IR TO POINT TO THE I-TH ROW.
	IR = I
	LEND = L
	MEND = M
	ASSIGN 20 TO I99996
	GO TO 450
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
20	ASSIGN 30 TO I99993
	GO TO 550
C
C     PERFORM-COL-INTERCHANGE
C
C     SET IC TO POINT TO I-TH COL.
30	IC = I
	ASSIGN 40 TO I99990
	GO TO 510
C
C     TEST-INDEP-OF-INCOMING-COL
40	IF (.NOT.(INDEP)) GO TO 110
C
C     ELIMINATE I-TH COL BELOW DIAG. USING MOD. GIVENS TRANSFORMATIONS
C     APPLIED TO (A B).
	J = M
	DO 100 JJ=IP1,M
	   JM1 = J - 1
	   JP = JM1
C     WHEN OPERATING NEAR THE ME LINE, USE THE LARGEST ELT.
C     ABOVE IT AS THE PIVOT.
	   IF (.NOT.(J.EQ.MEP1)) GO TO 80
	   IMAX = ME
	   AMAX = SCALE(ME)*W(ME,I)**2
50	   IF (.NOT.(JP.GE.I)) GO TO 70
	   T = SCALE(JP)*W(JP,I)**2
	   IF (.NOT.(T.GT.AMAX)) GO TO 60
	   IMAX = JP
	   AMAX = T
60	   JP = JP - 1
	   GO TO 50
70	   JP = IMAX
80	   IF (.NOT.(W(J,I).NE.ZERO)) GO TO 90
	   CALL SROTMG(SCALE(JP), SCALE(J), W(JP,I), W(J,I), SPARAM)
	   W(J,I) = ZERO
	   CALL SROTM(NP1-I, W(JP,IP1), MDW, W(J,IP1), MDW, SPARAM)
90	   J = JM1
100	CONTINUE
	KRANK = KRANK + 1
	GO TO 120
110	GO TO 140
120	I = IP1
	IP1 = IP1 + 1
	GO TO 10
130	KRANK = L1
140	CONTINUE
	KRP1 = KRANK + 1
	IF (.NOT.(KRANK.LT.ME)) GO TO 260
	FACTOR = ALSQ
C
C     DETERMINE THE RANK OF THE REMAINING EQUALITY CONSTRAINT
C     EQUATIONS.  REMOVE ANY REDUNDANT EQUATIONS.
	LP1 = L + 1
	RECALC = .TRUE.
	LB = MIN0(L+ME-KRANK,N)
	I = LP1
	IP1 = I + 1
150	IF (.NOT.(I.LE.LB)) GO TO 250
	IR = KRANK + I - L
	LEND = N
	MEND = ME
	ASSIGN 160 TO I99996
	GO TO 450
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
160	ASSIGN 170 TO I99993
	GO TO 550
C
C     PERFORM-COL-INTERCHANGE
C
C     ELIMINATE ELEMENTS IN THE I-TH COL.
170	J = ME
180	IF (.NOT.(J.GT.IR)) GO TO 200
	JM1 = J - 1
	IF (.NOT.(W(J,I).NE.ZERO)) GO TO 190
	CALL SROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM)
	W(J,I) = ZERO
	CALL SROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM)
190	J = JM1
	GO TO 180
C
C     SET IC=I=COL BEING ELIMINATED
200	IC = I
	ASSIGN 210 TO I99990
	GO TO 510
C
C     TEST-INDEP-OF-INCOMING-COL
210	IF (INDEP) GO TO 240
C
C     REMOVE ANY REDUNDANT OR DEPENDENT EQUALITY CONSTRAINTS.
	JJ = IR
220	IF (.NOT.(IR.LE.ME)) GO TO 230
	W(IR,1) = ZERO
	CALL SCOPY(N, W(IR,1), 0, W(IR,1), MDW)
	RNORM = RNORM + SCALE(IR)*W(IR,NP1)**2/ALSQ
	W(IR,NP1) = ZERO
	SCALE(IR) = ONE
	IR = IR + 1
	GO TO 220
C
C     REDUCE ME TO REFLECT ANY DISCOVERED DEPENDENT EQUALITY
C     CONSTRAINTS.
230	ME = JJ - 1
	MEP1 = ME + 1
	GO TO 270
240	I = IP1
	IP1 = IP1 + 1
	GO TO 150
250	CONTINUE
260	CONTINUE
270	CONTINUE
	IF (.NOT.(KRANK.LT.L1)) GO TO 410
C
C     TRY TO DETERMINE THE VARIABLES KRANK+1 THROUGH L1 FROM THE
C     LEAST SQUARES EQUATIONS.  CONTINUE THE TRIANGULARIZATION WITH
C     PIVOT ELEMENT W(MEP1,I).
C
	RECALC = .TRUE.
C
C     SET FACTOR=ALSQ TO REMOVE EFFECT OF HEAVY WEIGHT FROM
C     TEST FOR COL INDEPENDENCE.
	FACTOR = ALSQ
	KK = KRP1
	I = KK
	IP1 = I + 1
280	IF (.NOT.(I.LE.L1)) GO TO 400
C
C     SET IR TO POINT TO THE MEP1-ST ROW.
	IR = MEP1
	LEND = L
	MEND = M
	ASSIGN 290 TO I99996
	GO TO 450
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
290	ASSIGN 300 TO I99993
	GO TO 550
C
C     PERFORM-COL-INTERCHANGE
C
C     ELIMINATE I-TH COL BELOW THE IR-TH ELEMENT.
300	IRP1 = IR + 1
	J = M
	DO 320 JJ=IRP1,M
	   JM1 = J - 1
	   IF (.NOT.(W(J,I).NE.ZERO)) GO TO 310
	   CALL SROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM)
	   W(J,I) = ZERO
	   CALL SROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM)
310	   J = JM1
320	CONTINUE
C
C     TEST IF NEW PIVOT ELEMENT IS NONZERO. IF NOT, THE COL IS
C     DEPENDENT.
	T = SCALE(IR)*W(IR,I)**2
	IF (.NOT.(FAC*T+EPENSQ.EQ.EPENSQ)) GO TO 330
	INDEP = .FALSE.
	GO TO 360
C
C     COL TEST PASSED. NOW MUST PASS ROW NORM TEST TO BE CLASSIFIED
C     AS INDEPENDENT.
330	RN = ZERO
	DO 350 I1=IR,M
	   DO 340 J1=IP1,N
	      RN = AMAX1(RN,SCALE(I1)*W(I1,J1)**2)
340	   CONTINUE
350	CONTINUE
	INDEP = FAC*T + RN.NE.RN
C
C     IF INDEPENDENT, SWAP THE IR-TH AND KRP1-ST ROWS TO MAINTAIN THE
C     TRIANGULAR FORM.  UPDATE THE RANK INDICATOR KRANK AND THE
C     EQUALITY CONSTRAINT POINTER ME.
360	IF (.NOT.(INDEP)) GO TO 380
	IF (.NOT.(IR.NE.KRP1)) GO TO 370
	CALL SSWAP(NP1, W(KRP1,1), MDW, W(IR,1), MDW)
	CALL SSWAP(1, SCALE(KRP1), 1, SCALE(IR), 1)
	ITEMP = ITYPE(KRP1)
	ITYPE(KRP1) = ITYPE(IR)
	ITYPE(IR) = ITEMP
370	ME = MEP1
	MEP1 = ME + 1
	KRANK = KRP1
	KRP1 = KRANK + 1
	GO TO 390
380	GO TO 420
390	I = IP1
	IP1 = IP1 + 1
	GO TO 280
400	CONTINUE
410	CONTINUE
420	CONTINUE
C
C     IF PSEUDORANK IS LESS THAN L, APPLY HOUSEHOLDER TRANS.
C     FROM RIGHT.
	IF (.NOT.(KRANK.LT.L)) GO TO 440
	DO 430 I=1,KRANK
	   J = KRP1 - I
	   CALL H12(1, J, KRP1, L, W(J,1), MDW, H(J), W, MDW, 1, J-1)
430	CONTINUE
440	NIV = KRANK + NSOLN - L
	NIV1 = NIV + 1
	IF (L.EQ.N) DONE = .TRUE.
C
C  END OF INITIAL TRIANGULARIZATION.
	IDOPE(1) = ME
	IDOPE(2) = MEP1
	IDOPE(3) = KRANK
	IDOPE(4) = KRP1
	IDOPE(5) = NSOLN
	IDOPE(6) = NIV
	IDOPE(7) = NIV1
	IDOPE(8) = L1
	RETURN
450	CONTINUE
C
C     TO UPDATE-COL-SS-AND-FIND-PIVOT-COL
C
C     THE COL SS VECTOR WILL BE UPDATED AT EACH STEP. WHEN
C     NUMERICALLY NECESSARY, THESE VALUES WILL BE RECOMPUTED.
C
	IF (.NOT.(IR.NE.1 .AND. (.NOT.RECALC))) GO TO 470
C     UPDATE COL SS =SUM OF SQUARES.
	DO 460 J=I,LEND
	   H(J) = H(J) - SCALE(IR-1)*W(IR-1,J)**2
460	CONTINUE
C
C     TEST FOR NUMERICAL ACCURACY.
	MAX = ISAMAX(LEND-I+1,H(I),1) + I - 1
	RECALC = HBAR + TENM3*H(MAX).EQ.HBAR
C
C     IF REQUIRED, RECALCULATE COL SS, USING ROWS IR THROUGH MEND.
470	IF (.NOT.(RECALC)) GO TO 500
	DO 490 J=I,LEND
	   H(J) = ZERO
	   DO 480 K=IR,MEND
	      H(J) = H(J) + SCALE(K)*W(K,J)**2
480	   CONTINUE
490	CONTINUE
C
C     FIND COL WITH LARGEST SS.
	MAX = ISAMAX(LEND-I+1,H(I),1) + I - 1
	HBAR = H(MAX)
500	GO TO 590
510	CONTINUE
C
C     TO TEST-INDEP-OF-INCOMING-COL
C
C     TEST THE COL IC TO DETERMINE IF IT IS LINEARLY INDEPENDENT
C     OF THE COLS ALREADY IN THE BASIS.  IN THE INIT TRI
C     STEP, WE USUALLY WANT THE HEAVY WEIGHT ALAMDA TO
C     BE INCLUDED IN THE TEST FOR INDEPENDENCE.  IN THIS CASE THE
C     VALUE OF FACTOR WILL HAVE BEEN SET TO 1.E0 BEFORE THIS
C     PROCEDURE IS INVOKED.  IN THE POTENTIALLY RANK DEFICIENT
C     PROBLEM, THE VALUE OF FACTOR WILL HAVE BEEN
C     SET TO ALSQ=ALAMDA**2 TO REMOVE THE EFFECT OF THE HEAVY WEIGHT
C     FROM THE TEST FOR INDEPENDENCE.
C
C     WRITE NEW COL AS PARTITIONED VECTOR
C             (A1)  NUMBER OF COMPONENTS IN SOLN SO FAR = NIV
C             (A2)  M-NIV COMPONENTS
C     AND COMPUTE  SN = INVERSE WEIGHTED LENGTH OF A1
C                  RN = INVERSE WEIGHTED LENGTH OF A2
C     CALL THE COL INDEPENDENT WHEN RN .GT. TAU*SN
	SN = ZERO
	RN = ZERO
	DO 540 J=1,MEND
	   T = SCALE(J)
	   IF (J.LE.ME) T = T/FACTOR
	   T = T*W(J,IC)**2
	   IF (.NOT.(J.LT.IR)) GO TO 520
	   SN = SN + T
	   GO TO 530
520	   RN = RN + T
530	   CONTINUE
540	CONTINUE
	INDEP = RN.GT.TAU**2*SN
	GO TO 580
550	CONTINUE
C
C     TO PERFORM-COL-INTERCHANGE
C
	IF (.NOT.(MAX.NE.I)) GO TO 560
C     EXCHANGE ELEMENTS OF PERMUTED INDEX VECTOR AND PERFORM COL
C     INTERCHANGES.
	ITEMP = IPIVOT(I)
	IPIVOT(I) = IPIVOT(MAX)
	IPIVOT(MAX) = ITEMP
	CALL SSWAP(M, W(1,MAX), 1, W(1,I), 1)
	H(MAX) = H(I)
560	GO TO 570
570	GO TO I99993, (30, 170, 300)
580	GO TO I99990, (40, 210)
590	GO TO I99996, (20, 160, 290)
	END

	SUBROUTINE SROTMG (SD1,SD2,SX1,SY1,SPARAM)
	IMPLICIT INTEGER*4 (I-N)
C
C     CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
C     THE SECOND COMPONENT OF THE 2-VECTOR  (SQRT(SD1)*SX1,SQRT(SD2)*
C     SY2)**T.
C     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
C
C       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
C     H=(          )    (          )    (          )    (          )
C       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
C     LOCATIONS 2-4 OF SPARAM CONTAIN SH11,SH21,SH12, AND SH22
C     RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE
C     VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.)
C
C     THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
C     INEXACT.  THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
C     OF SD1 AND SD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
C
	DIMENSION SPARAM(5)
C
	DATA ZERO,ONE,TWO /0.E0,1.E0,2.E0/
	DATA GAM,GAMSQ,RGAMSQ/4096.E0,1.67772E7,5.96046E-8/
	IF(.NOT. SD1 .LT. ZERO) GO TO 10
C       GO ZERO-H-D-AND-SX1..
	    GO TO 60
10	CONTINUE
C     CASE-SD1-NONNEGATIVE
	SP2=SD2*SY1
	IF(.NOT. SP2 .EQ. ZERO) GO TO 20
	    SFLAG=-TWO
	    GO TO 260
C     REGULAR-CASE..
20	CONTINUE
	SP1=SD1*SX1
	SQ2=SP2*SY1
	SQ1=SP1*SX1
C
	IF(.NOT. ABS(SQ1) .GT. ABS(SQ2)) GO TO 40
	    SH21=-SY1/SX1
	    SH12=SP2/SP1
C
	    SU=ONE-SH12*SH21
C
	    IF(.NOT. SU .LE. ZERO) GO TO 30
C         GO ZERO-H-D-AND-SX1..
	         GO TO 60
30	    CONTINUE
	         SFLAG=ZERO
	         SD1=SD1/SU
	         SD2=SD2/SU
	         SX1=SX1*SU
C         GO SCALE-CHECK..
	         GO TO 100
40	CONTINUE
	    IF(.NOT. SQ2 .LT. ZERO) GO TO 50
C         GO ZERO-H-D-AND-SX1..
	         GO TO 60
50	    CONTINUE
	         SFLAG=ONE
	         SH11=SP1/SP2
	         SH22=SX1/SY1
	         SU=ONE+SH11*SH22
	         STEMP=SD2/SU
	         SD2=SD1/SU
	         SD1=STEMP
	         SX1=SY1*SU
C         GO SCALE-CHECK
	         GO TO 100
C     PROCEDURE..ZERO-H-D-AND-SX1..
60	CONTINUE
	    SFLAG=-ONE
	    SH11=ZERO
	    SH12=ZERO
	    SH21=ZERO
	    SH22=ZERO
C
	    SD1=ZERO
	    SD2=ZERO
	    SX1=ZERO
C         RETURN..
	    GO TO 220
C     PROCEDURE..FIX-H..
70	CONTINUE
	IF(.NOT. SFLAG .GE. ZERO) GO TO 90
C
	    IF(.NOT. SFLAG .EQ. ZERO) GO TO 80
	    SH11=ONE
	    SH22=ONE
	    SFLAG=-ONE
	    GO TO 90
80	    CONTINUE
	    SH21=-ONE
	    SH12=ONE
	    SFLAG=-ONE
90	CONTINUE
	GO TO IGO,(120,150,180,210)
C     PROCEDURE..SCALE-CHECK
100	CONTINUE
110	    CONTINUE
	    IF(.NOT. SD1 .LE. RGAMSQ) GO TO 130
	         IF(SD1 .EQ. ZERO) GO TO 160
	         ASSIGN 120 TO IGO
C              FIX-H..
	         GO TO 70
120	         CONTINUE
	         SD1=SD1*GAM**2
	         SX1=SX1/GAM
	         SH11=SH11/GAM
	         SH12=SH12/GAM
	    GO TO 110
130	CONTINUE
140	    CONTINUE
	    IF(.NOT. SD1 .GE. GAMSQ) GO TO 160
	         ASSIGN 150 TO IGO
C              FIX-H..
	         GO TO 70
150	         CONTINUE
	         SD1=SD1/GAM**2
	         SX1=SX1*GAM
	         SH11=SH11*GAM
	         SH12=SH12*GAM
	    GO TO 140
160	CONTINUE
170	    CONTINUE
	    IF(.NOT. ABS(SD2) .LE. RGAMSQ) GO TO 190
	         IF(SD2 .EQ. ZERO) GO TO 220
	         ASSIGN 180 TO IGO
C              FIX-H..
	         GO TO 70
180	         CONTINUE
	         SD2=SD2*GAM**2
	         SH21=SH21/GAM
	         SH22=SH22/GAM
	    GO TO 170
190	CONTINUE
200	    CONTINUE
	    IF(.NOT. ABS(SD2) .GE. GAMSQ) GO TO 220
	         ASSIGN 210 TO IGO
C              FIX-H..
	         GO TO 70
210	         CONTINUE
	         SD2=SD2/GAM**2
	         SH21=SH21*GAM
	         SH22=SH22*GAM
	    GO TO 200
220	CONTINUE
	    IF(SFLAG)250,230,240
230	    CONTINUE
	         SPARAM(3)=SH21
	         SPARAM(4)=SH12
	         GO TO 260
240	    CONTINUE
	         SPARAM(2)=SH11
	         SPARAM(5)=SH22
	         GO TO 260
250	    CONTINUE
	         SPARAM(2)=SH11
	         SPARAM(3)=SH21
	         SPARAM(4)=SH12
	         SPARAM(5)=SH22
260	CONTINUE
	    SPARAM(1)=SFLAG
	    RETURN
	END
