23 lines
446 B
Text
23 lines
446 B
Text
void ludec(real matrix a, real matrix l, real matrix u, real vector p) {
|
|
real scalar i,j,n,s
|
|
real vector js
|
|
|
|
l = a
|
|
n = rows(a)
|
|
p = 1::n
|
|
for (i=1; i<n; i++) {
|
|
maxindex(abs(l[i::n,i]), 1, js=., .)
|
|
j = js[1]+i-1
|
|
if (j!=i) {
|
|
l[(i\j),.] = l[(j\i),.]
|
|
p[(i\j)] = p[(j\i)]
|
|
}
|
|
for (j=i+1; j<=n; j++) {
|
|
l[j,i] = s = l[j,i]/l[i,i]
|
|
l[j,i+1..n] = l[j,i+1..n]-s*l[i,i+1..n]
|
|
}
|
|
}
|
|
|
|
u = uppertriangle(l)
|
|
l = lowertriangle(l, 1)
|
|
}
|