Truly parallel execution of pairwise_kernels and pairwise_distances

### Describe the workflow you want to enable

Both `pairwise_kernels` and `pairwise_distances` functions call `_parallel_pairwise` function, which is (contrary to its name) not parallel as it enforces the threading backend. Therefore, these functions are terribly slow, especially for computationally expensive user-defined metrics. I understand that the reasons for the threading backend are possibly large memory demands and data communication overhead but I suggest a different approach. Also, the documentation for these functions talks about parallel execution and processes which is currently simply not true.

### Describe your proposed solution

The memory and data communication issues can be reduced by a smarter distribution of the input data to individual processes. Right now, only `Y` is sliced in the `_parallel_pairwise` function which is suboptimal for parallel processing. Both `X` and `Y` should be sliced to lower the demands for multiprocessing. For example for 100x100 `X` and `Y` distributed to 100 processes, we have to copy 100+1 inputs to every process when slicing only `Y` while only 10+10 when slicing both `X` and `Y`. As a result, multiprocessing can be allowed. Also, joblib does automatic memmapping in some cases.

Alternatively, at least the documentation for  `pairwise_kernels` and `pairwise_distances` should be corrected.

### Describe alternatives you've considered, if relevant

_No response_

### Additional context

_No response_

COMMENT:
@ogrisel 
I totally get your point. I am afraid some complexity would have to be added to keep the performance of the threading approach for IO-bound tasks. Therefore, I will just create a PR to improve the documentation. I have added the final version of my function (with improved slicing and byte ordering) below just in case somebody is interested. As mentioned above, I use threading by default for predefined kernels and metrics and multiprocessing for user-defined functions. BTW `SharedMemory` from `multiprocessing` might also help with the overhead for multiprocessing.

Results using the function below and NOT using free threading:
<p float="left">
  <img src="https://github.com/user-attachments/assets/3ab025a4-5fee-4779-ac4c-f8f40c35b308" width="400" />
  <img src="https://github.com/user-attachments/assets/68e8d028-2288-46ed-8211-9c9873dea3f9" width="400" /> 
</p>

<p float="left">
  <img src="https://github.com/user-attachments/assets/2d942f9c-8d1b-432a-926b-f8ea12bb4930" width="400" />
  <img src="https://github.com/user-attachments/assets/4f09f531-a06d-4066-a506-c4c3174df304" width="400" /> 
</p>

```
def _parallel_pairwise4(X, Y, func, n_jobs, **kwds):
    """Break the pairwise matrix in n_jobs even slices
    and compute them in parallel."""

    if Y is None:
        Y = X

    if effective_n_jobs(n_jobs) == 1:
        return func(X, Y, **kwds)
    
    ret = np.empty((X.shape[0], Y.shape[0]), dtype=X.dtype)
    eff_jobs = effective_n_jobs(n_jobs)
    a = math.sqrt(eff_jobs/(len(X)*len(Y)))
    n_jobs1 = int(a*min(len(X), len(Y)))
    if n_jobs1 == 0:
        n_jobs1 = 1
    else:
        for i in range(n_jobs1, 0, -1):
            if eff_jobs%i==0:
                n_jobs1 = i
                break
    n_jobs2 = int(eff_jobs//n_jobs1)
    if len(Y) < len(X):
        n_jobs1, n_jobs2 = n_jobs2, n_jobs1
    slices1 = list(gen_even_slices(len(X), n_jobs1))
    slices2 = list(gen_even_slices(len(Y), n_jobs2))
    out = Parallel(n_jobs=n_jobs)(
        delayed(func)(X[s1], Y[s2], **kwds)
        for s1 in slices1 for s2 in slices2
    )
    i = 0
    for s1 in slices1:
        for s2 in slices2:
            ret[s1, s2] = out[i]
            i += 1

    return ret
```

Labels: New Feature