diff options
Diffstat (limited to 'math/llsq/original_f/nnls.f')
| -rw-r--r-- | math/llsq/original_f/nnls.f | 278 |
1 files changed, 278 insertions, 0 deletions
diff --git a/math/llsq/original_f/nnls.f b/math/llsq/original_f/nnls.f new file mode 100644 index 00000000..fb5aa767 --- /dev/null +++ b/math/llsq/original_f/nnls.f @@ -0,0 +1,278 @@ +c subroutine nnls (a,mda,m,n,b,x,rnorm,w,zz,index,mode) +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 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 + subroutine nnls (a,mda,m,n,b,x,rnorm,w,zz,index,mode) + 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 + 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) 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 + 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. +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) + 440 format (35h0 nnls quitting on iteration count.) + end |