Hi Phillip, Others may respond with more specific answers, but have you had the chance to read through the Julia performance tips in the Julia manual?
http://julia.readthedocs.org/en/latest/manual/performance-tips/ Cheers, Kevin On Monday, August 25, 2014, Phillip Berndt <phillip.ber...@googlemail.com> wrote: > Hi julia-users, > > I've recently stumbled over Julia and wanted to give it a try. > > To assess it's speed, I've implemented another micro-benchmark, namely a > version of Matlab's magic() function that generates magic squares. Since I > have no experience writing optimal Julia code, I started off with literal > translations of two different implementations - Matlab's and the one from > magic_square.py from PyPy, which is an optimized version for NumPy. I then > timed the calculation of all magic squares from N=3 to N=1000. The table > from Julia's homepage suggests that in most cases, it is significantly > faster than Python and Matlab. In my case, it's significantly slower, which > is somehow disappointing ;) My question now is: > > Can the implementation be optimized to outperform the other two? > > *The times:* > > Julia, Matlab version: elapsed time: 18.495374216 seconds (13404087428 > bytes allocated, 12.54% gc time) > Julia, Python version: elapsed time: 8.107275449 seconds (13532473792 > bytes allocated, 26.99% gc time) > Matlab: Elapsed time is 4.994960 seconds. > Python: 1 loops, best of 3: 2.09 s per loop > > My test machine is a 4 Core i7-4600 Notebook with 2.1 GHz and 8 GiB RAM, > running a current Linux Mint and Julia 0.3 stable. To be fair, Python does > not seem to gc during this loop (disabling gc doesn't alter the time here), > so one should compare with 8.1 s * (1.-.2699) = 5.91 s for Julia. That's > still much slower than Python. (By the way, even Octave only needs 4.46 > seconds.) If I translate the matrices in magic_python to account for > column-major storage, the execution time does not significantly improve. > > *The code:* > > Matlab: tic; arrayfun(@magic, 3:1000, 'UniformOutput', false); toc > IPython: import magic_square; %timeit [ magic_square.magic(x) for x in > range(3, 1001) ]; > Julia: I've uploaded the code to a Gist at > https://gist.github.com/phillipberndt/2db94bf5e0c16161dedc and will paste > a copy below this post. > > > Cheers, > Phillip > > > function magic_matlab(n::Int64) > # Works exactly as Matlab's magic.m > > if n % 2 == 1 > p = (1:n) > M = n * mod(broadcast(+, p', p - div(n+3, 2)), n) + > mod(broadcast(+, p', 2p - 2), n) + 1 > return M > elseif n % 4 == 0 > J = div([1:n] % 4, 2) > K = J' .== J > M = broadcast(+, [1:n:(n*n)]', [0:n-1]) > M[K] = n^2 + 1 - M[K] > return M > else > p = div(n, 2) > M = magic_matlab(p) > M = [M M+2p^2; M+3p^2 M+p^2] > if n == 2 > return M > end > i = (1:p) > k = (n-2)/4 > j = convert(Array{Int}, [(1:k); ((n-k+2):n)]) > M[[i; i+p],j] = M[[i+p; i],j] > i = k+1 > j = [1; i] > M[[i; i+p],j] = M[[i+p; i],j] > return M > end > end > @vectorize_1arg Int magic_matlab > > function magic_python(n::Int64) > # Works exactly as magic_square.py (from pypy) > > if n % 2 == 1 > m = (n >> 1) + 1 > b = n^2 + 1 > > M = reshape(repmat(1:n:b-n, 1, n+2)[m:end-m], n+1, n)[2:end, :] + > reshape(repmat(0:(n-1), 1, n+2), n+2, n)[2:end-1, :]' > return M > elseif n % 4 == 0 > b = n^2 + 1 > d = reshape(1:b-1, n, n) > > d[1:4:n, 1:4:n] = b - d[1:4:n, 1:4:n] > d[1:4:n, 4:4:n] = b - d[1:4:n, 4:4:n] > d[4:4:n, 1:4:n] = b - d[4:4:n, 1:4:n] > d[4:4:n, 4:4:n] = b - d[4:4:n, 4:4:n] > d[2:4:n, 2:4:n] = b - d[2:4:n, 2:4:n] > d[2:4:n, 3:4:n] = b - d[2:4:n, 3:4:n] > d[3:4:n, 2:4:n] = b - d[3:4:n, 2:4:n] > d[3:4:n, 3:4:n] = b - d[3:4:n, 3:4:n] > > return d > else > m = n >> 1 > k = m >> 1 > b = m^2 > > d = repmat(magic_python(m), 2, 2) > > d[1:m, 1:k] += 3*b > d[1+m:end, 1+k:m] += 3*b > d[1+k, 1+k] += 3*b > d[1+k, 1] -= 3*b > d[1+m+k, 1] += 3*b > d[1+m+k, 1+k] -= 3*b > d[1:m,1+m:n-k+1] += b+b > d[1+m:end, 1+m:n-k+1] += b > d[1:m, 1+n-k+1:end] += b > d[1+m:end, 1+n-k+1:end] += b+b > > return d > end > end > @vectorize_1arg Int magic_python > > print("Matlab version: ") > @time magic_matlab(3:1000) > > print("Python version: ") > @time magic_python(3:1000) > > >