Initial commit

1 files changed, 138 insertions, 0 deletions

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