Repatch (from linux) of OSX IRAF

1 files changed, 116 insertions, 0 deletions

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