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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <mach.h>
include <math/gsurfit.h>
include "dgsurfitdef.h"
# GSCHOFAC -- 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 dgschofac (matrix, nbands, nrows, matfac, ier)
double matrix[nbands, nrows] # data matrix
int nbands # number of bands
int nrows # number of rows
double matfac[nbands, nrows] # Cholesky factorization
int ier # error code
int i, n, j, imax, jmax
double ratio
begin
if (nrows == 1) {
if (matrix[1,1] > 0.)
matfac[1,1] = 1. / matrix[1,1]
return
}
# copy matrix into matfac
do n = 1, nrows {
do j = 1, nbands
matfac[j,n] = matrix[j,n]
}
do n = 1, nrows {
# test to see if matrix is singular
if (((matfac[1,n]+matrix[1,n])-matrix[1,n]) <= 1000/MAX_DOUBLE) {
do j = 1, nbands
matfac[j,n] = 0.
ier = SINGULAR
next
}
matfac[1,n] = 1. / 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
# GSCHOSLV -- Solve the matrix whose Cholesky factorization was calculated in
# GSCHOFAC for the coefficients. This routine was adapted from bchslv.f
# described in "A Practical Guide to Splines", by Carl de Boor (1978).
procedure dgschoslv (matfac, nbands, nrows, vector, coeff)
double matfac[nbands,nrows] # Cholesky factorization
int nbands # number of bands
int nrows # number of rows
double vector[nrows] # right side of matrix equation
double coeff[nrows] # coefficients
int i, n, j, jmax, nbndm1
begin
if (nrows == 1) {
coeff[1] = vector[1] * matfac[1,1]
return
}
# copy vector to coefficients
do i = 1, nrows
coeff[i] = vector[i]
# 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]
}
}
# back 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
|
