1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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
|
