Initial commit

1 files changed, 124 insertions, 0 deletions

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