Thank you so much !
With your function I was able, moving the matrix on cuda, to
compute pairwise distances 210/220 times faster with respect
to scipy.spatial.distance.pdist !
I checked this on matrices from 1000 to 10000 samples of flattened
cifar10 images.
There is a small error with respect to the scipy function, I estimated
the relative error (max / median, taken over each pairwise distance)
to be around 1e-5, which is fine for me.