From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/llsq/README | 33 +++++ math/llsq/bndacc.f | 74 +++++++++++ math/llsq/bndsol.f | 71 +++++++++++ math/llsq/diff.f | 6 + math/llsq/g1.f | 33 +++++ math/llsq/g2.f | 9 ++ math/llsq/gen.f | 31 +++++ math/llsq/h12.f | 80 ++++++++++++ math/llsq/hfti.f | 136 +++++++++++++++++++++ math/llsq/ldp.f | 79 ++++++++++++ math/llsq/mfeout.f | 64 ++++++++++ math/llsq/mkpkg | 23 ++++ math/llsq/nnls.f | 276 +++++++++++++++++++++++++++++++++++++++++ math/llsq/original_f/bndacc.f | 74 +++++++++++ math/llsq/original_f/bndsol.f | 70 +++++++++++ math/llsq/original_f/diff.f | 6 + math/llsq/original_f/g1.f | 33 +++++ math/llsq/original_f/g2.f | 9 ++ math/llsq/original_f/gen.f | 28 +++++ math/llsq/original_f/h12.f | 80 ++++++++++++ math/llsq/original_f/hfti.f | 136 +++++++++++++++++++++ math/llsq/original_f/ldp.f | 79 ++++++++++++ math/llsq/original_f/nnls.f | 278 ++++++++++++++++++++++++++++++++++++++++++ math/llsq/original_f/qrbd.f | 208 +++++++++++++++++++++++++++++++ math/llsq/original_f/sfeout.f | 64 ++++++++++ math/llsq/original_f/sva.f | 193 +++++++++++++++++++++++++++++ math/llsq/original_f/svdrs.f | 205 +++++++++++++++++++++++++++++++ math/llsq/progs/README | 5 + math/llsq/progs/band.x | 70 +++++++++++ math/llsq/progs/data4 | 15 +++ math/llsq/progs/lsq.x | 70 +++++++++++ math/llsq/progs/prog1.f | 124 +++++++++++++++++++ math/llsq/progs/prog2.f | 125 +++++++++++++++++++ math/llsq/progs/prog3.f | 138 +++++++++++++++++++++ math/llsq/progs/prog4.f | 22 ++++ math/llsq/progs/prog5.f | 146 ++++++++++++++++++++++ math/llsq/progs/prog6.f | 116 ++++++++++++++++++ math/llsq/qrbd.f | 208 +++++++++++++++++++++++++++++++ math/llsq/sva.f | 198 ++++++++++++++++++++++++++++++ math/llsq/svdrs.f | 205 +++++++++++++++++++++++++++++++ 40 files changed, 3820 insertions(+) create mode 100644 math/llsq/README create mode 100644 math/llsq/bndacc.f create mode 100644 math/llsq/bndsol.f create mode 100644 math/llsq/diff.f create mode 100644 math/llsq/g1.f create mode 100644 math/llsq/g2.f create mode 100644 math/llsq/gen.f create mode 100644 math/llsq/h12.f create mode 100644 math/llsq/hfti.f create mode 100644 math/llsq/ldp.f create mode 100644 math/llsq/mfeout.f create mode 100644 math/llsq/mkpkg create mode 100644 math/llsq/nnls.f create mode 100644 math/llsq/original_f/bndacc.f create mode 100644 math/llsq/original_f/bndsol.f create mode 100644 math/llsq/original_f/diff.f create mode 100644 math/llsq/original_f/g1.f create mode 100644 math/llsq/original_f/g2.f create mode 100644 math/llsq/original_f/gen.f create mode 100644 math/llsq/original_f/h12.f create mode 100644 math/llsq/original_f/hfti.f create mode 100644 math/llsq/original_f/ldp.f create mode 100644 math/llsq/original_f/nnls.f create mode 100644 math/llsq/original_f/qrbd.f create mode 100644 math/llsq/original_f/sfeout.f create mode 100644 math/llsq/original_f/sva.f create mode 100644 math/llsq/original_f/svdrs.f create mode 100644 math/llsq/progs/README create mode 100644 math/llsq/progs/band.x create mode 100644 math/llsq/progs/data4 create mode 100644 math/llsq/progs/lsq.x create mode 100644 math/llsq/progs/prog1.f create mode 100644 math/llsq/progs/prog2.f create mode 100644 math/llsq/progs/prog3.f create mode 100644 math/llsq/progs/prog4.f create mode 100644 math/llsq/progs/prog5.f create mode 100644 math/llsq/progs/prog6.f create mode 100644 math/llsq/qrbd.f create mode 100644 math/llsq/sva.f create mode 100644 math/llsq/svdrs.f (limited to 'math/llsq') diff --git a/math/llsq/README b/math/llsq/README new file mode 100644 index 00000000..bfaabb58 --- /dev/null +++ b/math/llsq/README @@ -0,0 +1,33 @@ +This directory contains a collection of routines for solving linear least +squares problems by the Singular Value Decomposition (SVD) method, as +described in "Solving Least Squares Problems", by Charles L. Lawson and +Richard J. Hanson, Prentice Hall, 1974. The appendix of this book contains +full listings of the Fortran codes, as well as a users guide. + +The numerical subroutines are in this directory. The directory "progs" +contains a number of examples of the use of these routines in Fortran programs. + +The numerical routines have been modified to eliminate the use of Fortran +i/o for error conditions. The integer status return IER has been added to +all such routines, and the Fortran write statement(s) removed. A successful +call returns zero in IER, while an unsucessful call returns a positive integer +error code, identifying the error. The original codes are in the directory +"original_f". + +The affected routines and the new calling sequences are as follows: + + subroutine BNDSOL (mode,g,mdg,nb,ip,ir,x,n,rnorm,ier) + subroutine LDP (g,mdg,m,n,h,x,xnorm,w,index,ier) + subroutine NNLS (a,mda,m,n,b,x,rnorm,w,zz,index,ier) + subroutine SVDRS (a,mda,mm,nn,b,mdb,nb,s,ier) + +The routines SVA and MFEOUT were not installed in the library, since they +do extensive i/o, but have been modified to reflect the changes to the +above subroutines. + +See "lsq.x" and "band.x" in progs/ for examples demonstrating the use +of these routines in IRAF spp programs. + +20Nov82 D.Tody + +Oct85 Added comma after P edit descriptor in sva.f diff --git a/math/llsq/bndacc.f b/math/llsq/bndacc.f new file mode 100644 index 00000000..80b9914c --- /dev/null +++ b/math/llsq/bndacc.f @@ -0,0 +1,74 @@ + subroutine bndacc (g,mdg,nb,ip,ir,mt,jt) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c sequential algorithm for banded least squares problem.. +c accumulation phase. for solution phase use bndsol. +c +c the calling program must set ir=1 and ip=1 before the first call +c to bndacc for a new case. +c +c the second subscript of g( ) must be dimensioned at least +c nb+1 in the calling program. + dimension g(mdg,1) + zero=0. +c +c alg. steps 1-4 are performed external to this subroutine. +c + nbp1=nb+1 + if (mt.le.0) return +c alg. step 5 + if (jt.eq.ip) go to 70 +c alg. steps 6-7 + if (jt.le.ir) go to 30 +c alg. steps 8-9 + do 10 i=1,mt + ig1=jt+mt-i + ig2=ir+mt-i + do 10 j=1,nbp1 + 10 g(ig1,j)=g(ig2,j) +c alg. step 10 + ie=jt-ir + do 20 i=1,ie + ig=ir+i-1 + do 20 j=1,nbp1 + 20 g(ig,j)=zero +c alg. step 11 + ir=jt +c alg. step 12 + 30 mu=min0(nb-1,ir-ip-1) + if (mu.eq.0) go to 60 +c alg. step 13 + do 50 l=1,mu +c alg. step 14 + k=min0(l,jt-ip) +c alg. step 15 + lp1=l+1 + ig=ip+l + do 40 i=lp1,nb + jg=i-k + 40 g(ig,jg)=g(ig,i) +c alg. step 16 + do 50 i=1,k + jg=nbp1-i + 50 g(ig,jg)=zero +c alg. step 17 + 60 ip=jt +c alg. steps 18-19 + 70 mh=ir+mt-ip + kh=min0(nbp1,mh) +c alg. step 20 + do 80 i=1,kh + 80 call h12 (1,i,max0(i+1,ir-ip+1),mh,g(ip,i),1,rho, + 1 g(ip,i+1),1,mdg,nbp1-i) +c alg. step 21 + ir=ip+kh +c alg. step 22 + if (kh.lt.nbp1) go to 100 +c alg. step 23 + do 90 i=1,nb + 90 g(ir-1,i)=zero +c alg. step 24 + 100 continue +c alg. step 25 + return + end diff --git a/math/llsq/bndsol.f b/math/llsq/bndsol.f new file mode 100644 index 00000000..4c4bd973 --- /dev/null +++ b/math/llsq/bndsol.f @@ -0,0 +1,71 @@ + subroutine bndsol (mode,g,mdg,nb,ip,ir,x,n,rnorm,ier) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c sequential solution of a banded least squares problem.. +c solution phase. for the accumulation phase use bndacc. +c +c mode = 1 solve r*x=y where r and y are in the g( ) array +c and x will be stored in the x( ) array. +c 2 solve (r**t)*x=y where r is in g( ), +c y is initially in x( ), and x replaces y in x( ), +c 3 solve r*x=y where r is in g( ). +c y is initially in x( ), and x replaces y in x( ). +c +c the second subscript of g( ) must be dimensioned at least +c nb+1 in the calling program. + dimension g(mdg,1),x(n) + zero=0. + ier = 0 +c + rnorm=zero + go to (10,90,50), mode +c ********************* mode = 1 +c algc step 26 + 10 do 20 j=1,n + 20 x(j)=g(j,nb+1) + rsq=zero + np1=n+1 + irm1=ir-1 + if (np1.gt.irm1) go to 40 + do 30 j=np1,irm1 + 30 rsq=rsq+g(j,nb+1)**2 + rnorm=sqrt(rsq) + 40 continue +c ********************* mode = 3 +c alg. step 27 + 50 do 80 ii=1,n + i=n+1-ii +c alg. step 28 + s=zero + l=max0(0,i-ip) +c alg. step 29 + if (i.eq.n) go to 70 +c alg. step 30 + ie=min0(n+1-i,nb) + do 60 j=2,ie + jg=j+l + ix=i-1+j + 60 s=s+g(i,jg)*x(ix) +c alg. step 31 + 70 if (g(i,l+1)) 80,130,80 + 80 x(i)=(x(i)-s)/g(i,l+1) +c alg. step 32 + return +c ********************* mode = 2 + 90 do 120 j=1,n + s=zero + if (j.eq.1) go to 110 + i1=max0(1,j-nb+1) + i2=j-1 + do 100 i=i1,i2 + l=j-i+1+max0(0,i-ip) + 100 s=s+x(i)*g(i,l) + 110 l=max0(0,j-ip) + if (g(j,l+1)) 120,130,120 + 120 x(j)=(x(j)-s)/g(j,l+1) + return +c +c error return: zero diagonal term in bndsol + 130 ier = 1 + return + end diff --git a/math/llsq/diff.f b/math/llsq/diff.f new file mode 100644 index 00000000..79dacb92 --- /dev/null +++ b/math/llsq/diff.f @@ -0,0 +1,6 @@ + function diff(x,y) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 7 +c to appear in 'solving least squares problems', prentice-hall, 1974 + diff=x-y + return + end diff --git a/math/llsq/g1.f b/math/llsq/g1.f new file mode 100644 index 00000000..5d44c50b --- /dev/null +++ b/math/llsq/g1.f @@ -0,0 +1,33 @@ + subroutine g1 (a,b,cos,sin,sig) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c +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 +c sig 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 diff --git a/math/llsq/g2.f b/math/llsq/g2.f new file mode 100644 index 00000000..6f01a582 --- /dev/null +++ b/math/llsq/g2.f @@ -0,0 +1,9 @@ + subroutine g2 (cos,sin,x,y) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c apply the rotation computed by g1 to (x,y). + xr=cos*x+sin*y + y=-sin*x+cos*y + x=xr + return + end diff --git a/math/llsq/gen.f b/math/llsq/gen.f new file mode 100644 index 00000000..e7933d6d --- /dev/null +++ b/math/llsq/gen.f @@ -0,0 +1,31 @@ + function gen(anoise) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c generate numbers for construction of test cases. + mi=891 + mj=457 + i=5 + j=7 + aj=0. +c + if (anoise) 10,30,20 + 10 gen=0. + return +c +c the sequence of values of j is bounded between 1 and 996 +c if initial j = 1,2,3,4,5,6,7,8, or 9, the period is 332 +c + 20 j=j*mj + j=j-997*(j/997) + aj=j-498 +c the sequence of values of i is bounded between 1 and 999 +c if initial i = 1,2,3,6,7, or 9, the period will be 50 +c if initial i = 4 or 8 the period will be 25 +c if initial i = 5 the period will be 10 +c + 30 i=i*mi + i=i-1000*(i/1000) + ai=i-500 + gen=ai+aj*anoise + return + end diff --git a/math/llsq/h12.f b/math/llsq/h12.f new file mode 100644 index 00000000..81daa629 --- /dev/null +++ b/math/llsq/h12.f @@ -0,0 +1,80 @@ +c subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +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 to +c zero elements indexed from l1 through m. if l1 gt. m +c 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 +c contain quantities defining the vector u of the +c householder transformation. on entry to h2 u() +c and up should contain quantities previously computed +c by h1. these will not be modified by h2. +c c() on entry to h1 or h2 c() contains a matrix which will be +c regarded as a set of vectors to which the householder +c transformation is to be applied. on exit c() contains the +c 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 + subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) + dimension u(iue,m), c(1) + 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 sngl. 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 diff --git a/math/llsq/hfti.f b/math/llsq/hfti.f new file mode 100644 index 00000000..62dcebe4 --- /dev/null +++ b/math/llsq/hfti.f @@ -0,0 +1,136 @@ + subroutine hfti (a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c solve least squares problem using algorithm, hfti. +c + dimension a(mda,n),b(mdb,nb),h(n),g(n),rnorm(nb) + integer ip(n) + double precision sm,dzero + szero=0. + dzero=0.d0 + factor=0.001 +c + k=0 + ldiag=min0(m,n) + if (ldiag.le.0) go to 270 + do 80 j=1,ldiag + if (j.eq.1) go to 20 +c +c update squared column lengths and find lmax +c .. + lmax=j + do 10 l=j,n + h(l)=h(l)-a(j-1,l)**2 + if (h(l).gt.h(lmax)) lmax=l + 10 continue + if(diff(hmax+factor*h(lmax),hmax)) 20,20,50 +c +c compute squared column lengths and find lmax +c .. + 20 lmax=j + do 40 l=j,n + h(l)=0. + do 30 i=j,m + 30 h(l)=h(l)+a(i,l)**2 + if (h(l).gt.h(lmax)) lmax=l + 40 continue + hmax=h(lmax) +c .. +c lmax has been determined +c +c do column interchanges if needed. +c .. + 50 continue + ip(j)=lmax + if (ip(j).eq.j) go to 70 + do 60 i=1,m + tmp=a(i,j) + a(i,j)=a(i,lmax) + 60 a(i,lmax)=tmp + h(lmax)=h(j) +c +c compute the j-th transformation and apply it to a and b. +c .. + 70 call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j) + 80 call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb) +c +c determine the pseudorank, k, using the tolerance, tau. +c .. + do 90 j=1,ldiag + if (abs(a(j,j)).le.tau) go to 100 + 90 continue + k=ldiag + go to 110 + 100 k=j-1 + 110 kp1=k+1 +c +c compute the norms of the residual vectors. +c + if (nb.le.0) go to 140 + do 130 jb=1,nb + tmp=szero + if (kp1.gt.m) go to 130 + do 120 i=kp1,m + 120 tmp=tmp+b(i,jb)**2 + 130 rnorm(jb)=sqrt(tmp) + 140 continue +c special for pseudorank = 0 + if (k.gt.0) go to 160 + if (nb.le.0) go to 270 + do 150 jb=1,nb + do 150 i=1,n + 150 b(i,jb)=szero + go to 270 +c +c if the pseudorank is less than n compute householder +c decomposition of first k rows. +c .. + 160 if (k.eq.n) go to 180 + do 170 ii=1,k + i=kp1-ii + 170 call h12 (1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1) + 180 continue +c +c + if (nb.le.0) go to 270 + do 260 jb=1,nb +c +c solve the k by k triangular system. +c .. + do 210 l=1,k + sm=dzero + i=kp1-l + if (i.eq.k) go to 200 + ip1=i+1 + do 190 j=ip1,k + 190 sm=sm+a(i,j)*dble(b(j,jb)) + 200 sm1=sm + 210 b(i,jb)=(b(i,jb)-sm1)/a(i,i) +c +c complete computation of solution vector. +c .. + if (k.eq.n) go to 240 + do 220 j=kp1,n + 220 b(j,jb)=szero + do 230 i=1,k + 230 call h12 (2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1) +c +c re-order the solution vector to compensate for the +c column interchanges. +c .. + 240 do 250 jj=1,ldiag + j=ldiag+1-jj + if (ip(j).eq.j) go to 250 + l=ip(j) + tmp=b(l,jb) + b(l,jb)=b(j,jb) + b(j,jb)=tmp + 250 continue + 260 continue +c .. +c the solution vectors, x, are now +c in the first n rows of the array b(,). +c + 270 krank=k + return + end diff --git a/math/llsq/ldp.f b/math/llsq/ldp.f new file mode 100644 index 00000000..a083f35b --- /dev/null +++ b/math/llsq/ldp.f @@ -0,0 +1,79 @@ + subroutine ldp (g,mdg,m,n,h,x,xnorm,w,index,ier) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c +c ********** least distance programming ********** +c + integer index(m) + dimension g(mdg,n), h(m), x(n), w(1) + zero=0. + one=1. + if (n.le.0) go to 120 + do 10 j=1,n + 10 x(j)=zero + xnorm=zero + if (m.le.0) go to 110 +c +c the declared dimension of w() must be at least (n+1)*(m+2)+2*m. +c +c first (n+1)*m locs of w() = matrix e for problem nnls. +c next n+1 locs of w() = vector f for problem nnls. +c next n+1 locs of w() = vector z for problem nnls. +c next m locs of w() = vector y for problem nnls. +c next m locs of w() = vector wdual for problem nnls. +c copy g**t into first n rows and m columns of e. +c copy h**t into row n+1 of e. +c + iw=0 + do 30 j=1,m + do 20 i=1,n + iw=iw+1 + 20 w(iw)=g(j,i) + iw=iw+1 + 30 w(iw)=h(j) + if=iw+1 +c store n zeros followed by a one into f. + do 40 i=1,n + iw=iw+1 + 40 w(iw)=zero + w(iw+1)=one +c + np1=n+1 + iz=iw+2 + iy=iz+np1 + iwdual=iy+m +c + call nnls (w,np1,np1,m,w(if),w(iy),rnorm,w(iwdual),w(iz),index, + * ier) +c use the following return if unsuccessful in nnls. + if (ier.ne.0) return + if (rnorm) 130,130,50 + 50 fac=one + iw=iy-1 + do 60 i=1,m + iw=iw+1 +c here we are using the solution vector y. + 60 fac=fac-h(i)*w(iw) +c + if (diff(one+fac,one)) 130,130,70 + 70 fac=one/fac + do 90 j=1,n + iw=iy-1 + do 80 i=1,m + iw=iw+1 +c here we are using the solution vector y. + 80 x(j)=x(j)+g(i,j)*w(iw) + 90 x(j)=x(j)*fac + do 100 j=1,n + 100 xnorm=xnorm+x(j)**2 + xnorm=sqrt(xnorm) +c successful return. + 110 ier=0 + return +c error return. n .le. 0. + 120 ier=2 + return +c returning with constraints not compatible. + 130 ier=4 + return + end diff --git a/math/llsq/mfeout.f b/math/llsq/mfeout.f new file mode 100644 index 00000000..fb439e40 --- /dev/null +++ b/math/llsq/mfeout.f @@ -0,0 +1,64 @@ + subroutine mfeout (a,mda,m,n,names,mode) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c subroutine for matrix output with labeling. +c +c a( ) matrix to be output +c mda first dimension of a array +c m no. of rows in a matrix +c n no. of cols in a matrix +c names() array of names. if names(1) = 1h , the rest +c of the names() array will be ignored. +c mode =1 for 4p8f15.0 format for v matrix. +c =2 for 8e15.8 format for candidate solutions. +c + dimension a(mda,01) + integer names(m),ihead(2) + logical notblk + data maxcol/8/, iblank/1h /,ihead(1)/4h col/,ihead(2)/4hsoln/ +c + notblk=names(1).ne.iblank + if (m.le.0.or.n.le.0) return +c + if (mode.eq.2) go to 10 + write (6,70) + go to 20 + 10 write (6,80) + 20 continue +c + nblock=n/maxcol + last=n-nblock*maxcol + ncol=maxcol + j1=1 +c +c main loop starts here +c + 30 if (nblock.gt.0) go to 40 + if (last.le.0) return + ncol=last + last=0 +c + 40 j2=j1+ncol-1 + write (6,90) (ihead(mode),j,j=j1,j2) +c + do 60 i=1,m + name=iblank + if (notblk) name=names(i) +c + if (mode.eq.2) go to 50 + write (6,100) i,name,(a(i,j),j=j1,j2) + go to 60 + 50 write (6,110) i,name,(a(i,j),j=j1,j2) + 60 continue +c + j1=j1+maxcol + nblock=nblock-1 + go to 30 +c + 70 format (45h0v-matrix of the singular value decomposition, + * 8h of a*d./47h (elements of v scaled up by a factor of 10**4)) + 80 format (35h0sequence of candidate solutions, x) + 90 format (1h0,11x,8(6x,a4,i4,1x)/1x) + 100 format (1x,i3,1x,a6,1x,4p8f15.0) + 110 format (1x,i3,1x,a6,1x,8e15.8) + end diff --git a/math/llsq/mkpkg b/math/llsq/mkpkg new file mode 100644 index 00000000..1d2cf4a5 --- /dev/null +++ b/math/llsq/mkpkg @@ -0,0 +1,23 @@ +# Update Lawson's and Hanson's linear least squares package (LIBLLSQ). + +$checkout libllsq.a lib$ +$update libllsq.a +$checkin libllsq.a lib$ +$exit + +libllsq.a: + bndacc.f + bndsol.f + diff.f + g1.f + g2.f + gen.f + h12.f + hfti.f + ldp.f + mfeout.f + nnls.f + qrbd.f + sva.f + svdrs.f + ; 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 diff --git a/math/llsq/original_f/bndacc.f b/math/llsq/original_f/bndacc.f new file mode 100644 index 00000000..80b9914c --- /dev/null +++ b/math/llsq/original_f/bndacc.f @@ -0,0 +1,74 @@ + subroutine bndacc (g,mdg,nb,ip,ir,mt,jt) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c sequential algorithm for banded least squares problem.. +c accumulation phase. for solution phase use bndsol. +c +c the calling program must set ir=1 and ip=1 before the first call +c to bndacc for a new case. +c +c the second subscript of g( ) must be dimensioned at least +c nb+1 in the calling program. + dimension g(mdg,1) + zero=0. +c +c alg. steps 1-4 are performed external to this subroutine. +c + nbp1=nb+1 + if (mt.le.0) return +c alg. step 5 + if (jt.eq.ip) go to 70 +c alg. steps 6-7 + if (jt.le.ir) go to 30 +c alg. steps 8-9 + do 10 i=1,mt + ig1=jt+mt-i + ig2=ir+mt-i + do 10 j=1,nbp1 + 10 g(ig1,j)=g(ig2,j) +c alg. step 10 + ie=jt-ir + do 20 i=1,ie + ig=ir+i-1 + do 20 j=1,nbp1 + 20 g(ig,j)=zero +c alg. step 11 + ir=jt +c alg. step 12 + 30 mu=min0(nb-1,ir-ip-1) + if (mu.eq.0) go to 60 +c alg. step 13 + do 50 l=1,mu +c alg. step 14 + k=min0(l,jt-ip) +c alg. step 15 + lp1=l+1 + ig=ip+l + do 40 i=lp1,nb + jg=i-k + 40 g(ig,jg)=g(ig,i) +c alg. step 16 + do 50 i=1,k + jg=nbp1-i + 50 g(ig,jg)=zero +c alg. step 17 + 60 ip=jt +c alg. steps 18-19 + 70 mh=ir+mt-ip + kh=min0(nbp1,mh) +c alg. step 20 + do 80 i=1,kh + 80 call h12 (1,i,max0(i+1,ir-ip+1),mh,g(ip,i),1,rho, + 1 g(ip,i+1),1,mdg,nbp1-i) +c alg. step 21 + ir=ip+kh +c alg. step 22 + if (kh.lt.nbp1) go to 100 +c alg. step 23 + do 90 i=1,nb + 90 g(ir-1,i)=zero +c alg. step 24 + 100 continue +c alg. step 25 + return + end diff --git a/math/llsq/original_f/bndsol.f b/math/llsq/original_f/bndsol.f new file mode 100644 index 00000000..05dd35d7 --- /dev/null +++ b/math/llsq/original_f/bndsol.f @@ -0,0 +1,70 @@ + subroutine bndsol (mode,g,mdg,nb,ip,ir,x,n,rnorm) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c sequential solution of a banded least squares problem.. +c solution phase. for the accumulation phase use bndacc. +c +c mode = 1 solve r*x=y where r and y are in the g( ) array +c and x will be stored in the x( ) array. +c 2 solve (r**t)*x=y where r is in g( ), +c y is initially in x( ), and x replaces y in x( ), +c 3 solve r*x=y where r is in g( ). +c y is initially in x( ), and x replaces y in x( ). +c +c the second subscript of g( ) must be dimensioned at least +c nb+1 in the calling program. + dimension g(mdg,1),x(n) + zero=0. +c + rnorm=zero + go to (10,90,50), mode +c ********************* mode = 1 +c algc step 26 + 10 do 20 j=1,n + 20 x(j)=g(j,nb+1) + rsq=zero + np1=n+1 + irm1=ir-1 + if (np1.gt.irm1) go to 40 + do 30 j=np1,irm1 + 30 rsq=rsq+g(j,nb+1)**2 + rnorm=sqrt(rsq) + 40 continue +c ********************* mode = 3 +c alg. step 27 + 50 do 80 ii=1,n + i=n+1-ii +c alg. step 28 + s=zero + l=max0(0,i-ip) +c alg. step 29 + if (i.eq.n) go to 70 +c alg. step 30 + ie=min0(n+1-i,nb) + do 60 j=2,ie + jg=j+l + ix=i-1+j + 60 s=s+g(i,jg)*x(ix) +c alg. step 31 + 70 if (g(i,l+1)) 80,130,80 + 80 x(i)=(x(i)-s)/g(i,l+1) +c alg. step 32 + return +c ********************* mode = 2 + 90 do 120 j=1,n + s=zero + if (j.eq.1) go to 110 + i1=max0(1,j-nb+1) + i2=j-1 + do 100 i=i1,i2 + l=j-i+1+max0(0,i-ip) + 100 s=s+x(i)*g(i,l) + 110 l=max0(0,j-ip) + if (g(j,l+1)) 120,130,120 + 120 x(j)=(x(j)-s)/g(j,l+1) + return +c + 130 write (6,140) mode,i,j,l + stop + 140 format (30h0zero diagonal term in bndsol.,12h mode,i,j,l=,4i6) + end diff --git a/math/llsq/original_f/diff.f b/math/llsq/original_f/diff.f new file mode 100644 index 00000000..79dacb92 --- /dev/null +++ b/math/llsq/original_f/diff.f @@ -0,0 +1,6 @@ + function diff(x,y) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 7 +c to appear in 'solving least squares problems', prentice-hall, 1974 + diff=x-y + return + end diff --git a/math/llsq/original_f/g1.f b/math/llsq/original_f/g1.f new file mode 100644 index 00000000..5d44c50b --- /dev/null +++ b/math/llsq/original_f/g1.f @@ -0,0 +1,33 @@ + subroutine g1 (a,b,cos,sin,sig) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c +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 +c sig 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 diff --git a/math/llsq/original_f/g2.f b/math/llsq/original_f/g2.f new file mode 100644 index 00000000..6f01a582 --- /dev/null +++ b/math/llsq/original_f/g2.f @@ -0,0 +1,9 @@ + subroutine g2 (cos,sin,x,y) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c apply the rotation computed by g1 to (x,y). + xr=cos*x+sin*y + y=-sin*x+cos*y + x=xr + return + end diff --git a/math/llsq/original_f/gen.f b/math/llsq/original_f/gen.f new file mode 100644 index 00000000..98181a93 --- /dev/null +++ b/math/llsq/original_f/gen.f @@ -0,0 +1,28 @@ + function gen(anoise) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c generate numbers for construction of test cases. + if (anoise) 10,30,20 + 10 mi=891 + mj=457 + i=5 + j=7 + aj=0. + gen=0. + return +c +c the sequence of values of j is bounded between 1 and 996 +c if initial j = 1,2,3,4,5,6,7,8, or 9, the period is 332 + 20 j=j*mj + j=j-997*(j/997) + aj=j-498 +c the sequence of values of i is bounded between 1 and 999 +c if initial i = 1,2,3,6,7, or 9, the period will be 50 +c if initial i = 4 or 8 the period will be 25 +c if initial i = 5 the period will be 10 + 30 i=i*mi + i=i-1000*(i/1000) + ai=i-500 + gen=ai+aj*anoise + return + end diff --git a/math/llsq/original_f/h12.f b/math/llsq/original_f/h12.f new file mode 100644 index 00000000..81daa629 --- /dev/null +++ b/math/llsq/original_f/h12.f @@ -0,0 +1,80 @@ +c subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +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 to +c zero elements indexed from l1 through m. if l1 gt. m +c 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 +c contain quantities defining the vector u of the +c householder transformation. on entry to h2 u() +c and up should contain quantities previously computed +c by h1. these will not be modified by h2. +c c() on entry to h1 or h2 c() contains a matrix which will be +c regarded as a set of vectors to which the householder +c transformation is to be applied. on exit c() contains the +c 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 + subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) + dimension u(iue,m), c(1) + 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 sngl. 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 diff --git a/math/llsq/original_f/hfti.f b/math/llsq/original_f/hfti.f new file mode 100644 index 00000000..62dcebe4 --- /dev/null +++ b/math/llsq/original_f/hfti.f @@ -0,0 +1,136 @@ + subroutine hfti (a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c solve least squares problem using algorithm, hfti. +c + dimension a(mda,n),b(mdb,nb),h(n),g(n),rnorm(nb) + integer ip(n) + double precision sm,dzero + szero=0. + dzero=0.d0 + factor=0.001 +c + k=0 + ldiag=min0(m,n) + if (ldiag.le.0) go to 270 + do 80 j=1,ldiag + if (j.eq.1) go to 20 +c +c update squared column lengths and find lmax +c .. + lmax=j + do 10 l=j,n + h(l)=h(l)-a(j-1,l)**2 + if (h(l).gt.h(lmax)) lmax=l + 10 continue + if(diff(hmax+factor*h(lmax),hmax)) 20,20,50 +c +c compute squared column lengths and find lmax +c .. + 20 lmax=j + do 40 l=j,n + h(l)=0. + do 30 i=j,m + 30 h(l)=h(l)+a(i,l)**2 + if (h(l).gt.h(lmax)) lmax=l + 40 continue + hmax=h(lmax) +c .. +c lmax has been determined +c +c do column interchanges if needed. +c .. + 50 continue + ip(j)=lmax + if (ip(j).eq.j) go to 70 + do 60 i=1,m + tmp=a(i,j) + a(i,j)=a(i,lmax) + 60 a(i,lmax)=tmp + h(lmax)=h(j) +c +c compute the j-th transformation and apply it to a and b. +c .. + 70 call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j) + 80 call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb) +c +c determine the pseudorank, k, using the tolerance, tau. +c .. + do 90 j=1,ldiag + if (abs(a(j,j)).le.tau) go to 100 + 90 continue + k=ldiag + go to 110 + 100 k=j-1 + 110 kp1=k+1 +c +c compute the norms of the residual vectors. +c + if (nb.le.0) go to 140 + do 130 jb=1,nb + tmp=szero + if (kp1.gt.m) go to 130 + do 120 i=kp1,m + 120 tmp=tmp+b(i,jb)**2 + 130 rnorm(jb)=sqrt(tmp) + 140 continue +c special for pseudorank = 0 + if (k.gt.0) go to 160 + if (nb.le.0) go to 270 + do 150 jb=1,nb + do 150 i=1,n + 150 b(i,jb)=szero + go to 270 +c +c if the pseudorank is less than n compute householder +c decomposition of first k rows. +c .. + 160 if (k.eq.n) go to 180 + do 170 ii=1,k + i=kp1-ii + 170 call h12 (1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1) + 180 continue +c +c + if (nb.le.0) go to 270 + do 260 jb=1,nb +c +c solve the k by k triangular system. +c .. + do 210 l=1,k + sm=dzero + i=kp1-l + if (i.eq.k) go to 200 + ip1=i+1 + do 190 j=ip1,k + 190 sm=sm+a(i,j)*dble(b(j,jb)) + 200 sm1=sm + 210 b(i,jb)=(b(i,jb)-sm1)/a(i,i) +c +c complete computation of solution vector. +c .. + if (k.eq.n) go to 240 + do 220 j=kp1,n + 220 b(j,jb)=szero + do 230 i=1,k + 230 call h12 (2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1) +c +c re-order the solution vector to compensate for the +c column interchanges. +c .. + 240 do 250 jj=1,ldiag + j=ldiag+1-jj + if (ip(j).eq.j) go to 250 + l=ip(j) + tmp=b(l,jb) + b(l,jb)=b(j,jb) + b(j,jb)=tmp + 250 continue + 260 continue +c .. +c the solution vectors, x, are now +c in the first n rows of the array b(,). +c + 270 krank=k + return + end diff --git a/math/llsq/original_f/ldp.f b/math/llsq/original_f/ldp.f new file mode 100644 index 00000000..65d4c77e --- /dev/null +++ b/math/llsq/original_f/ldp.f @@ -0,0 +1,79 @@ + subroutine ldp (g,mdg,m,n,h,x,xnorm,w,index,mode) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c +c ********** least distance programming ********** +c + integer index(m) + dimension g(mdg,n), h(m), x(n), w(1) + zero=0. + one=1. + if (n.le.0) go to 120 + do 10 j=1,n + 10 x(j)=zero + xnorm=zero + if (m.le.0) go to 110 +c +c the declared dimension of w() must be at least (n+1)*(m+2)+2*m. +c +c first (n+1)*m locs of w() = matrix e for problem nnls. +c next n+1 locs of w() = vector f for problem nnls. +c next n+1 locs of w() = vector z for problem nnls. +c next m locs of w() = vector y for problem nnls. +c next m locs of w() = vector wdual for problem nnls. +c copy g**t into first n rows and m columns of e. +c copy h**t into row n+1 of e. +c + iw=0 + do 30 j=1,m + do 20 i=1,n + iw=iw+1 + 20 w(iw)=g(j,i) + iw=iw+1 + 30 w(iw)=h(j) + if=iw+1 +c store n zeros followed by a one into f. + do 40 i=1,n + iw=iw+1 + 40 w(iw)=zero + w(iw+1)=one +c + np1=n+1 + iz=iw+2 + iy=iz+np1 + iwdual=iy+m +c + call nnls (w,np1,np1,m,w(if),w(iy),rnorm,w(iwdual),w(iz),index, + * mode) +c use the following return if unsuccessful in nnls. + if (mode.ne.1) return + if (rnorm) 130,130,50 + 50 fac=one + iw=iy-1 + do 60 i=1,m + iw=iw+1 +c here we are using the solution vector y. + 60 fac=fac-h(i)*w(iw) +c + if (diff(one+fac,one)) 130,130,70 + 70 fac=one/fac + do 90 j=1,n + iw=iy-1 + do 80 i=1,m + iw=iw+1 +c here we are using the solution vector y. + 80 x(j)=x(j)+g(i,j)*w(iw) + 90 x(j)=x(j)*fac + do 100 j=1,n + 100 xnorm=xnorm+x(j)**2 + xnorm=sqrt(xnorm) +c successful return. + 110 mode=1 + return +c error return. n .le. 0. + 120 mode=2 + return +c returning with constraints not compatible. + 130 mode=4 + return + end 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 diff --git a/math/llsq/original_f/qrbd.f b/math/llsq/original_f/qrbd.f new file mode 100644 index 00000000..d9a12a24 --- /dev/null +++ b/math/llsq/original_f/qrbd.f @@ -0,0 +1,208 @@ +c subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c qr algorithm for singular values of a bidiagonal matrix. +c +c the bidiagonal matrix +c +c (q1,e2,0... ) +c ( q2,e3,0... ) +c d= ( . ) +c ( . 0) +c ( .en) +c ( 0,qn) +c +c is pre and post multiplied by +c elementary rotation matrices +c ri and pi so that +c +c rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn) +c +c to within working accuracy. +c +c 1. ei and qi occupy e(i) and q(i) as input. +c +c 2. rm...r1*c replaces 'c' in storage as output. +c +c 3. v*p1**(t)...pm**(t) replaces 'v' in storage as output. +c +c 4. si occupies q(i) as output. +c +c 5. the si's are nonincreasing and nonnegative. +c +c this code is based on the paper and 'algol' code.. +c ref.. +c 1. reinsch,c.h. and golub,g.h. 'singular value decomposition +c and least squares solutions' (numer. math.), vol. 14,(1970). +c + subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) + logical wntv ,havers,fail + dimension q(nn),e(nn),v(mdv,nn),c(mdc,ncc) + zero=0. + one=1. + two=2. +c + n=nn + ipass=1 + if (n.le.0) return + n10=10*n + wntv=nrv.gt.0 + havers=ncc.gt.0 + fail=.false. + nqrs=0 + e(1)=zero + dnorm=zero + do 10 j=1,n + 10 dnorm=amax1(abs(q(j))+abs(e(j)),dnorm) + do 200 kk=1,n + k=n+1-kk +c +c test for splitting or rank deficiencies.. +c first make test for last diagonal term, q(k), being small. + 20 if(k.eq.1) go to 50 + if(diff(dnorm+q(k),dnorm)) 50,25,50 +c +c since q(k) is small we will make a special pass to +c transform e(k) to zero. +c + 25 cs=zero + sn=-one + do 40 ii=2,k + i=k+1-ii + f=-sn*e(i+1) + e(i+1)=cs*e(i+1) + call g1 (q(i),f,cs,sn,q(i)) +c transformation constructed to zero position (i,k). +c + if (.not.wntv) go to 40 + do 30 j=1,nrv + 30 call g2 (cs,sn,v(j,i),v(j,k)) +c accumulate rt. transformations in v. +c + 40 continue +c +c the matrix is now bidiagonal, and of lower order +c since e(k) .eq. zero.. +c + 50 do 60 ll=1,k + l=k+1-ll + if(diff(dnorm+e(l),dnorm)) 55,100,55 + 55 if(diff(dnorm+q(l-1),dnorm)) 60,70,60 + 60 continue +c this loop can't complete since e(1) = zero. +c + go to 100 +c +c cancellation of e(l), l.gt.1. + 70 cs=zero + sn=-one + do 90 i=l,k + f=-sn*e(i) + e(i)=cs*e(i) + if(diff(dnorm+f,dnorm)) 75,100,75 + 75 call g1 (q(i),f,cs,sn,q(i)) + if (.not.havers) go to 90 + do 80 j=1,ncc + 80 call g2 (cs,sn,c(i,j),c(l-1,j)) + 90 continue +c +c test for convergence.. + 100 z=q(k) + if (l.eq.k) go to 170 +c +c shift from bottom 2 by 2 minor of b**(t)*b. + x=q(l) + y=q(k-1) + g=e(k-1) + h=e(k) + f=((y-z)*(y+z)+(g-h)*(g+h))/(two*h*y) + g=sqrt(one+f**2) + if (f.lt.zero) go to 110 + t=f+g + go to 120 + 110 t=f-g + 120 f=((x-z)*(x+z)+h*(y/t-h))/x +c +c next qr sweep.. + cs=one + sn=one + lp1=l+1 + do 160 i=lp1,k + g=e(i) + y=q(i) + h=sn*g + g=cs*g + call g1 (f,h,cs,sn,e(i-1)) + f=x*cs+g*sn + g=-x*sn+g*cs + h=y*sn + y=y*cs + if (.not.wntv) go to 140 +c +c accumulate rotations (from the right) in 'v' + do 130 j=1,nrv + 130 call g2 (cs,sn,v(j,i-1),v(j,i)) + 140 call g1 (f,h,cs,sn,q(i-1)) + f=cs*g+sn*y + x=-sn*g+cs*y + if (.not.havers) go to 160 + do 150 j=1,ncc + 150 call g2 (cs,sn,c(i-1,j),c(i,j)) +c apply rotations from the left to +c right hand sides in 'c'.. +c + 160 continue + e(l)=zero + e(k)=f + q(k)=x + nqrs=nqrs+1 + if (nqrs.le.n10) go to 20 +c return to 'test for splitting'. +c + fail=.true. +c .. +c cutoff for convergence failure. 'nqrs' will be 2*n usually. + 170 if (z.ge.zero) go to 190 + q(k)=-z + if (.not.wntv) go to 190 + do 180 j=1,nrv + 180 v(j,k)=-v(j,k) + 190 continue +c convergence. q(k) is made nonnegative.. +c + 200 continue + if (n.eq.1) return + do 210 i=2,n + if (q(i).gt.q(i-1)) go to 220 + 210 continue + if (fail) ipass=2 + return +c .. +c every singular value is in order.. + 220 do 270 i=2,n + t=q(i-1) + k=i-1 + do 230 j=i,n + if (t.ge.q(j)) go to 230 + t=q(j) + k=j + 230 continue + if (k.eq.i-1) go to 270 + q(k)=q(i-1) + q(i-1)=t + if (.not.havers) go to 250 + do 240 j=1,ncc + t=c(i-1,j) + c(i-1,j)=c(k,j) + 240 c(k,j)=t + 250 if (.not.wntv) go to 270 + do 260 j=1,nrv + t=v(j,i-1) + v(j,i-1)=v(j,k) + 260 v(j,k)=t + 270 continue +c end of ordering algorithm. +c + if (fail) ipass=2 + return + end diff --git a/math/llsq/original_f/sfeout.f b/math/llsq/original_f/sfeout.f new file mode 100644 index 00000000..fb439e40 --- /dev/null +++ b/math/llsq/original_f/sfeout.f @@ -0,0 +1,64 @@ + subroutine mfeout (a,mda,m,n,names,mode) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c subroutine for matrix output with labeling. +c +c a( ) matrix to be output +c mda first dimension of a array +c m no. of rows in a matrix +c n no. of cols in a matrix +c names() array of names. if names(1) = 1h , the rest +c of the names() array will be ignored. +c mode =1 for 4p8f15.0 format for v matrix. +c =2 for 8e15.8 format for candidate solutions. +c + dimension a(mda,01) + integer names(m),ihead(2) + logical notblk + data maxcol/8/, iblank/1h /,ihead(1)/4h col/,ihead(2)/4hsoln/ +c + notblk=names(1).ne.iblank + if (m.le.0.or.n.le.0) return +c + if (mode.eq.2) go to 10 + write (6,70) + go to 20 + 10 write (6,80) + 20 continue +c + nblock=n/maxcol + last=n-nblock*maxcol + ncol=maxcol + j1=1 +c +c main loop starts here +c + 30 if (nblock.gt.0) go to 40 + if (last.le.0) return + ncol=last + last=0 +c + 40 j2=j1+ncol-1 + write (6,90) (ihead(mode),j,j=j1,j2) +c + do 60 i=1,m + name=iblank + if (notblk) name=names(i) +c + if (mode.eq.2) go to 50 + write (6,100) i,name,(a(i,j),j=j1,j2) + go to 60 + 50 write (6,110) i,name,(a(i,j),j=j1,j2) + 60 continue +c + j1=j1+maxcol + nblock=nblock-1 + go to 30 +c + 70 format (45h0v-matrix of the singular value decomposition, + * 8h of a*d./47h (elements of v scaled up by a factor of 10**4)) + 80 format (35h0sequence of candidate solutions, x) + 90 format (1h0,11x,8(6x,a4,i4,1x)/1x) + 100 format (1x,i3,1x,a6,1x,4p8f15.0) + 110 format (1x,i3,1x,a6,1x,8e15.8) + end diff --git a/math/llsq/original_f/sva.f b/math/llsq/original_f/sva.f new file mode 100644 index 00000000..469a90ef --- /dev/null +++ b/math/llsq/original_f/sva.f @@ -0,0 +1,193 @@ + subroutine sva (a,mda,m,n,mdata,b,sing,names,iscale,d) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c singular value analysis printout +c +c iscale set by user to 1, 2, or 3 to select column scaling +c option. +c 1 subr will use identity scaling and ignore the d() +c array. +c 2 subr will scale nonzero cols to have unit euclide +c length and will store reciprocal lengths of +c original nonzero cols in d(). +c 3 user supplies col scale factors in d(). subr +c will mult col j by d(j) and remove the scaling +c from the soln at the end. +c + dimension a(mda,n), b(m), sing(01),d(n) +c sing(3*n) + dimension names(n) + logical test + double precision sb, dzero + dzero=0.d0 + one=1. + zero=0. + if (m.le.0 .or. n.le.0) return + np1=n+1 + write (6,270) + write (6,260) m,n,mdata + go to (60,10,10), iscale +c +c apply column scaling to a +c + 10 do 50 j=1,n + a1=d(j) + go to (20,20,40), iscale + 20 sb=dzero + do 30 i=1,m + 30 sb=sb+dble(a(i,j))**2 + a1=dsqrt(sb) + if (a1.eq.zero) a1=one + a1=one/a1 + d(j)=a1 + 40 do 50 i=1,m + 50 a(i,j)=a(i,j)*a1 + write (6,280) iscale,(d(j),j=1,n) + go to 70 + 60 continue + write (6,290) + 70 continue +c +c obtain sing. value decomp. of scaled matrix. +c +c********************************************************** + call svdrs (a,mda,m,n,b,1,1,sing) +c********************************************************** +c +c print the v matrix. +c + call mfeout (a,mda,n,n,names,1) + if (iscale.eq.1) go to 90 +c +c replace v by d*v in the array a() +c + do 80 i=1,n + do 80 j=1,n + 80 a(i,j)=d(i)*a(i,j) +c +c g now in b array. v now in a array. +c +c +c obtain summary output. +c + 90 continue + write (6,220) +c +c compute cumulative sums of squares of components of +c g and store them in sing(i), i=minmn+1,...,2*minmn+1 +c + sb=dzero + minmn=min0(m,n) + minmn1=minmn+1 + if (m.eq.minmn) go to 110 + do 100 i=minmn1,m + 100 sb=sb+dble(b(i))**2 + 110 sing(2*minmn+1)=sb + do 120 jj=1,minmn + j=minmn+1-jj + sb=sb+dble(b(j))**2 + js=minmn+j + 120 sing(js)=sb + a3=sing(minmn+1) + a4=sqrt(a3/float(max0(1,mdata))) + write (6,230) a3,a4 +c + nsol=0 +c +c +c + do 160 k=1,minmn + if (sing(k).eq.zero) go to 130 + nsol=k + pi=b(k)/sing(k) + a1=one/sing(k) + a2=b(k)**2 + a3=sing(minmn1+k) + a4=sqrt(a3/float(max0(1,mdata-k))) + test=sing(k).ge.100..or.sing(k).lt..001 + if (test) write (6,240) k,sing(k),pi,a1,b(k),a2,a3,a4 + if (.not.test) write (6,250) k,sing(k),pi,a1,b(k),a2,a3,a4 + go to 140 + 130 write (6,240) k,sing(k) + pi=zero + 140 do 150 i=1,n + a(i,k)=a(i,k)*pi + 150 if (k.gt.1) a(i,k)=a(i,k)+a(i,k-1) + 160 continue +c +c compute and print values of ynorm and rnorm. +c + write (6,300) + j=0 + ysq=zero + go to 180 + 170 j=j+1 + ysq=ysq+(b(j)/sing(j))**2 + 180 ynorm=sqrt(ysq) + js=minmn1+j + rnorm=sqrt(sing(js)) + yl=-1000. + if (ynorm.gt.0.) yl=alog10(ynorm) + rl=-1000. + if (rnorm.gt.0.) rl=alog10(rnorm) + write (6,310) j,ynorm,rnorm,yl,rl + if (j.lt.nsol) go to 170 +c +c compute values of xnorm and rnorm for a sequence of values of +c the levenberg-marquardt parameter. +c + if (sing(1).eq.zero) go to 210 + el=alog10(sing(1))+one + el2=alog10(sing(nsol))-one + del=(el2-el)/20. + ten=10. + aln10=alog(ten) + write (6,320) + do 200 ie=1,21 +c compute alamb=10.**el + alamb=exp(aln10*el) + ys=0. + js=minmn1+nsol + rs=sing(js) + do 190 i=1,minmn + sl=sing(i)**2+alamb**2 + ys=ys+(b(i)*sing(i)/sl)**2 + rs=rs+(b(i)*(alamb**2)/sl)**2 + 190 continue + ynorm=sqrt(ys) + rnorm=sqrt(rs) + rl=-1000. + if (rnorm.gt.zero) rl=alog10(rnorm) + yl=-1000. + if (ynorm.gt.zero) yl=alog10(ynorm) + write (6,330) alamb,ynorm,rnorm,el,yl,rl + el=el+del + 200 continue +c +c print candidate solutions. +c + 210 if (nsol.ge.1) call mfeout (a,mda,n,nsol,names,2) + return + 220 format (42h0 index sing. value p coef ,48hrecip. s. + 1v. g coef g**2 ,39h c.s.s. + 2 n.s.r.c.s.s.) + 230 format (1h ,5x,1h0,88x,1pe15.4,1pe17.4) + 240 format (1h ,i6,e12.4,1p(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4, + 12x,e15.4)) + 250 format (1h ,i6,f12.4,1p(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4, + 12x,e15.4)) + 260 format (5h0m = ,i6,8h, n = ,i4,12h, mdata = ,i8) + 270 format (45h0singular value analysis of the least squares,42h probl + 1em, a*x=b, scaled as (a*d)*y=b .) + 280 format (19h0scaling option no.,i2,18h. d is a diagonal,46h matrix + 1 with the following diagonal elements../(5x,10e12.4)) + 290 format (50h0scaling option no. 1. d is the identity matrix./1x) + 300 format (6h0index,13x,28h ynorm rnorm,14x,28h log1 + 10(ynorm) log10(rnorm)/1x) + 310 format (1h ,i4,14x,2e14.5,14x,2f14.5) + 320 format (54h0norms of solution and residual vectors for a range of, + 144h values of the levenberg-marquardt parameter,9h, lambda./1h0,4x + 2,42h lambda ynorm rnorm,42h log10(lambda) + 3log10(ynorm) log10(rnorm)) + 330 format (5x,3e14.5,3f14.5) + end diff --git a/math/llsq/original_f/svdrs.f b/math/llsq/original_f/svdrs.f new file mode 100644 index 00000000..66e488dc --- /dev/null +++ b/math/llsq/original_f/svdrs.f @@ -0,0 +1,205 @@ +c subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c singular value decomposition also treating right side vector. +c +c the array s occupies 3*n cells. +c a occupies m*n cells +c b occupies m*nb cells. +c +c special singular value decomposition subroutine. +c we have the m x n matrix a and the system a*x=b to solve. +c either m .ge. n or m .lt. n is permitted. +c the singular value decomposition +c a = u*s*v**(t) is made in such a way that one gets +c (1) the matrix v in the first n rows and columns of a. +c (2) the diagonal matrix of ordered singular values in +c the first n cells of the array s(ip), ip .ge. 3*n. +c (3) the matrix product u**(t)*b=g gets placed back in b. +c (4) the user must complete the solution and do his own +c singular value analysis. +c ******* +c give special +c treatment to rows and columns which are entirely zero. this +c causes certain zero sing. vals. to appear as exact zeros rather +c than as about eta times the largest sing. val. it similarly +c cleans up the associated columns of u and v. +c method.. +c 1. exchange cols of a to pack nonzero cols to the left. +c set n = no. of nonzero cols. +c use locations a(1,nn),a(1,nn-1),...,a(1,n+1) to record the +c col permutations. +c 2. exchange rows of a to pack nonzero rows to the top. +c quit packing if find n nonzero rows. make same row exchanges +c in b. set m so that all nonzero rows of the permuted a +c are in first m rows. if m .le. n then all m rows are +c nonzero. if m .gt. n then the first n rows are known +c to be nonzero,and rows n+1 thru m may be zero or nonzero. +c 3. apply original algorithm to the m by n problem. +c 4. move permutation record from a(,) to s(i),i=n+1,...,nn. +c 5. build v up from n by n to nn by nn by placing ones on +c the diagonal and zeros elsewhere. this is only partly done +c explicitly. it is completed during step 6. +c 6. exchange rows of v to compensate for col exchanges of step 2. +c 7. place zeros in s(i),i=n+1,nn to represent zero sing vals. +c + subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s) + dimension a(mda,nn),b(mdb,nb),s(nn,3) + zero=0. + one=1. +c +c begin.. special for zero rows and cols. +c +c pack the nonzero cols to the left +c + n=nn + if (n.le.0.or.mm.le.0) return + j=n + 10 continue + do 20 i=1,mm + if (a(i,j)) 50,20,50 + 20 continue +c +c col j is zero. exchange it with col n. +c + if (j.eq.n) go to 40 + do 30 i=1,mm + 30 a(i,j)=a(i,n) + 40 continue + a(1,n)=j + n=n-1 + 50 continue + j=j-1 + if (j.ge.1) go to 10 +c if n=0 then a is entirely zero and svd +c computation can be skipped + ns=0 + if (n.eq.0) go to 240 +c pack nonzero rows to the top +c quit packing if find n nonzero rows + i=1 + m=mm + 60 if (i.gt.n.or.i.ge.m) go to 150 + if (a(i,i)) 90,70,90 + 70 do 80 j=1,n + if (a(i,j)) 90,80,90 + 80 continue + go to 100 + 90 i=i+1 + go to 60 +c row i is zero +c exchange rows i and m + 100 if(nb.le.0) go to 115 + do 110 j=1,nb + t=b(i,j) + b(i,j)=b(m,j) + 110 b(m,j)=t + 115 do 120 j=1,n + 120 a(i,j)=a(m,j) + if (m.gt.n) go to 140 + do 130 j=1,n + 130 a(m,j)=zero + 140 continue +c exchange is finished + m=m-1 + go to 60 +c + 150 continue +c end.. special for zero rows and columns +c begin.. svd algorithm +c method.. +c (1) reduce the matrix to upper bidiagonal form with +c householder transformations. +c h(n)...h(1)aq(1)...q(n-2) = (d**t,0)**t +c where d is upper bidiagonal. +c +c (2) apply h(n)...h(1) to b. here h(n)...h(1)*b replaces b +c in storage. +c +c (3) the matrix product w= q(1)...q(n-2) overwrites the first +c n rows of a in storage. +c +c (4) an svd for d is computed. here k rotations ri and pi are +c computed so that +c rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sm) +c to working accuracy. the si are nonnegative and nonincreasing. +c here rk...r1*b overwrites b in storage while +c a*p1**(t)...pk**(t) overwrites a in storage. +c +c (5) it follows that,with the proper definitions, +c u**(t)*b overwrites b, while v overwrites the first n row and +c columns of a. +c + l=min0(m,n) +c the following loop reduces a to upper bidiagonal and +c also applies the premultiplying transformations to b. +c + do 170 j=1,l + if (j.ge.m) go to 160 + call h12 (1,j,j+1,m,a(1,j),1,t,a(1,j+1),1,mda,n-j) + call h12 (2,j,j+1,m,a(1,j),1,t,b,1,mdb,nb) + 160 if (j.ge.n-1) go to 170 + call h12 (1,j+1,j+2,n,a(j,1),mda,s(j,3),a(j+1,1),mda,1,m-j) + 170 continue +c +c copy the bidiagonal matrix into the array s() for qrbd. +c + if (n.eq.1) go to 190 + do 180 j=2,n + s(j,1)=a(j,j) + 180 s(j,2)=a(j-1,j) + 190 s(1,1)=a(1,1) +c + ns=n + if (m.ge.n) go to 200 + ns=m+1 + s(ns,1)=zero + s(ns,2)=a(m,m+1) + 200 continue +c +c construct the explicit n by n product matrix, w=q1*q2*...*ql*i +c in the array a(). +c + do 230 k=1,n + i=n+1-k + if(i.gt.min0(m,n-2)) go to 210 + call h12 (2,i+1,i+2,n,a(i,1),mda,s(i,3),a(1,i+1),1,mda,n-i) + 210 do 220 j=1,n + 220 a(i,j)=zero + 230 a(i,i)=one +c +c compute the svd of the bidiagonal matrix +c + call qrbd (ipass,s(1,1),s(1,2),ns,a,mda,n,b,mdb,nb) +c + go to (240,310), ipass + 240 continue + if (ns.ge.n) go to 260 + nsp1=ns+1 + do 250 j=nsp1,n + 250 s(j,1)=zero + 260 continue + if (n.eq.nn) return + np1=n+1 +c move record of permutations +c and store zeros + do 280 j=np1,nn + s(j,1)=a(1,j) + do 270 i=1,n + 270 a(i,j)=zero + 280 continue +c permute rows and set zero singular values. + do 300 k=np1,nn + i=s(k,1) + s(k,1)=zero + do 290 j=1,nn + a(k,j)=a(i,j) + 290 a(i,j)=zero + a(i,k)=one + 300 continue +c end.. special for zero rows and columns + return + 310 write (6,320) + stop + 320 format (49h convergence failure in qr bidiagonal svd routine) + end diff --git a/math/llsq/progs/README b/math/llsq/progs/README new file mode 100644 index 00000000..e0c5a142 --- /dev/null +++ b/math/llsq/progs/README @@ -0,0 +1,5 @@ + +This directory contains both Fortran and IRAF preprocessor programs +demonstrating the use of the Lawson-Hanson llsq procedures. The Fortran +programs have not been modified to reflect the removal of Fortran i/o +from the library procedures as installed in the IRAF llsq.a library. diff --git a/math/llsq/progs/band.x b/math/llsq/progs/band.x new file mode 100644 index 00000000..f65290ba --- /dev/null +++ b/math/llsq/progs/band.x @@ -0,0 +1,70 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + + +define MAXPTS 1024 + +task band + +# This procedure solves for the natural cubic spline which interpolates +# the curve y[i] = 100., for i = 1 to n. The matrix is of length n+2, +# and has a bandwidth of 3. Lawsons and Hansons routines BNDACC and +# BNDSOL (modified to return an error code IER) are used to solve the +# system. The computations are carried out in single precision. + + +procedure band() + +real g[MAXPTS,4], c[MAXPTS], rnorm +int ip, ir, i, n, ier, geti() +real marktime, cptime() + +begin + n = min (MAXPTS, geti ("npts")) + marktime = cptime() + + ip = 1 + ir = 1 + + g[ir,1] = 6. # first row + g[ir,2] = -12. + g[ir,3] = 6. + g[ir,4] = 0. + call bndacc (g, MAXPTS, 3, ip, ir, 1, 1) + + do i = 2, n-1 { # tridiagonal elements + g[ir,1] = 1. + g[ir,2] = 4. + g[ir,3] = 1. + g[ir,4] = 100. + call bndacc (g, MAXPTS, 3, ip, ir, 1, i-1) + } + + g[ir,1] = 6. # last row + g[ir,2] = -12. + g[ir,3] = 6. + g[ir,4] = 0. + call bndacc (g, MAXPTS, 3, ip, ir, 1, n-2) + + call printf ("matrix accumulation took %8.2f cpu seconds\n") + call pargr (cptime() - marktime) + marktime = cptime() + + # solve for the b-spline coeff C + call bndsol (1, g, MAXPTS, 3, ip, ir, c, n, rnorm, ier) + + call printf ("solution took %8.2f cpu seconds (rnorm=%g, ier=%d)\n") + call pargr (cptime() - marktime) + call pargr (rnorm) + call pargi (ier) + + call printf ("selected coefficients:\n") + for (i=1; i <= n;) { # print the first and last 4 coeff + call printf ("%8d%15.6f\n") + call pargi (i) + call pargr (c[i]) + if (i == 4) + i = max(i+1, n-3) + else + i = i + 1 + } +end diff --git a/math/llsq/progs/data4 b/math/llsq/progs/data4 new file mode 100644 index 00000000..33bf1b08 --- /dev/null +++ b/math/llsq/progs/data4 @@ -0,0 +1,15 @@ + -.13405547 -.20162827 -.16930778 -.18971990 -.17387234 -.4361 + -.10379475 -.15766336 -.13346256 -.14848550 -.13597690 -.3437 + -.08779597 -.12883867 -.10683007 -.12011796 -.10932972 -.2657 + .02058554 .00335331 -.01641270 .00078606 .00271659 -.0392 + -.03248093 -.01876799 .00410639 -.01405894 -.01384391 .0193 + .05967662 .06667714 .04352153 .05740438 .05024962 .0747 + .06712457 .07352437 .04489770 .06471862 .05876455 .0935 + .08687186 .09368296 .05672327 .08141043 .07302320 .1079 + .02149662 .06222662 .07213486 .06200069 .05570931 .1930 + .06687407 .10344506 .09153849 .09508223 .08393667 .2058 + .15879069 .18088339 .11540692 .16160727 .14796479 .2606 + .17642887 .20361830 .13057860 .18385729 .17005549 .3142 + .11414080 .17259611 .14816471 .16007466 .14374096 .3529 + .07846038 .14669563 .14365800 .14003842 .12571177 .3615 + .10803175 .16994623 .14971519 .15885312 .14301547 .3647 diff --git a/math/llsq/progs/lsq.x b/math/llsq/progs/lsq.x new file mode 100644 index 00000000..ae002a16 --- /dev/null +++ b/math/llsq/progs/lsq.x @@ -0,0 +1,70 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + + +define M 256 +define N 256 + +task lsq + +# This procedure fits a natural cubic spline to an array of n data points. +# The system being solved is a tridiagonal matrix of n+2 rows. The system +# is solved by Lawsons and Hansons routine HFTI, which solves a general +# m by n linear system of equations. This is enormous overkill for this +# problem (see "band.x"), but serves to give timing estimates for the code. + + +procedure lsq() + +real a[M,N], b[M], tau, rnorm, h[N], g[N] +int krank, ip[N] +int i, j, m, n, geti() +real marktime, cptime() + +begin + m = min (M, geti ("npts")) # size of matrix + n = min (N, m) + tau = 1e-6 + + do j = 1, n # set up b-spline matrix + do i = 1, m + a[i,j] = 0. + + a[1,1] = 6. # first row + a[1,2] = -12. + a[1,3] = 6. + + a[m,n] = 6. # last row + a[m,n-1] = -12. + a[m,n-2] = 6. + + do j = 2, m-1 { # tridiagonal elements + a[j,j-1] = 1. + a[j,j] = 4. + a[j,j+1] = 1. + } + + b[1] = 0. # natural spline bndry conditions + b[m] = 0. + + do i = 2, m-1 # set up data vector + b[i] = 100. + + marktime = cptime() + call hfti (a,M,m,n, b,1,1, tau, krank,rnorm, h,g,ip) + + call printf ("took %8.2f cpu seconds (krank=%d, rnorm=%g)\n") + call pargr (cptime() - marktime) + call pargi (krank) + call pargr (rnorm) + + call printf ("selected coefficients:\n") + for (i=1; i <= m;) { # print first, last 4 coeff + call printf ("%8d%15.5f\n") + call pargi (i) + call pargr (b[i]) + if (i == 4) + i = max(i+1, m-3) + else + i = i + 1 + } +end diff --git a/math/llsq/progs/prog1.f b/math/llsq/progs/prog1.f new file mode 100644 index 00000000..13f64929 --- /dev/null +++ b/math/llsq/progs/prog1.f @@ -0,0 +1,124 @@ +c prog1 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate algorithms hft and hs1 for solving least squares +c problems and algorithm cov for omputing the associated covariance +c matrices. +c + dimension a(8,8),h(8),b(8) + real gen,anoise + double precision sm + data mda/8/ +c + do 180 noise=1,2 + anoise=0. + if (noise.eq.2) anoise=1.e-4 + write (6,230) + write (6,240) anoise +c initialize the data generation function +c .. + dummy=gen(-1.) + do 180 mn1=1,6,5 + mn2=mn1+2 + do 180 m=mn1,mn2 + do 180 n=mn1,mn2 + np1=n+1 + write (6,250) m,n +c generate data +c .. + do 10 i=1,m + do 10 j=1,n + 10 a(i,j)=gen(anoise) + do 20 i=1,m + 20 b(i)=gen(anoise) + if(m .lt. n) go to 180 +c +c ****** begin algorithm hft ****** +c .. + do 30 j=1,n + 30 call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j) +c .. +c the algorithm 'hft' is completed. +c +c ****** begin algorithm hs1 ****** +c apply the transformations q(n)...q(1)=q to b +c replacing the previous contents of the array, b . +c .. + do 40 j=1,n + 40 call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,1,1) +c solve the triangular system for the solution x. +c store x in the array, b . +c .. + do 80 k=1,n + i=np1-k + sm=0.d0 + if (i.eq.n) go to 60 + ip1=i+1 + do 50 j=ip1,n + 50 sm=sm+a(i,j)*dble(b(j)) + 60 if (a(i,i)) 80,70,80 + 70 write (6,260) + go to 180 + 80 b(i)=(b(i)-sm)/a(i,i) +c compute length of residual vector. +c .. + srsmsq=0. + if (n.eq.m) go to 100 + mmn=m-n + do 90 j=1,mmn + npj=n+j + 90 srsmsq=srsmsq+b(npj)**2 + srsmsq=sqrt(srsmsq) +c ****** begin algorithm cov ****** +c compute unscaled covariance matrix ((a**t)*a)**(-1) +c .. + 100 do 110 j=1,n + 110 a(j,j)=1./a(j,j) + if (n.eq.1) go to 140 + nm1=n-1 + do 130 i=1,nm1 + ip1=i+1 + do 130 j=ip1,n + jm1=j-1 + sm=0.d0 + do 120 l=i,jm1 + 120 sm=sm+a(i,l)*dble(a(l,j)) + 130 a(i,j)=-sm*a(j,j) +c .. +c the upper triangle of a has been inverted +c upon itself. + 140 do 160 i=1,n + do 160 j=i,n + sm=0.d0 + do 150 l=j,n + 150 sm=sm+a(i,l)*dble(a(j,l)) + 160 a(i,j)=sm +c .. +c the upper triangular part of the +c symmetric matrix (a**t*a)**(-1) has +c replaced the upper triangular part of +c the a array. + write (6,200) (i,b(i),i=1,n) + write (6,190) srsmsq + write (6,210) + do 170 i=1,n + 170 write (6,220) (i,j,a(i,j),j=i,n) + 180 continue + stop + 190 format (1h0,8x,17hresidual length =,e12.4) + 200 format (1h0,8x,34hestimated parameters, x=a**(+)*b,,22h computed + 1by 'hft,hs1'//(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8)) + 210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa + 1rameters.,19h computed by 'cov'./1x) + 220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8) + 230 format (52h1 prog1. this program demonstrates the algorithms,19 + 1h hft, hs1, and cov.//40h caution.. 'prog1' does no checking for , + 252hnear rank deficient matrices. results in such cases,20h may be + 3 meaningless.,/34h such cases are handled by 'prog2',11h or 'pro + 4') + 240 format (1h0,54hthe relative noise level of the generated data will + 1 be,e16.4) + 250 format (1h0////9h0 m n/1x,2i4) + 260 format (1h0,8x,36h****** terminating this case due to,37h a divis + 1or being exactly zero. ******) + end diff --git a/math/llsq/progs/prog2.f b/math/llsq/progs/prog2.f new file mode 100644 index 00000000..5fb601cc --- /dev/null +++ b/math/llsq/progs/prog2.f @@ -0,0 +1,125 @@ +c prog2 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate algorithm hfti for solving least squares problems +c and algorithm cov for computing the associated unscaled +c covariance matrix. +c + dimension a(8,8),h(8),b(8),g(8) + real gen,anoise + integer ip(8) + double precision sm + data mda /8/ +c + do 180 noise=1,2 + anorm=500. + anoise=0. + tau=.5 + if (noise.eq.1) go to 10 + anoise=1.e-4 + tau=anorm*anoise*10. + 10 continue +c initialize the data generation function +c .. + dummy=gen(-1.) + write (6,230) + write (6,240) anoise,anorm,tau +c + do 180 mn1=1,6,5 + mn2=mn1+2 + do 180 m=mn1,mn2 + do 180 n=mn1,mn2 + write (6,250) m,n +c generate data +c .. + do 20 i=1,m + do 20 j=1,n + 20 a(i,j)=gen(anoise) + do 30 i=1,m + 30 b(i)=gen(anoise) +c +c ****** call hfti ****** +c + call hfti(a,mda,m,n,b,1,1,tau,krank,srsmsq,h,g,ip) +c +c + write (6,260) krank + write (6,200) (i,b(i),i=1,n) + write (6,190) srsmsq + if (krank.lt.n) go to 180 +c ****** algorithm cov bigins here ****** +c .. + do 40 j=1,n + 40 a(j,j)=1./a(j,j) + if (n.eq.1) go to 70 + nm1=n-1 + do 60 i=1,nm1 + ip1=i+1 + do 60 j=ip1,n + jm1=j-1 + sm=0.d0 + do 50 l=i,jm1 + 50 sm=sm+a(i,l)*dble(a(l,j)) + 60 a(i,j)=-sm*a(j,j) +c .. +c the upper triangle of a has been inverted +c upon itself. + 70 do 90 i=1,n + do 90 j=i,n + sm=0.d0 + do 80 l=j,n + 80 sm=sm+a(i,l)*dble(a(j,l)) + 90 a(i,j)=sm + if (n.lt.2) go to 160 + do 150 ii=2,n + i=n+1-ii + if (ip(i).eq.i) go to 150 + k=ip(i) + tmp=a(i,i) + a(i,i)=a(k,k) + a(k,k)=tmp + if (i.eq.1) go to 110 + do 100 l=2,i + tmp=a(l-1,i) + a(l-1,i)=a(l-1,k) + 100 a(l-1,k)=tmp + 110 ip1=i+1 + km1=k-1 + if (ip1.gt.km1) go to 130 + do 120 l=ip1,km1 + tmp=a(i,l) + a(i,l)=a(l,k) + 120 a(l,k)=tmp + 130 if (k.eq.n) go to 150 + kp1=k+1 + do 140 l=kp1,n + tmp=a(i,l) + a(i,l)=a(k,l) + 140 a(k,l)=tmp + 150 continue + 160 continue +c .. +c covariance has been computed and repermuted. +c the upper triangular part of the +c symmetric matrix (a**t*a)**(-1) has +c replaced the upper triangular part of +c the a array. + write (6,210) + do 170 i=1,n + 170 write (6,220) (i,j,a(i,j),j=i,n) + 180 continue + stop + 190 format (1h0,8x,17hresidual length =,e12.4) + 200 format (1h0,8x,34hestimated parameters, x=a**(+)*b,,22h computed + 1by 'hfti' //(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8)) + 210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa + 1rameters.,19h computed by 'cov'./1x) + 220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8) + 230 format (52h1 prog2. this program demonstates the algorithms,16 + 1h hfti and cov.) + 240 format (1h0,54hthe relative noise level of the generated data will + 1 be,e16.4/33h0the matrix norm is approximately,e12.4/43h0the absol + 2ute pseudorank tolerance, tau, is,e12.4) + 250 format (1h0////9h0 m n/1x,2i4) + 260 format (1h0,8x,12hpseudorank =,i4) + end diff --git a/math/llsq/progs/prog3.f b/math/llsq/progs/prog3.f new file mode 100644 index 00000000..290f8161 --- /dev/null +++ b/math/llsq/progs/prog3.f @@ -0,0 +1,138 @@ +c prog3 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate the use of subroutine svdrs to compute the +c singular value decomposition of a matrix, a, and solve a least +c squares problem, a*x=b. +c +c the s.v.d. a= u*(s,0)**t*v**t is +c computed so that.. +c (1) u**t*b replaces the m+1 st. col. of b. +c +c (2) u**t replaces the m by m identity in +c the first m cols. of b. +c +c (3) v replaces the first n rows and cols. +c of a. +c +c (4) the diagonal entries of the s matrix +c replace the first n entries of the array s. +c +c the array s( ) must be dimensioned at least 3*n . +c + dimension a(8,8),b(8,9),s(24),x(8),aa(8,8) + real gen,anoise + double precision sm + data mda,mdb/8,8/ +c + do 150 noise=1,2 + anoise = 0. + rho = 1.e-3 + if(noise .eq. 1) go to 5 + anoise = 1.e-4 + rho = 10. * anoise + 5 continue + write(6,230) + write(6,240) anoise,rho +c initialize data generation function +c .. + dummy= gen(-1.) +c + do 150 mn1=1,6,5 + mn2=mn1+2 + do 150 m=mn1,mn2 + do 150 n=mn1,mn2 + write (6,160) m,n + do 20 i=1,m + do 10 j=1,m + 10 b(i,j)=0. + b(i,i)=1. + do 20 j=1,n + a(i,j)= gen(anoise) + 20 aa(i,j)=a(i,j) + do 30 i=1,m + 30 b(i,m+1)= gen(anoise) +c +c the arrays are now filled in.. +c compute the s.v.d. +c ****************************************************** + call svdrs (a,mda,m,n,b,mdb,m+1,s) +c ****************************************************** +c + write (6,170) + write (6,220) (i,s(i),i=1,n) + write (6,180) + write (6,220) (i,b(i,m+1),i=1,m) +c +c test for disparity of ratio of singular values. +c .. + krank=n + tau=rho*s(1) + do 40 i=1,n + if (s(i).le.tau) go to 50 + 40 continue + go to 55 + 50 krank=i-1 + 55 write(6,250) tau, krank +c compute solution vector assuming pseudorank is krank +c .. + 60 do 70 i=1,krank + 70 b(i,m+1)=b(i,m+1)/s(i) + do 90 i=1,n + sm=0.d0 + do 80 j=1,krank + 80 sm=sm+a(i,j)*dble(b(j,m+1)) + 90 x(i)=sm +c compute predicted norm of residual vector. +c .. + srsmsq=0. + if (krank.eq.m) go to 110 + kp1=krank+1 + do 100 i=kp1,m + 100 srsmsq=srsmsq+b(i,m+1)**2 + srsmsq=sqrt(srsmsq) +c + 110 continue + write (6,190) + write (6,220) (i,x(i),i=1,n) + write (6,200) srsmsq +c compute the frobenius norm of a**t- v*(s,0)*u**t. +c +c compute v*s first. +c + minmn=min0(m,n) + do 120 j=1,minmn + do 120 i=1,n + 120 a(i,j)=a(i,j)*s(j) + dn=0. + do 140 j=1,m + do 140 i=1,n + sm=0.d0 + do 130 l=1,minmn + 130 sm=sm+a(i,l)*dble(b(l,j)) +c computed difference of (i,j) th entry +c of a**t-v*(s,0)*u**t. +c .. + t=aa(j,i)-sm + 140 dn=dn+t**2 + dn=sqrt(dn)/(sqrt(float(n))*s(1)) + write (6,210) dn + 150 continue + stop + 160 format (1h0////9h0 m n/1x,2i4) + 170 format (1h0,8x,25hsingular values of matrix) + 180 format (1h0,8x,30htransformed right side, u**t*b) + 190 format (1h0,8x,33hestimated parameters, x=a**(+)*b,21h computed b + 1y 'svdrs' ) + 200 format (1h0,8x,24hresidual vector length =,e12.4) + 210 format (1h0,8x,43hfrobenius norm (a-u*(s,0)**t*v**t)/(sqrt(n), + *22h*spectral norm of a) =,e12.4) + 220 format (9x,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8) + 230 format(51h1prog3. this program demonstrates the algorithm, + * 9h, svdrs. ) + 240 format(55h0the relative noise level of the generated data will be, + *e16.4/44h0the relative tolerance, rho, for pseudorank, + *17h determination is,e16.4) + 250 format(1h0,8x,36habsolute pseudorank tolerance, tau =, + *e12.4,10x,12hpseudorank =,i5) + end diff --git a/math/llsq/progs/prog4.f b/math/llsq/progs/prog4.f new file mode 100644 index 00000000..090b1865 --- /dev/null +++ b/math/llsq/progs/prog4.f @@ -0,0 +1,22 @@ +c prog4 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate singular value analysis. +c + dimension a(15,5),b(15),sing(15) + data names/1h / +c + read (5,10) ((a(i,j),j=1,5),b(i),i=1,15) + write (6,20) + write (6,30) ((a(i,j),j=1,5),b(i),i=1,15) + write (6,40) +c + call sva (a,15,15,5,15,b,sing,names,1,d) +c + stop + 10 format (6f12.0) + 20 format (46h1prog4. demonstrate singular value analysis/53h list + 1ing of input matrix, a, and vector, b, follows..) + 30 format (1h /(5f12.8,f20.4)) + 40 format (1h1) + end diff --git a/math/llsq/progs/prog5.f b/math/llsq/progs/prog5.f new file mode 100644 index 00000000..5484b846 --- /dev/null +++ b/math/llsq/progs/prog5.f @@ -0,0 +1,146 @@ +c prog5 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c example of the use of subroutines bndacc and bndsol to solve +c sequentially the banded least squares problem which arises in +c spline curve fitting. +c + dimension x(12),y(12),b(10),g(12,5),c(12),q(4),cov(12) +c +c define functions p1 and p2 to be used in constructing +c cubic splines over uniformly spaced breakpoints. +c + p1(t)=.25*t**2*t + p2(t)=-(1.-t)**2*(1.+t)*.75+1. + zero=0. +c + write (6,230) + mdg=12 + nband=4 + m=12 +c set ordinates of data + y(1)=2.2 + y(2)=4.0 + y(3)=5.0 + y(4)=4.6 + y(5)=2.8 + y(6)=2.7 + y(7)=3.8 + y(8)=5.1 + y(9)=6.1 + y(10)=6.3 + y(11)=5.0 + y(12)=2.0 +c set abcissas of data + do 10 i=1,m + 10 x(i)=2*i +c +c begin loop thru cases using increasing nos of breakpoints. +c + do 150 nbp=5,10 + nc=nbp+2 +c set breakpoints + b(1)=x(1) + b(nbp)=x(m) + h=(b(nbp)-b(1))/float(nbp-1) + if (nbp.le.2) go to 30 + do 20 i=3,nbp + 20 b(i-1)=b(i-2)+h + 30 continue + write (6,240) nbp,(b(i),i=1,nbp) +c +c initialize ir and ip before first call to bndacc. +c + ir=1 + ip=1 + i=1 + jt=1 + 40 mt=0 + 50 continue + if (x(i).gt.b(jt+1)) go to 60 +c +c set row for ith data point +c + u=(x(i)-b(jt))/h + ig=ir+mt + g(ig,1)=p1(1.-u) + g(ig,2)=p2(1.-u) + g(ig,3)=p2(u) + g(ig,4)=p1(u) + g(ig,5)=y(i) + mt=mt+1 + if (i.eq.m) go to 60 + i=i+1 + go to 50 +c +c send block of data to processor +c + 60 continue + call bndacc (g,mdg,nband,ip,ir,mt,jt) + if (i.eq.m) go to 70 + jt=jt+1 + go to 40 +c compute solution c() + 70 continue + call bndsol (1,g,mdg,nband,ip,ir,c,nc,rnorm) +c write solution coefficients + write (6,180) (c(l),l=1,nc) + write (6,210) rnorm +c +c compute and print x,y,yfit,r=y-yfit +c + write (6,160) + jt=1 + do 110 i=1,m + 80 if (x(i).le.b(jt+1)) go to 90 + jt=jt+1 + go to 80 +c + 90 u=(x(i)-b(jt))/h + q(1)=p1(1.-u) + q(2)=p2(1.-u) + q(3)=p2(u) + q(4)=p1(u) + yfit=zero + do 100 l=1,4 + ic=jt-1+l + 100 yfit=yfit+c(ic)*q(l) + r=y(i)-yfit + write (6,170) i,x(i),y(i),yfit,r + 110 continue +c +c compute residual vector norm. +c + if (m.le.nc) go to 150 + sigsq=(rnorm**2)/float(m-nc) + sigfac=sqrt(sigsq) + write (6,220) sigfac + write (6,200) +c +c compute and print cols. of covariance. +c + do 140 j=1,nc + do 120 i=1,nc + 120 cov(i)=zero + cov(j)=1. + call bndsol (2,g,mdg,nband,ip,ir,cov,nc,rdummy) + call bndsol (3,g,mdg,nband,ip,ir,cov,nc,rdummy) +c +c compute the jth col. of the covariance matrix. + do 130 i=1,nc + 130 cov(i)=cov(i)*sigsq + 140 write (6,190) (l,j,cov(l),l=1,nc) + 150 continue + stop + 160 format (4h0 i,8x,1hx,10x,1hy,6x,4hyfit,4x,8hr=y-yfit/1x) + 170 format (1x,i3,4x,f6.0,4x,f6.2,4x,f6.2,4x,f8.4) + 180 format (4h0c =,10f10.5/(4x,10f10.5)) + 190 format (4(02x,2i5,e15.7)) + 200 format (46h0covariance matrix of the spline coefficients.) + 210 format (9h0rnorm =,e15.8) + 220 format (9h0sigfac =,e15.8) + 230 format (50h1prog5. execute a sequence of cubic spline fits,27h + 1to a discrete set of data.) + 240 format (1x////,11h0new case..,/47h0number of breakpoints, includin + 1g endpoints, is,i5/14h0breakpoints =,10f10.5,/(14x,10f10.5)) + end diff --git a/math/llsq/progs/prog6.f b/math/llsq/progs/prog6.f new file mode 100644 index 00000000..71364ec4 --- /dev/null +++ b/math/llsq/progs/prog6.f @@ -0,0 +1,116 @@ +c prog6 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 15 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate the use of the subroutine ldp for least +c distance programming by solving the constrained line data fitting +c problem of chapter 23. +c + dimension e(4,2), f(4), g(3,2), h(3), g2(3,2), h2(3), x(2), z(2), + 1w(4) + dimension wldp(21), s(6), t(4) + integer index(3) +c + write (6,110) + mde=4 + mdgh=3 +c + me=4 + mg=3 + n=2 +c define the least squares and constraint matrices. + t(1)=0.25 + t(2)=0.5 + t(3)=0.5 + t(4)=0.8 +c + w(1)=0.5 + w(2)=0.6 + w(3)=0.7 + w(4)=1.2 +c + do 10 i=1,me + e(i,1)=t(i) + e(i,2)=1. + 10 f(i)=w(i) +c + g(1,1)=1. + g(1,2)=0. + g(2,1)=0. + g(2,2)=1. + g(3,1)=-1. + g(3,2)=-1. +c + h(1)=0. + h(2)=0. + h(3)=-1. +c +c compute the singular value decomposition of the matrix, e. +c + call svdrs (e,mde,me,n,f,1,1,s) +c + write (6,120) ((e(i,j),j=1,n),i=1,n) + write (6,130) f,(s(j),j=1,n) +c +c generally rank determination and levenberg-marquardt +c stabilization could be inserted here. +c +c define the constraint matrix for the z coordinate system. + do 30 i=1,mg + do 30 j=1,n + sm=0. + do 20 l=1,n + 20 sm=sm+g(i,l)*e(l,j) + 30 g2(i,j)=sm/s(j) +c define constraint rt side for the z coordinate system. + do 50 i=1,mg + sm=0. + do 40 j=1,n + 40 sm=sm+g2(i,j)*f(j) + 50 h2(i)=h(i)-sm +c + write (6,140) ((g2(i,j),j=1,n),i=1,mg) + write (6,150) h2 +c +c solve the constrained problem in z-coordinates. +c + call ldp (g2,mdgh,mg,n,h2,z,znorm,wldp,index,mode) +c + write (6,200) mode,znorm + write (6,160) z +c +c transform back from z-coordinates to x-coordinates. + do 60 j=1,n + 60 z(j)=(z(j)+f(j))/s(j) + do 80 i=1,n + sm=0. + do 70 j=1,n + 70 sm=sm+e(i,j)*z(j) + 80 x(i)=sm + res=znorm**2 + np1=n+1 + do 90 i=np1,me + 90 res=res+f(i)**2 + res=sqrt(res) +c compute the residuals. + do 100 i=1,me + 100 f(i)=w(i)-x(1)*t(i)-x(2) + write (6,170) (x(j),j=1,n) + write (6,180) (i,f(i),i=1,me) + write (6,190) res + stop +c + 110 format (46h0prog6.. example of constrained curve fitting,26h usin + 1g the subroutine ldp.,/43h0related intermediate quantities are giv + 2en.) + 120 format (10h0v = ,2f10.5/(10x,2f10.5)) + 130 format (10h0f tilda =,4f10.5/10h0s = ,2f10.5) + 140 format (10h0g tilda =,2f10.5/(10x,2f10.5)) + 150 format (10h0h tilda =,3f10.5) + 160 format (10h0z = ,2f10.5) + 170 format (52h0the coeficients of the fitted line f(t)=x(1)*t+x(2),12 + 1h are x(1) = ,f10.5,14h and x(2) = ,f10.5) + 180 format (30h0the consecutive residuals are/1x,4(i10,f10.5)) + 190 format (23h0the residuals norm is ,f10.5) + 200 format (18h0mode (from ldp) =,i3,2x,7hznorm =,f10.5) +c + end diff --git a/math/llsq/qrbd.f b/math/llsq/qrbd.f new file mode 100644 index 00000000..d9a12a24 --- /dev/null +++ b/math/llsq/qrbd.f @@ -0,0 +1,208 @@ +c subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c qr algorithm for singular values of a bidiagonal matrix. +c +c the bidiagonal matrix +c +c (q1,e2,0... ) +c ( q2,e3,0... ) +c d= ( . ) +c ( . 0) +c ( .en) +c ( 0,qn) +c +c is pre and post multiplied by +c elementary rotation matrices +c ri and pi so that +c +c rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn) +c +c to within working accuracy. +c +c 1. ei and qi occupy e(i) and q(i) as input. +c +c 2. rm...r1*c replaces 'c' in storage as output. +c +c 3. v*p1**(t)...pm**(t) replaces 'v' in storage as output. +c +c 4. si occupies q(i) as output. +c +c 5. the si's are nonincreasing and nonnegative. +c +c this code is based on the paper and 'algol' code.. +c ref.. +c 1. reinsch,c.h. and golub,g.h. 'singular value decomposition +c and least squares solutions' (numer. math.), vol. 14,(1970). +c + subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) + logical wntv ,havers,fail + dimension q(nn),e(nn),v(mdv,nn),c(mdc,ncc) + zero=0. + one=1. + two=2. +c + n=nn + ipass=1 + if (n.le.0) return + n10=10*n + wntv=nrv.gt.0 + havers=ncc.gt.0 + fail=.false. + nqrs=0 + e(1)=zero + dnorm=zero + do 10 j=1,n + 10 dnorm=amax1(abs(q(j))+abs(e(j)),dnorm) + do 200 kk=1,n + k=n+1-kk +c +c test for splitting or rank deficiencies.. +c first make test for last diagonal term, q(k), being small. + 20 if(k.eq.1) go to 50 + if(diff(dnorm+q(k),dnorm)) 50,25,50 +c +c since q(k) is small we will make a special pass to +c transform e(k) to zero. +c + 25 cs=zero + sn=-one + do 40 ii=2,k + i=k+1-ii + f=-sn*e(i+1) + e(i+1)=cs*e(i+1) + call g1 (q(i),f,cs,sn,q(i)) +c transformation constructed to zero position (i,k). +c + if (.not.wntv) go to 40 + do 30 j=1,nrv + 30 call g2 (cs,sn,v(j,i),v(j,k)) +c accumulate rt. transformations in v. +c + 40 continue +c +c the matrix is now bidiagonal, and of lower order +c since e(k) .eq. zero.. +c + 50 do 60 ll=1,k + l=k+1-ll + if(diff(dnorm+e(l),dnorm)) 55,100,55 + 55 if(diff(dnorm+q(l-1),dnorm)) 60,70,60 + 60 continue +c this loop can't complete since e(1) = zero. +c + go to 100 +c +c cancellation of e(l), l.gt.1. + 70 cs=zero + sn=-one + do 90 i=l,k + f=-sn*e(i) + e(i)=cs*e(i) + if(diff(dnorm+f,dnorm)) 75,100,75 + 75 call g1 (q(i),f,cs,sn,q(i)) + if (.not.havers) go to 90 + do 80 j=1,ncc + 80 call g2 (cs,sn,c(i,j),c(l-1,j)) + 90 continue +c +c test for convergence.. + 100 z=q(k) + if (l.eq.k) go to 170 +c +c shift from bottom 2 by 2 minor of b**(t)*b. + x=q(l) + y=q(k-1) + g=e(k-1) + h=e(k) + f=((y-z)*(y+z)+(g-h)*(g+h))/(two*h*y) + g=sqrt(one+f**2) + if (f.lt.zero) go to 110 + t=f+g + go to 120 + 110 t=f-g + 120 f=((x-z)*(x+z)+h*(y/t-h))/x +c +c next qr sweep.. + cs=one + sn=one + lp1=l+1 + do 160 i=lp1,k + g=e(i) + y=q(i) + h=sn*g + g=cs*g + call g1 (f,h,cs,sn,e(i-1)) + f=x*cs+g*sn + g=-x*sn+g*cs + h=y*sn + y=y*cs + if (.not.wntv) go to 140 +c +c accumulate rotations (from the right) in 'v' + do 130 j=1,nrv + 130 call g2 (cs,sn,v(j,i-1),v(j,i)) + 140 call g1 (f,h,cs,sn,q(i-1)) + f=cs*g+sn*y + x=-sn*g+cs*y + if (.not.havers) go to 160 + do 150 j=1,ncc + 150 call g2 (cs,sn,c(i-1,j),c(i,j)) +c apply rotations from the left to +c right hand sides in 'c'.. +c + 160 continue + e(l)=zero + e(k)=f + q(k)=x + nqrs=nqrs+1 + if (nqrs.le.n10) go to 20 +c return to 'test for splitting'. +c + fail=.true. +c .. +c cutoff for convergence failure. 'nqrs' will be 2*n usually. + 170 if (z.ge.zero) go to 190 + q(k)=-z + if (.not.wntv) go to 190 + do 180 j=1,nrv + 180 v(j,k)=-v(j,k) + 190 continue +c convergence. q(k) is made nonnegative.. +c + 200 continue + if (n.eq.1) return + do 210 i=2,n + if (q(i).gt.q(i-1)) go to 220 + 210 continue + if (fail) ipass=2 + return +c .. +c every singular value is in order.. + 220 do 270 i=2,n + t=q(i-1) + k=i-1 + do 230 j=i,n + if (t.ge.q(j)) go to 230 + t=q(j) + k=j + 230 continue + if (k.eq.i-1) go to 270 + q(k)=q(i-1) + q(i-1)=t + if (.not.havers) go to 250 + do 240 j=1,ncc + t=c(i-1,j) + c(i-1,j)=c(k,j) + 240 c(k,j)=t + 250 if (.not.wntv) go to 270 + do 260 j=1,nrv + t=v(j,i-1) + v(j,i-1)=v(j,k) + 260 v(j,k)=t + 270 continue +c end of ordering algorithm. +c + if (fail) ipass=2 + return + end diff --git a/math/llsq/sva.f b/math/llsq/sva.f new file mode 100644 index 00000000..ee7f893d --- /dev/null +++ b/math/llsq/sva.f @@ -0,0 +1,198 @@ + subroutine sva (a,mda,m,n,mdata,b,sing,names,iscale,d) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c singular value analysis printout +c +c iscale set by user to 1, 2, or 3 to select column scaling +c option. +c 1 subr will use identity scaling and ignore the d() +c array. +c 2 subr will scale nonzero cols to have unit euclide +c length and will store reciprocal lengths of +c original nonzero cols in d(). +c 3 user supplies col scale factors in d(). subr +c will mult col j by d(j) and remove the scaling +c from the soln at the end. +c + dimension a(mda,n), b(m), sing(01),d(n) +c sing(3*n) + dimension names(n) + logical test + double precision sb, dzero + dzero=0.d0 + one=1. + zero=0. + if (m.le.0 .or. n.le.0) return + np1=n+1 + write (6,270) + write (6,260) m,n,mdata + go to (60,10,10), iscale +c +c apply column scaling to a +c + 10 do 50 j=1,n + a1=d(j) + go to (20,20,40), iscale + 20 sb=dzero + do 30 i=1,m + 30 sb=sb+dble(a(i,j))**2 + a1=dsqrt(sb) + if (a1.eq.zero) a1=one + a1=one/a1 + d(j)=a1 + 40 do 50 i=1,m + 50 a(i,j)=a(i,j)*a1 + write (6,280) iscale,(d(j),j=1,n) + go to 70 + 60 continue + write (6,290) + 70 continue +c +c obtain sing. value decomp. of scaled matrix. +c +c********************************************************** + call svdrs (a,mda,m,n,b,1,1,sing,ier) + if (ier.eq.0) goto 75 + write (6,295) + return + 75 continue +c********************************************************** +c +c print the v matrix. +c + call mfeout (a,mda,n,n,names,1) + if (iscale.eq.1) go to 90 +c +c replace v by d*v in the array a() +c + do 80 i=1,n + do 80 j=1,n + 80 a(i,j)=d(i)*a(i,j) +c +c g now in b array. v now in a array. +c +c +c obtain summary output. +c + 90 continue + write (6,220) +c +c compute cumulative sums of squares of components of +c g and store them in sing(i), i=minmn+1,...,2*minmn+1 +c + sb=dzero + minmn=min0(m,n) + minmn1=minmn+1 + if (m.eq.minmn) go to 110 + do 100 i=minmn1,m + 100 sb=sb+dble(b(i))**2 + 110 sing(2*minmn+1)=sb + do 120 jj=1,minmn + j=minmn+1-jj + sb=sb+dble(b(j))**2 + js=minmn+j + 120 sing(js)=sb + a3=sing(minmn+1) + a4=sqrt(a3/float(max0(1,mdata))) + write (6,230) a3,a4 +c + nsol=0 +c +c +c + do 160 k=1,minmn + if (sing(k).eq.zero) go to 130 + nsol=k + pi=b(k)/sing(k) + a1=one/sing(k) + a2=b(k)**2 + a3=sing(minmn1+k) + a4=sqrt(a3/float(max0(1,mdata-k))) + test=sing(k).ge.100..or.sing(k).lt..001 + if (test) write (6,240) k,sing(k),pi,a1,b(k),a2,a3,a4 + if (.not.test) write (6,250) k,sing(k),pi,a1,b(k),a2,a3,a4 + go to 140 + 130 write (6,240) k,sing(k) + pi=zero + 140 do 150 i=1,n + a(i,k)=a(i,k)*pi + 150 if (k.gt.1) a(i,k)=a(i,k)+a(i,k-1) + 160 continue +c +c compute and print values of ynorm and rnorm. +c + write (6,300) + j=0 + ysq=zero + go to 180 + 170 j=j+1 + ysq=ysq+(b(j)/sing(j))**2 + 180 ynorm=sqrt(ysq) + js=minmn1+j + rnorm=sqrt(sing(js)) + yl=-1000. + if (ynorm.gt.0.) yl=alog10(ynorm) + rl=-1000. + if (rnorm.gt.0.) rl=alog10(rnorm) + write (6,310) j,ynorm,rnorm,yl,rl + if (j.lt.nsol) go to 170 +c +c compute values of xnorm and rnorm for a sequence of values of +c the levenberg-marquardt parameter. +c + if (sing(1).eq.zero) go to 210 + el=alog10(sing(1))+one + el2=alog10(sing(nsol))-one + del=(el2-el)/20. + ten=10. + aln10=alog(ten) + write (6,320) + do 200 ie=1,21 +c compute alamb=10.**el + alamb=exp(aln10*el) + ys=0. + js=minmn1+nsol + rs=sing(js) + do 190 i=1,minmn + sl=sing(i)**2+alamb**2 + ys=ys+(b(i)*sing(i)/sl)**2 + rs=rs+(b(i)*(alamb**2)/sl)**2 + 190 continue + ynorm=sqrt(ys) + rnorm=sqrt(rs) + rl=-1000. + if (rnorm.gt.zero) rl=alog10(rnorm) + yl=-1000. + if (ynorm.gt.zero) yl=alog10(ynorm) + write (6,330) alamb,ynorm,rnorm,el,yl,rl + el=el+del + 200 continue +c +c print candidate solutions. +c + 210 if (nsol.ge.1) call mfeout (a,mda,n,nsol,names,2) + return + 220 format (42h0 index sing. value p coef ,48hrecip. s. + 1v. g coef g**2 , + 234h c.s.s. n.s.r.c.s.s.) + 230 format (1h ,5x,1h0,88x,1pe15.4,1pe17.4) + 240 format (1h ,i6,e12.4,1p,(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4 + 1,2x,e15.4)) + 250 format (1h ,i6,f12.4,1p,(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4 + 1,2x,e15.4)) + 260 format (5h0m = ,i6,8h, n = ,i4,12h, mdata = ,i8) + 270 format (45h0singular value analysis of the least squares,42h probl + 1em, a*x=b, scaled as (a*d)*y=b .) + 280 format (19h0scaling option no.,i2,18h. d is a diagonal,46h matrix + 1 with the following diagonal elements../(5x,10e12.4)) + 290 format (50h0scaling option no. 1. d is the identity matrix./1x) + 295 format (49h convergence failure in qr bidiagonal svd routine) + 300 format (6h0index,13x,28h ynorm rnorm,14x,28h log1 + 10(ynorm) log10(rnorm)/1x) + 310 format (1h ,i4,14x,2e14.5,14x,2f14.5) + 320 format (54h0norms of solution and residual vectors for a range of, + 144h values of the levenberg-marquardt parameter,9h, lambda./1h0,4x + 2,42h lambda ynorm rnorm, + 345h log10(lambda) log10(ynorm) log10(rnorm)) + 330 format (5x,3e14.5,3f14.5) + end diff --git a/math/llsq/svdrs.f b/math/llsq/svdrs.f new file mode 100644 index 00000000..46dba02d --- /dev/null +++ b/math/llsq/svdrs.f @@ -0,0 +1,205 @@ +c subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s,ier) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c singular value decomposition also treating right side vector. +c +c the array s occupies 3*n cells. +c a occupies m*n cells +c b occupies m*nb cells. +c +c special singular value decomposition subroutine. +c we have the m x n matrix a and the system a*x=b to solve. +c either m .ge. n or m .lt. n is permitted. +c the singular value decomposition +c a = u*s*v**(t) is made in such a way that one gets +c (1) the matrix v in the first n rows and columns of a. +c (2) the diagonal matrix of ordered singular values in +c the first n cells of the array s(ip), ip .ge. 3*n. +c (3) the matrix product u**(t)*b=g gets placed back in b. +c (4) the user must complete the solution and do his own +c singular value analysis. +c ******* +c give special +c treatment to rows and columns which are entirely zero. this +c causes certain zero sing. vals. to appear as exact zeros rather +c than as about eta times the largest sing. val. it similarly +c cleans up the associated columns of u and v. +c method.. +c 1. exchange cols of a to pack nonzero cols to the left. +c set n = no. of nonzero cols. +c use locations a(1,nn),a(1,nn-1),...,a(1,n+1) to record the +c col permutations. +c 2. exchange rows of a to pack nonzero rows to the top. +c quit packing if find n nonzero rows. make same row exchanges +c in b. set m so that all nonzero rows of the permuted a +c are in first m rows. if m .le. n then all m rows are +c nonzero. if m .gt. n then the first n rows are known +c to be nonzero,and rows n+1 thru m may be zero or nonzero. +c 3. apply original algorithm to the m by n problem. +c 4. move permutation record from a(,) to s(i),i=n+1,...,nn. +c 5. build v up from n by n to nn by nn by placing ones on +c the diagonal and zeros elsewhere. this is only partly done +c explicitly. it is completed during step 6. +c 6. exchange rows of v to compensate for col exchanges of step 2. +c 7. place zeros in s(i),i=n+1,nn to represent zero sing vals. +c + subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s,ier) + dimension a(mda,nn),b(mdb,nb),s(nn,3) + zero=0. + ier = 0 + one=1. +c +c begin.. special for zero rows and cols. +c +c pack the nonzero cols to the left +c + n=nn + if (n.le.0.or.mm.le.0) return + j=n + 10 continue + do 20 i=1,mm + if (a(i,j)) 50,20,50 + 20 continue +c +c col j is zero. exchange it with col n. +c + if (j.eq.n) go to 40 + do 30 i=1,mm + 30 a(i,j)=a(i,n) + 40 continue + a(1,n)=j + n=n-1 + 50 continue + j=j-1 + if (j.ge.1) go to 10 +c if n=0 then a is entirely zero and svd +c computation can be skipped + ns=0 + if (n.eq.0) go to 240 +c pack nonzero rows to the top +c quit packing if find n nonzero rows + i=1 + m=mm + 60 if (i.gt.n.or.i.ge.m) go to 150 + if (a(i,i)) 90,70,90 + 70 do 80 j=1,n + if (a(i,j)) 90,80,90 + 80 continue + go to 100 + 90 i=i+1 + go to 60 +c row i is zero +c exchange rows i and m + 100 if(nb.le.0) go to 115 + do 110 j=1,nb + t=b(i,j) + b(i,j)=b(m,j) + 110 b(m,j)=t + 115 do 120 j=1,n + 120 a(i,j)=a(m,j) + if (m.gt.n) go to 140 + do 130 j=1,n + 130 a(m,j)=zero + 140 continue +c exchange is finished + m=m-1 + go to 60 +c + 150 continue +c end.. special for zero rows and columns +c begin.. svd algorithm +c method.. +c (1) reduce the matrix to upper bidiagonal form with +c householder transformations. +c h(n)...h(1)aq(1)...q(n-2) = (d**t,0)**t +c where d is upper bidiagonal. +c +c (2) apply h(n)...h(1) to b. here h(n)...h(1)*b replaces b +c in storage. +c +c (3) the matrix product w= q(1)...q(n-2) overwrites the first +c n rows of a in storage. +c +c (4) an svd for d is computed. here k rotations ri and pi are +c computed so that +c rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sm) +c to working accuracy. the si are nonnegative and nonincreasing. +c here rk...r1*b overwrites b in storage while +c a*p1**(t)...pk**(t) overwrites a in storage. +c +c (5) it follows that,with the proper definitions, +c u**(t)*b overwrites b, while v overwrites the first n row and +c columns of a. +c + l=min0(m,n) +c the following loop reduces a to upper bidiagonal and +c also applies the premultiplying transformations to b. +c + do 170 j=1,l + if (j.ge.m) go to 160 + call h12 (1,j,j+1,m,a(1,j),1,t,a(1,j+1),1,mda,n-j) + call h12 (2,j,j+1,m,a(1,j),1,t,b,1,mdb,nb) + 160 if (j.ge.n-1) go to 170 + call h12 (1,j+1,j+2,n,a(j,1),mda,s(j,3),a(j+1,1),mda,1,m-j) + 170 continue +c +c copy the bidiagonal matrix into the array s() for qrbd. +c + if (n.eq.1) go to 190 + do 180 j=2,n + s(j,1)=a(j,j) + 180 s(j,2)=a(j-1,j) + 190 s(1,1)=a(1,1) +c + ns=n + if (m.ge.n) go to 200 + ns=m+1 + s(ns,1)=zero + s(ns,2)=a(m,m+1) + 200 continue +c +c construct the explicit n by n product matrix, w=q1*q2*...*ql*i +c in the array a(). +c + do 230 k=1,n + i=n+1-k + if(i.gt.min0(m,n-2)) go to 210 + call h12 (2,i+1,i+2,n,a(i,1),mda,s(i,3),a(1,i+1),1,mda,n-i) + 210 do 220 j=1,n + 220 a(i,j)=zero + 230 a(i,i)=one +c +c compute the svd of the bidiagonal matrix +c + call qrbd (ipass,s(1,1),s(1,2),ns,a,mda,n,b,mdb,nb) +c + go to (240,310), ipass + 240 continue + if (ns.ge.n) go to 260 + nsp1=ns+1 + do 250 j=nsp1,n + 250 s(j,1)=zero + 260 continue + if (n.eq.nn) return + np1=n+1 +c move record of permutations +c and store zeros + do 280 j=np1,nn + s(j,1)=a(1,j) + do 270 i=1,n + 270 a(i,j)=zero + 280 continue +c permute rows and set zero singular values. + do 300 k=np1,nn + i=s(k,1) + s(k,1)=zero + do 290 j=1,nn + a(k,j)=a(i,j) + 290 a(i,j)=zero + a(i,k)=one + 300 continue +c end.. special for zero rows and columns + return +c convergence failure in qr bidiagonal svd routine + 310 ier = 1 + end -- cgit