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
118
119
120
121
122
123
124
125
126
127
128
129
130
|
include <mach.h>
include <math/nlfit.h>
$if (datatype == r)
include "nlfitdefr.h"
$else
include "nlfitdefd.h"
$endif
# NL_CHFAC -- Routine to calculate the Cholesky factorization of a
# symmetric, positive semi-definite banded matrix. This routines was
# adapted from the bchfac.f routine described in "A Practical Guide
# to Splines", Carl de Boor (1978).
procedure nl_chfac$t (matrix, nbands, nrows, matfac, ier)
PIXEL matrix[nbands, nrows] # data matrix
int nbands # number of bands
int nrows # number of rows
PIXEL matfac[nbands, nrows] # Cholesky factorization
int ier # error code
int i, n, j, imax, jmax
PIXEL ratio
begin
# Test for a single element matrix.
if (nrows == 1) {
if (matrix[1,1] > PIXEL (0.0))
matfac[1,1] = 1. / matrix[1,1]
return
}
# Copy the original matrix into matfac.
do n = 1, nrows {
do j = 1, nbands
matfac[j,n] = matrix[j,n]
}
# Compute the factorization of the matrix.
do n = 1, nrows {
# Test to see if matrix is singular.
if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= EPSILON$T) {
#if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= PIXEL(0.0)) {
do j = 1, nbands
matfac[j,n] = PIXEL (0.0)
ier = SINGULAR
next
}
matfac[1,n] = PIXEL (1.0) / matfac[1,n]
imax = min (nbands - 1, nrows - n)
if (imax < 1)
next
jmax = imax
do i = 1, imax {
ratio = matfac[i+1,n] * matfac[1,n]
do j = 1, jmax
matfac[j,n+i] = matfac[j,n+i] - matfac[j+i,n] * ratio
jmax = jmax - 1
matfac[i+1,n] = ratio
}
}
end
# NL_CHSLV -- Solve the matrix whose Cholesky factorization was calculated in
# NL_CHFAC for the coefficients. This routine was adapted from bchslv.f
# described in "A Practical Guide to Splines", by Carl de Boor (1978).
procedure nl_chslv$t (matfac, nbands, nrows, vector, coeff)
PIXEL matfac[nbands,nrows] # Cholesky factorization
int nbands # number of bands
int nrows # number of rows
PIXEL vector[nrows] # right side of matrix equation
PIXEL coeff[nrows] # coefficients
int i, n, j, jmax, nbndm1
begin
# Test for a single element matrix.
if (nrows == 1) {
coeff[1] = vector[1] * matfac[1,1]
return
}
# Copy input vector to coefficients vector.
do i = 1, nrows
coeff[i] = vector[i]
# Perform forward substitution.
nbndm1 = nbands - 1
do n = 1, nrows {
jmax = min (nbndm1, nrows - n)
if (jmax >= 1) {
do j = 1, jmax
coeff[j+n] = coeff[j+n] - matfac[j+1,n] * coeff[n]
}
}
# Perform backward substitution.
for (n = nrows; n >= 1; n = n - 1) {
coeff[n] = coeff[n] * matfac[1,n]
jmax = min (nbndm1, nrows - n)
if (jmax >= 1) {
do j = 1, jmax
coeff[n] = coeff[n] - matfac[j+1,n] * coeff[j+n]
}
}
end
# NL_DAMP -- Procedure to add damping to matrix
procedure nl_damp$t (inmatrix, outmatrix, constant, nbands, nrows)
PIXEL inmatrix[nbands,ARB] # input matrix
PIXEL outmatrix[nbands,ARB] # output matrix
PIXEL constant # damping constant
int nbands, nrows # dimensions of matrix
int i
begin
do i = 1, nrows
outmatrix[1,i] = inmatrix[1,i] * constant
end
|
