50 lines
1.1 KiB
Text
50 lines
1.1 KiB
Text
fcn make_array(n,m,v){ (m).pump(List.createLong(m).write,v)*n }
|
|
fcn eye(n){ // Creates a nxn identity matrix.
|
|
I:=make_array(n,n,0.0);
|
|
foreach j in (n){ I[j][j]=1.0 }
|
|
I
|
|
}
|
|
|
|
// Creates the pivoting matrix for A.
|
|
fcn pivotize(A){
|
|
n:=A.len(); // rows
|
|
P:=eye(n);
|
|
foreach i in (n){
|
|
max,row:=A[i][i],i;
|
|
foreach j in ([i..n-1]){
|
|
if(A[j][i]>max) max,row=A[j][i],j;
|
|
}
|
|
if(i!=row) P.swap(i,row);
|
|
}
|
|
// Return P.
|
|
P
|
|
}
|
|
|
|
// Decomposes a square matrix A by PA=LU and returns L, U and P.
|
|
fcn lu(A){
|
|
n:=A.len();
|
|
L:=eye(n);
|
|
U:=make_array(n,n,0.0);
|
|
P:=pivotize(A);
|
|
A=matMult(P,A);
|
|
|
|
foreach j in (n){
|
|
foreach i in (j+1){
|
|
U[i][j]=A[i][j] - (i).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0);
|
|
}
|
|
foreach i in ([j..n-1]){
|
|
L[i][j]=( A[i][j] -
|
|
(j).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0) ) /
|
|
U[j][j];
|
|
}
|
|
}
|
|
// Return L, U and P.
|
|
return(L,U,P);
|
|
}
|
|
|
|
fcn matMult(a,b){
|
|
n,m,p:=a[0].len(),a.len(),b[0].len();
|
|
ans:=make_array(n,m,0.0);
|
|
foreach i,j,k in (m,p,n){ ans[i][j]+=a[i][k]*b[k][j]; }
|
|
ans
|
|
}
|