Initial commit

1 files changed, 70 insertions, 0 deletions

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