aboutsummaryrefslogtreecommitdiff
path: root/math/llsq/progs
diff options
context:
space:
mode:
Diffstat (limited to 'math/llsq/progs')
-rw-r--r--math/llsq/progs/README5
-rw-r--r--math/llsq/progs/band.x70
-rw-r--r--math/llsq/progs/data415
-rw-r--r--math/llsq/progs/lsq.x70
-rw-r--r--math/llsq/progs/prog1.f124
-rw-r--r--math/llsq/progs/prog2.f125
-rw-r--r--math/llsq/progs/prog3.f138
-rw-r--r--math/llsq/progs/prog4.f22
-rw-r--r--math/llsq/progs/prog5.f146
-rw-r--r--math/llsq/progs/prog6.f116
10 files changed, 831 insertions, 0 deletions
diff --git a/math/llsq/progs/README b/math/llsq/progs/README
new file mode 100644
index 00000000..e0c5a142
--- /dev/null
+++ b/math/llsq/progs/README
@@ -0,0 +1,5 @@
+
+This directory contains both Fortran and IRAF preprocessor programs
+demonstrating the use of the Lawson-Hanson llsq procedures. The Fortran
+programs have not been modified to reflect the removal of Fortran i/o
+from the library procedures as installed in the IRAF llsq.a library.
diff --git a/math/llsq/progs/band.x b/math/llsq/progs/band.x
new file mode 100644
index 00000000..f65290ba
--- /dev/null
+++ b/math/llsq/progs/band.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+define MAXPTS 1024
+
+task band
+
+# This procedure solves for the natural cubic spline which interpolates
+# the curve y[i] = 100., for i = 1 to n. The matrix is of length n+2,
+# and has a bandwidth of 3. Lawsons and Hansons routines BNDACC and
+# BNDSOL (modified to return an error code IER) are used to solve the
+# system. The computations are carried out in single precision.
+
+
+procedure band()
+
+real g[MAXPTS,4], c[MAXPTS], rnorm
+int ip, ir, i, n, ier, geti()
+real marktime, cptime()
+
+begin
+ n = min (MAXPTS, geti ("npts"))
+ marktime = cptime()
+
+ ip = 1
+ ir = 1
+
+ g[ir,1] = 6. # first row
+ g[ir,2] = -12.
+ g[ir,3] = 6.
+ g[ir,4] = 0.
+ call bndacc (g, MAXPTS, 3, ip, ir, 1, 1)
+
+ do i = 2, n-1 { # tridiagonal elements
+ g[ir,1] = 1.
+ g[ir,2] = 4.
+ g[ir,3] = 1.
+ g[ir,4] = 100.
+ call bndacc (g, MAXPTS, 3, ip, ir, 1, i-1)
+ }
+
+ g[ir,1] = 6. # last row
+ g[ir,2] = -12.
+ g[ir,3] = 6.
+ g[ir,4] = 0.
+ call bndacc (g, MAXPTS, 3, ip, ir, 1, n-2)
+
+ call printf ("matrix accumulation took %8.2f cpu seconds\n")
+ call pargr (cptime() - marktime)
+ marktime = cptime()
+
+ # solve for the b-spline coeff C
+ call bndsol (1, g, MAXPTS, 3, ip, ir, c, n, rnorm, ier)
+
+ call printf ("solution took %8.2f cpu seconds (rnorm=%g, ier=%d)\n")
+ call pargr (cptime() - marktime)
+ call pargr (rnorm)
+ call pargi (ier)
+
+ call printf ("selected coefficients:\n")
+ for (i=1; i <= n;) { # print the first and last 4 coeff
+ call printf ("%8d%15.6f\n")
+ call pargi (i)
+ call pargr (c[i])
+ if (i == 4)
+ i = max(i+1, n-3)
+ else
+ i = i + 1
+ }
+end
diff --git a/math/llsq/progs/data4 b/math/llsq/progs/data4
new file mode 100644
index 00000000..33bf1b08
--- /dev/null
+++ b/math/llsq/progs/data4
@@ -0,0 +1,15 @@
+ -.13405547 -.20162827 -.16930778 -.18971990 -.17387234 -.4361
+ -.10379475 -.15766336 -.13346256 -.14848550 -.13597690 -.3437
+ -.08779597 -.12883867 -.10683007 -.12011796 -.10932972 -.2657
+ .02058554 .00335331 -.01641270 .00078606 .00271659 -.0392
+ -.03248093 -.01876799 .00410639 -.01405894 -.01384391 .0193
+ .05967662 .06667714 .04352153 .05740438 .05024962 .0747
+ .06712457 .07352437 .04489770 .06471862 .05876455 .0935
+ .08687186 .09368296 .05672327 .08141043 .07302320 .1079
+ .02149662 .06222662 .07213486 .06200069 .05570931 .1930
+ .06687407 .10344506 .09153849 .09508223 .08393667 .2058
+ .15879069 .18088339 .11540692 .16160727 .14796479 .2606
+ .17642887 .20361830 .13057860 .18385729 .17005549 .3142
+ .11414080 .17259611 .14816471 .16007466 .14374096 .3529
+ .07846038 .14669563 .14365800 .14003842 .12571177 .3615
+ .10803175 .16994623 .14971519 .15885312 .14301547 .3647
diff --git a/math/llsq/progs/lsq.x b/math/llsq/progs/lsq.x
new file mode 100644
index 00000000..ae002a16
--- /dev/null
+++ b/math/llsq/progs/lsq.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+define M 256
+define N 256
+
+task lsq
+
+# This procedure fits a natural cubic spline to an array of n data points.
+# The system being solved is a tridiagonal matrix of n+2 rows. The system
+# is solved by Lawsons and Hansons routine HFTI, which solves a general
+# m by n linear system of equations. This is enormous overkill for this
+# problem (see "band.x"), but serves to give timing estimates for the code.
+
+
+procedure lsq()
+
+real a[M,N], b[M], tau, rnorm, h[N], g[N]
+int krank, ip[N]
+int i, j, m, n, geti()
+real marktime, cptime()
+
+begin
+ m = min (M, geti ("npts")) # size of matrix
+ n = min (N, m)
+ tau = 1e-6
+
+ do j = 1, n # set up b-spline matrix
+ do i = 1, m
+ a[i,j] = 0.
+
+ a[1,1] = 6. # first row
+ a[1,2] = -12.
+ a[1,3] = 6.
+
+ a[m,n] = 6. # last row
+ a[m,n-1] = -12.
+ a[m,n-2] = 6.
+
+ do j = 2, m-1 { # tridiagonal elements
+ a[j,j-1] = 1.
+ a[j,j] = 4.
+ a[j,j+1] = 1.
+ }
+
+ b[1] = 0. # natural spline bndry conditions
+ b[m] = 0.
+
+ do i = 2, m-1 # set up data vector
+ b[i] = 100.
+
+ marktime = cptime()
+ call hfti (a,M,m,n, b,1,1, tau, krank,rnorm, h,g,ip)
+
+ call printf ("took %8.2f cpu seconds (krank=%d, rnorm=%g)\n")
+ call pargr (cptime() - marktime)
+ call pargi (krank)
+ call pargr (rnorm)
+
+ call printf ("selected coefficients:\n")
+ for (i=1; i <= m;) { # print first, last 4 coeff
+ call printf ("%8d%15.5f\n")
+ call pargi (i)
+ call pargr (b[i])
+ if (i == 4)
+ i = max(i+1, m-3)
+ else
+ i = i + 1
+ }
+end
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
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
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
diff --git a/math/llsq/progs/prog4.f b/math/llsq/progs/prog4.f
new file mode 100644
index 00000000..090b1865
--- /dev/null
+++ b/math/llsq/progs/prog4.f
@@ -0,0 +1,22 @@
+c prog4
+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 singular value analysis.
+c
+ dimension a(15,5),b(15),sing(15)
+ data names/1h /
+c
+ read (5,10) ((a(i,j),j=1,5),b(i),i=1,15)
+ write (6,20)
+ write (6,30) ((a(i,j),j=1,5),b(i),i=1,15)
+ write (6,40)
+c
+ call sva (a,15,15,5,15,b,sing,names,1,d)
+c
+ stop
+ 10 format (6f12.0)
+ 20 format (46h1prog4. demonstrate singular value analysis/53h list
+ 1ing of input matrix, a, and vector, b, follows..)
+ 30 format (1h /(5f12.8,f20.4))
+ 40 format (1h1)
+ end
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
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
generated by cgit (git 2.34.1) at 2025-09-20 10:28:37 -0400