      subroutine wnnls(w,mdw,me,ma,n,l,prgopt,x,rnorm,mode,iwork,work)
c***begin prologue  wnnls
c***date written   790701   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  k1a2a
c***keywords  constrained least squares,curve fitting,data fitting,
c             equality constraints,inequality constraints,
c             nonnegativity constraints,quadratic programming
c***author  hanson, r. j., (snla)
c           haskell, k. h., (snla)
c***purpose  solve a linearly constrained least squares problem with
c            equality constraints and nonnegativity constraints on
c            selected variables.
c***description
c
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(srelpr).  this quantity
c                  can be no smaller than srelpr, the arithmetic-
c                  storage precision.  the quantity used
c                  here is internally restricted to be at
c                  least srelpr.  the data set for this option
c                  is the new tolerance.
c
c                  key=9
c                         change the blow-up parameter from the
c                  nominal value of sqrt(srelpr).  the reciprocal of
c                  this parameter is used in rejecting solution
c                  components as too large when a variable is
c                  first brought into the active set.  too large
c                  means that the proposed component times the
c                  reciprocal of the parameter is not less than
c                  the ratio of the norms of the right-side
c                  vector and the data matrix.
c                  this parameter can be no smaller than srelpr,
c                  the arithmetic-storage precision.
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     iwork(1),    the amounts of working storage actually allocated
c     iwork(2)     for the working arrays work(*) and iwork(*),
c                  respectively.  these quantities are compared with
c                  the actual amounts of storage needed for wnnls( ).
c                  insufficient storage allocated for either work(*)
c                  or iwork(*) is considered an error.  this feature
c                  was included in wnnls( ) because miscalculating
c                  the storage formulas for work(*) and iwork(*)
c                  might very well lead to subtle and hard-to-find
c                  execution errors.
c
c                  the length of work(*) must be at least
c
c                  lw = me+ma+5*n
c                  this test will not be made if iwork(1).le.0.
c
c                  the length of iwork(*) must be at least
c
c                  liw = me+ma+n
c                  this test will not be made if iwork(2).le.0.
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***references  k.h. haskell and r.j. hanson, *an algorithm for
c                 linear least squares problems with equality and
c                 nonnegativity constraints*, sand77-0552, june 1978.
c               k.h. haskell and r.j. hanson, *selected algorithms for
c                 the linearly constrained least squares problem--
c                 a users guide*, sand78-1290, august 1979.
c               k.h. haskell and r.h. hanson, *an algorithm for
c                 linear least squares problems with equality and
c                 nonnegativity constraints*, math. prog. 21 (1981),
c                 pp. 98-118.
c               r.j. hanson and k.h. haskell, *two algorithms for the
c                 linearly constrained least squares problem*, acm
c                 trans. on math. software, sept. 1982.
c***routines called  wnlsm,xerror,xerrwv
c***end prologue  wnnls
c
c     the editing required to convert this subroutine from single to
c     double precision involves the following character string changes.
c     use an editing command (change) /string-1/(to)string-2/.
c     (start at line with c++ in cols. 1-3.)
c     /real (12 blanks)/double precision/,/, dummy/,sngl(dummy)/
c
c     written by karen h. haskell, sandia laboratories,
c     and r.j. hanson, sandia laboratories.
c     revised feb.25, 1982.
c
c     subroutines called by wnnls( )
c
c++
c     wnlsm         companion subroutine to wnnls( ), where
c                   most of the computation takes place.
c
c     xerror,xerrwv 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        toms, v. 5, no. 3, p. 308.  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, and math. programming, vol. 21, (1981), p. 98-118.
c
c     4. slatec common math. library error handling
c        package.  by r. e. jones.  available as sandia
c        technical report sand78-1189.
c
      real              dummy, w(mdw,1), prgopt(1), x(1), work(1), rnorm
      integer iwork(*)
c
c
c***first executable statement  wnnls
      mode = 0
       iwork(1)=mdw*6
       iwork(2)=mdw*2
      if (ma+me.le.0 .or. n.le.0) return
      if (.not.(iwork(1).gt.0)) go to 20
      lw = me + ma + 5*n
      if (.not.(iwork(1).lt.lw)) go to 10
      nerr = 2
      iopt = 1
      print*,'work array',iwork(1),lw
      call xerrwv( 'wnnls( ), insufficient storage allocated for work(*)
     1, need lw=i1 below', 70, nerr, iopt, 1, lw, 0, 0, dummy, dummy)
      mode = 2
      return
   10 continue
   20 if (.not.(iwork(2).gt.0)) go to 40
      liw = me + ma + n
      if (.not.(iwork(2).lt.liw)) go to 30
      nerr = 2
      iopt = 1
      call xerrwv( 'wnnls( ), insufficient storage allocated for iwork(*
     1), need liw=i1 below', 72, nerr, iopt, 1, liw, 0, 0, dummy, dummy)
      mode = 2
      return
   30 continue
   40 if (.not.(mdw.lt.me+ma)) go to 50
      nerr = 1
      iopt = 1
      call xerror( 'wnnls( ), the value mdw.lt.me+ma is an error', 44,
     1 nerr, iopt)
      mode = 2
      return
   50 if (0.le.l .and. l.le.n) go to 60
      nerr = 2
      iopt = 1
      call xerror( 'wnnls( ), l.le.0.and.l.le.n is required', 39, nerr,
     1 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
   60 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,
     1 iwork(l1), work(1), work(l1), work(l2), work(l3), work(l4),
     2 work(l5))
      return
      end
