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)
}