From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 math/llsq/progs/band.x | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 70 insertions(+)
 create mode 100644 math/llsq/progs/band.x

(limited to 'math/llsq/progs/band.x')

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
-- 
cgit