Instantly share code, notes, and snippets. jfpuget/LU_decomposition.ipynb Created Jan 27, 2016

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. for k = 1:N y *= 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) Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

dilawar commented Dec 31, 2016

 I would be cool to see how PyPy performs.

unyty commented Mar 20, 2018

 A perfect exercise for beginners

t184256 commented May 6, 2018 • edited

 But why do you compile Cython code without optimizations?
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.