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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/surfit.h>
include "surfitdef.h"
# ISLACCUM -- Procedure to add points on a line to the data set.
# The inner products of the non-zero x basis functions are stored
# in the SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The
# main diagonal is stored in the first row of XMATRIX. Successive
# non-zero diagonals are stored in the succeeding rows. This method
# of storage is particularly efficient for the large symmetric
# banded matrices produced during spline fits. The inner products
# of the data ordinates and the non-zero x basis functions are
# caculated and stored in the lineno-th row of the SF_NXCOEFF(sf)
# by SF_NLINES(sf) matrix XCOEFF.
procedure islaccum (sf, cols, lineno, z, w, ncols, wtflag)
pointer sf # pointer to surface descriptor
int cols[ncols] # column values
int lineno # lineno of data being accumulated
real z[ncols] # surface values on lineno at cols
real w[ncols] # weight of the data points
int ncols # number of data points
int wtflag # type of weighting desired
int i, ii, j, k
pointer xbzptr, xbptr
pointer xlzptr
pointer xmzptr, xmindex
pointer xczptr, xcindex
pointer bw, rows, left
pointer sp
begin
# count the number of points
SF_NXPTS(sf) = SF_NXPTS(sf) + ncols
# calculate the weights, default is uniform weighting
switch (wtflag) {
case SF_UNIFORM:
call amovkr (1.0, w, ncols)
case SF_USER:
# do not alter weights
default:
call amovkr (1.0, w, ncols)
}
# set up temporary storage
call smark (sp)
call salloc (bw, ncols, TY_REAL)
call salloc (left, ncols, TY_INT)
call salloc (rows, ncols, TY_INT)
# set up the pointers
xbzptr = SF_XBASIS(sf) - 1
xmzptr = SF_XMATRIX(sf)
xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
# accumulate the line
switch (SF_TYPE(sf)) {
case SF_LEGENDRE, SF_CHEBYSHEV:
do i = 1, SF_XORDER(sf) {
do j = 1, ncols
Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
xcindex = xczptr + i
do j = 1, ncols
XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
xbptr = xbzptr
ii = 0
do k = i, SF_XORDER(sf) {
xmindex = xmzptr + ii
do j = 1, ncols
XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
XBASIS(xbptr+cols[j])
ii = ii + 1
xbptr = xbptr + SF_NCOLS(sf)
}
xbzptr = xbzptr + SF_NCOLS(sf)
xmzptr = xmzptr + SF_XORDER(sf)
}
case SF_SPLINE3, SF_SPLINE1:
xlzptr = SF_XLEFT(sf) - 1
do j = 1, ncols
Memi[left+j-1] = XLEFT(xlzptr+cols[j])
call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
call aaddki (Memi[left], xczptr, Memi[left], ncols)
do i = 1, SF_XORDER(sf) {
do j = 1, ncols {
Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
xcindex = Memi[left+j-1] + i
XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
}
xbptr = xbzptr
ii = 0
do k = i, SF_XORDER(sf) {
do j = 1, ncols {
xmindex = Memi[rows+j-1] + ii
XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
XBASIS(xbptr+cols[j])
}
ii = ii + 1
xbptr = xbptr + SF_NCOLS(sf)
}
xbzptr = xbzptr + SF_NCOLS(sf)
call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
}
}
# release space
call sfree (sp)
end
|
