<jupyter_start><jupyter_text>A Speed Comparison Of C, Julia, Python, Numba, and Cython on LU Factorization [<NAME>](https://www.ibm.com/developerworks/community/blogs/jfp?lang=en)The code in this notebook has been used to generate the results discussed in [A Speed Comparison Of C, Julia, Python, Numba, and Cython on LU Factorization](https://www.ibm.com/developerworks/community/blogs/jfp/entry/a_comparison_of_c_julia_python_numba_cython_scipy_and_blas_on_lu_factorization)This study was inspired by <NAME>'s [Speed comparision Numba vs C vs pure Python at the example of the LU factorization](http://expdesign.iwr.uni-heidelberg.de/people/swalter/blog/python_numba/index.html)Timingsheavily depend on the machine you use. We start with an import of useful packages.<jupyter_code>import numpy as np
import pandas as pd
import time
from numpy.testing import *
from scipy.linalg import *

import gc
from matplotlib import pyplot as plt
%matplotlib inline

from numba import njit, void, double<jupyter_output><empty_output><jupyter_text>Pure Python<jupyter_code>def det_by_lu(y, x):
    y[0] = 1.

    N = x.shape[0]
    for k in range(N):
        y[0] *= x[k,k]
        for i in range(k+1, N):
            x[i,k] /= x[k,k]
            for j in range(k+1, N):
                x[i,j] -= x[i,k] * x[k,j]
def run_python(A,B,y,N):

    # check that result is correct
    np.copyto(B,A)
    det_by_lu(y, B)

    L = np.tril(B, -1) + np.eye(N)
    U = np.triu(B)
    assert_almost_equal( L.dot(U),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + (100000 // (N*N))
    
    for l in range(loops):
        np.copyto(B,A)
        det_by_lu(y, B)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Numpy<jupyter_code>def numpy_det_by_lu(y, x):
    y[0] = 1.

    N = x.shape[0]
    with np.errstate(invalid='ignore'):
        for k in range(N):
            y[0] *= x[k,k]
            xk = x[k]
            for i in range(k+1, N):
                xi = x[i]
                xi[k] /= xk[k]
                xi[k+1:] -= xi[k] * xk[k+1:]
def run_numpy(A,B,y,N):

    # check that result is correct
    np.copyto(B,A)
    numpy_det_by_lu(y, B)

    L = np.tril(B, -1) + np.eye(N)
    U = np.triu(B)
    assert_almost_equal( L.dot(U),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + (100000 // (N*N))
    
    for l in range(loops):
        np.copyto(B,A)
        numpy_det_by_lu(y, B)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Numba<jupyter_code>def numba_det_by_lu(y, x):
    y[0] = 1.

    N = x.shape[0]
    for k in range(N):
        y[0] *= x[k,k]
        for i in range(k+1, N):
            x[i,k] /= x[k,k]
            for j in range(k+1, N):
                x[i,j] -= x[i,k] * x[k,j]
            

fastdet_by_lu = njit(void(double[:], double[:,:]))(numba_det_by_lu)
def run_numba(A,B,y,N):

    # check that result is correct
    np.copyto(B,A)
    fastdet_by_lu(y, B)

    L = np.tril(B, -1) + np.eye(N)
    U = np.triu(B)
    assert_almost_equal( L.dot(U),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + min(1000000 // (N*N), 20000)
    
    for l in range(loops):
        np.copyto(B,A)
        fastdet_by_lu(y, B)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Cython<jupyter_code>%load_ext cython
%%cython --verbose --compile-args=-march=native --compile-args=-ffast-math

import cython

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef cython_det_by_lu(double[:] y, double[:,:] x):
    y[0] = 1.

    cdef int N = x.shape[0]
    cdef int i,j,k
    
    for k in range(N):
        y[0] *= x[k,k]
        for i in range(k+1, N):
            x[i,k] /= x[k,k]
            for j in range(k+1, N):
                x[i,j] -= x[i,k] * x[k,j]
def run_cython(A,B,y,N):

    # check that result is correct
    np.copyto(B,A)
    cython_det_by_lu(y, B)

    L = np.tril(B, -1) + np.eye(N)
    U = np.triu(B)
    assert_almost_equal( L.dot(U),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + min(1000000 // (N*N), 20000)
    for l in range(loops):
        np.copyto(B,A)
        cython_det_by_lu(y, B)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>C<jupyter_code>!gcc -shared -O3 -march=native -ffast-math lu.c -o lu.so
from cffi import FFI
ffi = FFI()
ffi.cdef('void det_by_lu(double *y, double *B, int N);')
C = ffi.dlopen("lu.so")
c_det_by_lu = C.det_by_lu
def run_c(A,B,y,N):
    # run c code
    #B = numpy.zeros((N,N), order='F')
    #B[:,:] = A
    np.copyto(B,A)
    c_det_by_lu(ffi.cast("double *", y.ctypes.data),
                ffi.cast("double *", B.ctypes.data),
                ffi.cast("int", N))
    
    # check that result is correct
    L = np.tril(B, -1) + np.eye(N)
    U = np.triu(B)
    assert_almost_equal( L.dot(U),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + min(1000000 // (N*N), 20000)
    
    for l in range(loops):
        np.copyto(B,A)
        c_det_by_lu(ffi.cast("double *", y.ctypes.data),
                    ffi.cast("double *", B.ctypes.data),
                    ffi.cast("int", N))
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Scipy<jupyter_code>def run_scipy(A,B,y,N):

    # check that result is correct
    np.copyto(B,A)
    (P,L,U) = lu(B,overwrite_a=True)

    assert_almost_equal( P.dot(L.dot(U)),  A)
    
    gc.disable()
    st = time.time()
    
    loops = 1 + min(1000000 // (N*N), 20000)

    for l in range(loops):
        np.copyto(B,A)
        lu(B,overwrite_a=True)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Lapack<jupyter_code>def run_lapack(A,B,y,N):

    # check that result is correct
    
    gc.disable()
    st = time.time()
    
    loops = 1 + min(1000000 // (N*N), 20000)

    for l in range(loops):
        np.copyto(B,A)
        lu_factor(B,overwrite_a=True)
        
    et = time.time()
    gc.enable()
    
    return  (et - st)/loops<jupyter_output><empty_output><jupyter_text>Timings<jupyter_code>def timings(n=7, \
            series=['pure python', 'c', 'numba', 'numpy', \
                    'cython', 'scipy', 'lapack', 'julia']):
    
    Ns = np.array([5,10,30,100, 200, 300, 400, 600, 1000, 2000, 4000, 8000])
    Fs = [run_python, run_c, run_numba, run_numpy, \
          run_cython, run_scipy, run_lapack]
    times = pd.DataFrame(np.zeros((n, len(Fs)+1)), index = Ns[:n], columns = series)
    for i,N in enumerate(Ns[:n]):
        print ('N =', N, end=" ")
        A = np.random.random((N,N))
        B = np.empty(A.shape)
        y = np.zeros(1)
        for j,label in enumerate(series[:-1]):
            if label != '': 
                print(j, end=" ")
                times.loc[N,label] = Fs[j](A,B,y,N)
        print('')
    return times<jupyter_output><empty_output><jupyter_text>Running the code<jupyter_code>times = timings()<jupyter_output>N = 5 0 1 2 3 4 5 6 
N = 10 0 1 2 3 4 5 6 
N = 30 0 1 2 3 4 5 6 
N = 100 0 1 2 3 4 5 6 
N = 200 0 1 2 3 4 5 6 
N = 300 0 1 2 3 4 5 6 
N = 400 0 1 2 3 4 5 6<jupyter_text>Checking results<jupyter_code>times<jupyter_output><empty_output><jupyter_text>Adding Julia timings.Here is the output of the julia code<jupyter_code>!julia lu.jl
julia_times = [1.75764045e-7, 6.4414374e-7, 1.270633365e-5, 0.00042416499000000004, 0.003296228324, 0.010970920214999999, 0.02566115192, 0.08423355495599999, 0.3971115121]<jupyter_output><empty_output><jupyter_text>We add these into our timings data frame<jupyter_code>M = min(len(julia_times),times.shape[0])
times['julia'] = julia_times[:M]<jupyter_output><empty_output><jupyter_text>Checking results again<jupyter_code>times<jupyter_output><empty_output><jupyter_text>Display<jupyter_code>def plot_times(times, 
               cols = [],\
               name="runtimes.png"):
    plt.figure(figsize=(10,6))
    if cols == []:
        cols = times.columns
    for i,label in enumerate(cols):
        if label != '':
            plt.loglog(times.index, times[label], label=label)
    plt.xlabel("N (matrix size)")
    plt.ylabel("runtime [sec]")
    plt.grid()
    plt.legend(loc=2)
    plt.savefig(name)
    plt.show()
plot_times(times, cols=['pure python', 'numpy', 'scipy', 'lapack', 'c', 'cython',  'numba', 'julia'], name='runtimes_1')
plot_times(times, cols=['pure python', 'c', 'cython'], name='runtimes_2')
plot_times(times, cols=['pure python', 'c', 'numpy'], name='runtimes_3')
plot_times(times, cols=['pure python', 'c', 'scipy', 'lapack'], name='runtimes_4')
plot_times(times, cols=['pure python', 'c', 'julia'], name='runtimes_5')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15871054700662698, 'fraction_numerical': 0.03931259126137257, 'mean_word_length': 2.603399433427762, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 62, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '346572', 'license': 'Apache-2.0', 'n_tokens_mistral': 3995, 'n_tokens_neox': 3680, 'n_words': 753, 'path': 'benchmarks/lu_decomposition/LU_decomposition.ipynb', 'repo_name': 'ashwinvis/sandbox', 'size': '8929'}