Initial commit

1 files changed, 146 insertions, 0 deletions

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