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
|
define TINY 1.e-20
# ludcmp -- lower-upper decomposition
# This routine decomposes a matrix (in-place) into lower and upper
# triangular portions. This is the same as the Numerical Recipes version
# except that memory is allocated instead of using the fixed array VV.
#
# Phil Hodge, 28-Oct-1988 Subroutine copied from Numerical Recipes.
procedure ludcmp (a, n, np, indx, d)
real a[np,np] # io: input a, output decomposed a
int n # i: logical size of a is n x n
int np # i: space allocated for a
int indx[n] # o: index to be used by xt_lubksb
real d # o: +1 or -1
#--
pointer sp
pointer vv # scratch space
real aamax
real sum
real dum
int i, j, k
int imax
begin
call smark (sp)
call salloc (vv, n, TY_REAL)
d = 1.
do i = 1, n {
aamax = 0.
do j = 1, n
if (abs(a[i,j]) > aamax)
aamax = abs(a[i,j])
if (aamax == 0.)
call error (0, "singular matrix")
Memr[vv+i-1] = 1. / aamax
}
do j = 1, n {
if (j > 1) {
do i = 1, j-1 {
sum = a[i,j]
if (i > 1) {
do k = 1, i-1
sum = sum - a[i,k] * a[k,j]
a[i,j] = sum
}
}
}
aamax = 0.
do i = j, n {
sum = a[i,j]
if (j > 1) {
do k = 1, j-1
sum = sum - a[i,k] * a[k,j]
a[i,j] = sum
}
dum = Memr[vv+i-1] * abs[sum]
if (dum >= aamax) {
imax = i
aamax = dum
}
}
if (j != imax) {
do k = 1, n {
dum = a[imax,k]
a[imax,k] = a[j,k]
a[j,k] = dum
}
d = -d
Memr[vv+imax-1] = Memr[vv+j-1]
}
indx[j] = imax
if (j != n) {
if (a[j,j] == 0.)
a[j,j] = TINY
dum = 1. / a[j,j]
do i = j+1, n
a[i,j] = a[i,j] * dum
}
}
if (a[n,n] == 0.)
a[n,n] = TINY
call sfree (sp)
end
|
