      SUBROUTINE WNNLS(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,
     * IWORK, WORK)
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 THEN MODE
C                  IS SET TO 2 AND THE SUBROUTINE TERMINATES.
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. EITHER ONE OF
C                            THE RESTRICTIONS MDW.GE.ME+MA AND
C                            0.LE.L.LE.N IS NOT SATISFIED, OR
C                            THE OPTION VECTOR PRGOPT(*) IS NOT
C                            PROPERLY DEFINED.
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.
C
C     SPMPAR        FUNCTION TO COMPUTE THE RELATIVE MACHINE
C                   PRECISION.
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     REVISED SEPTEMBER 6, 1979.
      DIMENSION W(MDW,1), PRGOPT(1), X(1), IWORK(1), WORK(1)
C
      MODE = 0
      M = ME+MA
      IF (M.LE.0 .OR. N.LE.0) RETURN
      IF (MDW.GE.M) GO TO 10
      MODE = 2
      RETURN
   10 IF (0.LE.L .AND. L.LE.N) GO TO 20
      MODE = 2
      RETURN
C
C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS REQUIRED
C     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
