4 ms·

I think Jake was pretty straight forward in saying in the post that he isn't a fortran expert and was looking for someone to offer up a better version if he was

• StefanKarpinski 13y ago
  The text is definitely clear about it if you read everything, but the chart could easily be misread as saying that Numba and Cython are 2x faster than hand-written Fortran. That's my main objection. I also think that writing a simple C version of this code would serve as a much better baseline than f2py-generated Fortran code. That was definitely an inappropriate labeling of the Python benchmarks on the Julia home page. We fixed it as soon as it was pointed out and I'm grateful to whoever did point that out. I'm just doing the same thing here. Dahua's explanation is very cool and he's really created an incredibly high-quality, comprehensive package for distance calculations. It's great work.
• mitmatt 13y ago
  That computation is probably something you were already familiar with in disguise: it's just the law of cosines [1][2] along with the fact that a matrix of inner products can be computed with a matrix-matrix multiplication [3]. The first claim is just that for vectors x and y in an inner-product space ||x-y||^2 = <x-y,x-y> = ||x||^2 + ||y||^2 - 2<x,y>, and the second claim is just that if you collect a bunch of vectors into the columns of X and another bunch into the columns of Y, then the (i,j) element of X'Y is <x_i,y_j>, i.e. the inner product of the i'th vector in X with the j'th vector in Y. (Bonus fact: you can use that relationship to translate the statement "this set of pairwise distances can be embedded in a Euclidean space" to an equivalent condition that some simple matrix is negative semidefinite on a subspace orthogonal to the all-ones vector. Euclidean metric embedding is very useful!) [1] http://en.wikipedia.org/wiki/Law_of_cosines#Vector_formulation http://en.wikipedia.org/wiki/Law_of_cosines#Vector_formulati... [2] http://en.wikipedia.org/wiki/Polarization_identity http://en.wikipedia.org/wiki/Polarization_identity [3] http://en.wikipedia.org/wiki/Gramian_matrix http://en.wikipedia.org/wiki/Gramian_matrix
• gus_massa 13y ago
  At least, he is traveling through the Fortran arrays in the wrong order. In Fortran arrays are stored in memory in the opposite order of C, Fortran order: A(1,1); A(2,1); A(1,2); A(2,2); C order: A[1,1]; A[1,2]; A[2,1]; A[2,2]; ( See for example: http://www.xlsoft.com/jp/products/intel/cvf/docs/vf-html_e/pg/pguaracc.htm http://www.xlsoft.com/jp/products/intel/cvf/docs/vf-html_e/p... ) Compiling in gfortran with -O3 sometimes solves that automatically for you, but you must try to use the correct order. And with less optimization the difference in the order of the index can produce very big difference in time dew to the cache problems. Also, the C/phyton program uses a few trick like += and tmp variables. I don't know how they interact with the optimizations, -O3 does very strange things and perhaps they generate the same code.