c......this subroutine is for positivitity constraints
      subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode)
c
c...........nonnegative least squares............
c...........from C.L.Lawson and R.J. Hanson,'solving for least 
c          squares problems'
c        Given an m by n matrix, a, and an m-vector, b, computer
c        an n-vector, x, which solves the least squares problem
c
c               a * x = b    subject to x .ge. 0.0
c     a(),mda,m,n   mda is the first dimensioning parameter for
c                   the array, a().  on entry a() contains the m
c                   by n matrix, a.    on exit a() contains the 
c                   product matrix, q*a, where q is an m by m
c                   orthogonal matrix generated implicitly by
c                   this subroutine.
c     b()   on entry b() contains the m-vector, b.  on exit b()
c            contains q*b.
c     x()   on entry x() need not be initialized.   on exit x()
c            will contain the solution vector.
c     rnorm on exit rnorm contains the euclidean norm of the
c           residual vector.
c     w()   an n-array of working space. on exit w() will contain
c           the dual solution vectors.  w will satisfy w(i)=0.
c           for all i in set p and w(i) .le. 0. for all i in set z
c     zz()  an m-array of working space.
c     index()     an integer working array of length at least n.
c                 on exit the contents of this array define the 
c                 sets p and z as follows:
c                   index(1)  thru index(nsetp) = set p
c                   index(iz1) thru index(iz2)  = set z.
c                   iz1 = nsetp + 1 = npp1
c                   iz2 = n
c     mode   this is a success-failure flag with the following
c            meanings:
c                   1. the solution has been computed successfully.
c                   2. the dimensions of the problem are bad.
c                       either m .le. 0  or n .le. 0.
c                   3. iteration count exceeded. more than 3*n iterations
c
      dimension a(mda,n),b(m),x(n),w(n),zz(m)
      integer index(n)
c
       zero=0.
       one=1.
       two=2.
       factor=0.01
c
       mode=1
       if (m.gt.0.and.n.gt.0) go to 10
       mode=2
       return
  10   iter=0
       itmax=3*n
c
c..... initialize the arrays index() and x().
c
        do 20 i=1,n
           x(i)=zero
  20       index(i)=i
c
       iz2=n
       iz1=1
       nsetp=0
       npp1=1
c            ****** main loop begins here ******
  30   continue 
c      quit if all coefficients are already in the solution.
c      or if m cols of a have been triangularized.
c
       if (iz1.gt.iz2.or.nsetp.ge.m) go to 350
c
c          compute components of the dual (negative gradient)
c            vector w()
       do 50 iz=iz1,iz2
          j=index(iz)
          sm=zero
              do 40 l=npp1,m
  40          sm=sm+a(l,j)*b(l)
  50      w(j)=sm
c                     find largest positive w(j)
  60   wmax=zero
            do 70 iz=iz1,iz2
            j=index(iz)
            if (w(j).le.wmax) go to 70
            wmax=w(j)
            izmax=iz
  70        continue
c
c          if wmax.le.0. go to termination.
c          this indicates satisfaction of the Kuhn-Tucker conditions.
c
       if (wmax) 350,350,80
  80   iz=izmax
       j=index(iz)
c
c     the sign of w(j) is ok for j to be moved to set p.
c     begin the fransformation and check new diagonal element
c     to avoid near linear dependence.
c
       asave=a(npp1,j)
       call h12(1,npp1,npp1+1,m,a(1,j),1,up,dummy,1,1,0)
       unorm=zero
       if (nsetp.eq.0) go to 100
          do 90 l=1,nsetp
  90      unorm=unorm+a(l,j)**2
  100  unorm=sqrt(unorm)
       if (diff(unorm+abs(a(npp1,j))*factor,unorm)) 130,130,110
c
c    col j is sufficiently independent. copy b into zz, update
c    zz and solve for ztest ( =proposed new value for x(j) )
  110  do 120 l=1,m
  120    zz(l)=b(l)
       call h12(2,npp1,npp1+1,m,a(1,j),1,up,zz,1,1,1)
       ztest=zz(npp1)/a(npp1,j)
c
c            see if ztest is positive
c
       if (ztest) 130,130,140
c
c   reject j as a candidate to be moved from set z to p.
c   restore a(npp1,j), set w(j)=0., and back to test 
c   dual coeffs again.
c
  130  a(npp1,j)= asave
         w(j)=zero
         go to 60
c
c   the index j=index(iz) has been selected to be moved from set 
c   z to set p.  update b, update indices, apply householder 
c   transformations to cols in new set z, zero subdiagonal elts 
c   in col j, set w(j)=0.
c
  140  do 150 l=1,m
  150  b(l)=zz(l)
c
       index(iz)=index(iz1)
       index(iz1)=j
       iz1=iz1+1
       nsetp=npp1
       npp1=npp1+1
c
       if (iz1.gt.iz2) go to 170
         do 160 jz=iz1,iz2
         jj=index(jz)
  160    call h12(2,nsetp,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
  170  continue
c
       if (nsetp.eq.m) go to 190
          do 180 l=npp1,m
  180     a(l,j)=zero
  190  continue
c
       w(j)=zero
c                   solve the triangular system
c                  store the solution temporarily in zz()
       assign 200 to next
       go to 400
  200  continue
c              ******** secondary loop begins here *******
c                  iteration counter.
c
  210  iter=iter+1
       if (iter.le.itmax) go to 220
       mode=3
       write(6,440)
       go to 350
  220  continue
c
c   see if all new constrained coeffs are feasible.
c            if not compute alpha
       alpha=two
         do 240 ip=1,nsetp
         l=index(ip)
         if (zz(ip)) 230,230,240
c
  230    t=-x(l)/(zz(ip)-x(l))
         if (alpha.le.t) go to 240
         alpha=t
         jj=ip
  240    continue
c
c     if all new constrained coeffs are feasible then alpha
c     will atill =2.  if so exit from secondary loop to main 
c         loop.
        if (alpha.eq.two) go to 330
c
c          otherwise use alpha which will be between 0. and 1. 
c          to interpolate between the old x and the new zz.
c
        do 250 ip=1,nsetp
        l=index(ip)
  250   x(l)=x(l)+alpha*(zz(ip)-x(l))
c
c    modify a and b and the index arrays to move coefficient 
c     i from set p to set z. 
c
        i=index(jj)
  260   x(i)=zero
c
        if (jj.eq.nsetp) go to 290
        jj=jj+1
           do 280 j=jj,nsetp
           ii=index(j)
           index(j-1)=ii
           call g1(a(j-1,ii),a(j,ii),cc,ss,a(j-1,ii))
           a(j,ii)=zero
              do 270 l=1,n
              if (l.ne.ii) call g2(cc,ss,a(j-1,l),a(j,l))
  270         continue
  280      call g2(cc,ss,b(j-1),b(j))
  290   npp1=nsetp
        nsetp=nsetp-1
        iz1=iz1-1
        index(iz1)=i
c
c   see if the remaining coeffs in set p are feasible. they
c   should be because of the way alpha was determined. if 
c   any are infeasible it is due to round-off error. any that
c   are nonpositive will be set to zero and moved from set p
c   to set z.
c
        do 300 jj=1,nsetp
        i=index(jj)
        if (x(i)) 260,260,300
  300   continue
c
c     copy b() into zz(). then solve again and loop back.
c
        do 310 i=1,m
  310   zz(i)=b(i)
        assign 320 to next
        go to 400
  320    continue
        go to 210
c           ***** end of secondary loop ******
c
  330      do 340 ip=1,nsetp
           i=index(ip)
  340      x(i)=zz(ip)
c     all new coeffs are positive. loop back to beginning.
        go to 30
c
c       ***** end of main loop *****
c
c          come to here for termination.
c       compute the norm of the final residual vector.
c
  350   sm=zero
        if (npp1.gt.m) go to 370
           do 360 i=npp1,m
  360      sm=sm+b(i)**2
        go to 390
  370      do 380 j=1,n
  380      w(j)=zero
  390   rnorm=sqrt(sm)
        return
c
c   the following block of code is used as an internal subroutine
c   to solve the triangular system. putting the solution in zz().
c
  400   do 430 l=1,nsetp
        ip=nsetp+1-l
        if (l.eq.1) go to 420
            do 410 ii=1,ip
  410       zz(ii)=zz(ii)-a(ii,jj)*zz(ip+1)
  420   jj=index(ip)
  430   zz(ip)=zz(ip)/a(ip,jj)
        go to next, (200,320)
  440   format(35h0 nnls quitting on iteration count.)
        end
c
c*****************************************************
c
       subroutine h12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
c
c      construction and/or application of a single
c      householder transformation.  q = i + u*(u**t)/b
c
c     mode      =1 or 2  to select algorithm h1 or h2.
c     lpivot   is the index of the pivot element.
c     l1,m   if l1.le.m   the transformation will be constructed 
c            to zero elements indexed from l1 through m.   if
c            l1.gt.m the subroutine does an identity transformation.
c     u(),iue,up    on entry to h1 u() contains the pivot vector.
c                   iue is the storage increment between elements.
c                     on exit from h1 u() and up contain quantities 
c                   defining the vector u of the householder 
c                   transformation. on entry to h2 u() and up should
c                   contain quantities previously computed by h1. 
c                   this will not be modified by h2.
c     c()    on entry to h1 or h2 c() contains a matrix which will
c            be regarded as a set of vectors to which the householder 
c            transformation is to be applied. on exit c() contains 
c            the set of transformed vectors.
c     ice    storage increment between elements of vectors in c()
c     icv    storage increment between vectors in c().
c     ncv    number of vectors in c() to be transformed. if ncv.le.0
c            no operations will be done on c().
c
       dimension u(iue,m),c(*)
       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 single 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
c
c*******************************************
c
       subroutine g1(a,b,cos,sin,sig)
c
c   compute orthogonal rotation matrix
c   compute.. matrix (c,s) so that (c,s)(a) = (sqrt(a**2+b**2))
c                    (-s,c)       (-s,c)(b)   (   0          )
c   compute sig = sqrt(a**2+b**2)
c   sig is computed last to allow for the possibility that sig
c   may be in the same location as a or b.
c
       zero=0.
       one=1.
       if (abs(a).le.abs(b)) go to 10
       xr=b/a
       yr=sqrt(one+xr**2)
       cos=sign(one/yr,a)
       sin=cos*xr
       sig=abs(a)*yr
       return
  10   if (b) 20,30,20
  20   xr=a/b
       yr=sqrt(one+xr**2)
       sin=sign(one/yr,b)
       cos=sin*xr
       sig=abs(b)*yr
       return
  30   sig=zero
       cos=zero
       sin=one
       return
       end
c
c******************************************
c
       subroutine g2(cos,sin,x,y)
c
c      apply the rotation computed by g1 to (x,y).
c
       xr=cos*x+sin*y
       y=-sin*x+cos*y
       x=xr
       return
       end
c
c****************************************************
c
       function diff(x,y)
c
       diff=x-y
       return
       end
c
c**************************************************