Initial commit

1 files changed, 276 insertions, 0 deletions

diff --git a/math/llsq/nnls.f b/math/llsq/nnls.f
new file mode 100644
index 00000000..58337147
--- /dev/null
+++ b/math/llsq/nnls.f
@@ -0,0 +1,276 @@
+c     subroutine nnls  (a,mda,m,n,b,x,rnorm,w,zz,index,ier)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c	  **********   nonnegative least squares   **********
+c
+c     given an m by n matrix, a, and an m-vector, b,  compute an
+c     n-vector, x, which solves the least squares problem
+c
+c		       a * x = b  subject to x .ge. 0
+c
+c     a(),mda,m,n     mda is the first dimensioning parameter for the
+c		      array, a().   on entry a() contains the m by n
+c		      matrix, a.	   on exit a() contains
+c		      the product matrix, q*a , where q is an
+c		      m by m orthogonal matrix generated implicitly by
+c		      this subroutine.
+c     b()     on entry b() contains the m-vector, b.   on exit b() con-
+c	      tains q*b.
+c     x()     on entry x() need not be initialized.  on exit x() will
+c	      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 vector.   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 sets
+c		  p and z as follows..
+c
+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     ier     this is a success-failure flag with the following
+c	      meanings.
+c	      0     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
+      subroutine nnls (a,mda,m,n,b,x,rnorm,w,zz,index,ier)
+      dimension a(mda,n), b(m), x(n), w(n), zz(m)
+      integer index(n)
+      zero=0.
+      one=1.
+      two=2.
+      factor=0.01
+c
+      ier=0
+      if (m.gt.0.and.n.gt.0) go to 10
+      ier=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) vector w().
+c
+	  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 transformation and check new diagonal element to avoid
+c     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 zz and
+c   > solve for ztest ( = proposed new value for x(j) ).
+c
+  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     reject j as a candidate to be moved from set z to set p.
+c     restore a(npp1,j), set w(j)=0., and loop back to test dual
+c
+      if (ztest) 130,130,140
+c
+c     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
+c     set z to set p.	 update b,  update indices,  apply householder
+c     transformations to cols in new set z,  zero subdiagonal elts in
+c     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
+c			******	secondary loop begins here ******
+c
+c			   iteration counter.
+c
+  210 iter=iter+1
+      if (iter.le.itmax) go to 220
+         ier=3
+         go to 350
+  220 continue
+c
+c		     see if all new constrained coeffs are feasible.
+c				   if not compute alpha.
+c
+      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 will
+c	   still = 2.	 if so exit from secondary loop to main loop.
+c
+      if (alpha.eq.two) go to 330
+c
+c	   otherwise use alpha which will be between 0. and 1. to
+c	   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 i
+c	 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 should
+c	 be because of the way alpha was determined.
+c	 if any are infeasible it is due to round-off error.  any
+c	 that are nonpositive will be set to zero
+c	 and moved from set p 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)
+      end