Skip to content

Instantly share code, notes, and snippets.

@ashwinvis
Forked from jfpuget/LU_decomposition.ipynb
Created March 9, 2020 12:29
Show Gist options
  • Save ashwinvis/473166210a79cf4fdfd53928d96fc319 to your computer and use it in GitHub Desktop.
Save ashwinvis/473166210a79cf4fdfd53928d96fc319 to your computer and use it in GitHub Desktop.
A Speed Comparison Of C, Julia, Python, Numba, and Cython on LU Factorization
// lu.c
inline int _(int row, int col, int rows){
return row*rows + col;
}
void det_by_lu(double *y, double *x, int N){
int i,j,k;
*y = 1.;
for(k = 0; k < N; ++k){
*y *= x[_(k,k,N)];
for(i = k+1; i < N; ++i){
x[_(i,k,N)] /= x[_(k,k,N)];
}
for(i = k+1; i < N; ++i){
#pragma omp simd
for(j = k+1; j < N; ++j){
x[_(i,j,N)] -= x[_(i,k,N)] * x[_(k,j,N)];
}
}
}
}
function det_by_lu(y, x, N)
y[1] = 1.
for k = 1:N
y[1] *= x[k,k]
for i = k+1:N
x[i,k] /= x[k,k]
end
for j = k+1:N
for i = k+1:N
x[i,j] -= x[i,k] * x[k,j]
end
end
end
end
function run_julia(y,A,B,N)
loops = max(10000000 // (N*N), 1)
print(loops)
for l in 1:loops
B[:,:] = A
det_by_lu(y, B, N)
end
end
y = [0.0]
N=5
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=5
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=10
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=30
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=100
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=200
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=300
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=400
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=600
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
N=1000
A = rand(N,N)
B = zeros(N,N)
@time run_julia(y,A,B,N)
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment