Initial commit

1 files changed, 125 insertions, 0 deletions

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