I haven't done a proper benchmark yet but I saw much better runtime on my laptop. I'll do a benchmark today and post the results here.

I ran 10 times the following benchmark on a Google Cloud instance with 64 vCPUs:

from sklearn.datasets import make_blobs
from sklearn.neighbors import NearestNeighbors
from sklearn.externals.joblib import cpu_count
import time

d = make_blobs(100000, 100)[0]
nn = NearestNeighbors().fit(d)

for n_jobs in range(1, cpu_count() + 1):
    nn.n_jobs = n_jobs
    start = time.time()
    nn.radius_neighbors()
    end = time.time()
    print('{},{}'.format(n_jobs, end - start))

Although the scaling is not linear it definitely runs faster:

One thing that was not obvious to me but made a huge difference was that in _query_radius_single doing this:

-                        raise ValueError("Fatal: count out of range. "
-                                         "This should never happen.")
+                        with gil:
+                            raise ValueError("Fatal: count out of range. "
+                                             "This should never happen.")

instead of this:

-                        raise ValueError("Fatal: count out of range. "
-                                         "This should never happen.")
+                        return -1

have a completely different scaling behavior.

In both case the function is nogil and I expected that in the first case the GIL would be acquired only if the with gil: statement is reached. But as explained here any function with a with gil: statement will have to acquire the GIL on return, regardless of whether with gil: was reached or not.

@jnothman Thank you for your comments. I think this is ready for review.

Nice work ! The code seems reasonable, though I am a novice in C memory allocation.

As a general comment, your code could benefit from more comments, especially near non-standard function like np.PyArray_SimpleNewFromData, np.PyArray_UpdateFlags, or memcpy.

For instance, you could add comments with the equivalent python expressions:
# equivalent to: distances[i] = np_dist_arr[:counts[i]].copy()

@TomDLT Thank you for your comment. I added some comments as you suggested.

@TomDLT thank you for the review. I updated the release notes.

Merge when green