Inconsistent results with same random seed

### Describe the bug

I'm not sure this is a bug and I couldn't find anything in the issues archive, but I'm seeing an inconsistency between consecutive runs of kmeans.fit, even when setting the same random seed via np.random.seed or random_state. I pinned it down to multithreading: if I force self._n_threads to 1 inside _kmeans.py/fit the inconsistency goes away. See test code below. Is this expected behavior, and is there a way to eliminate multithreading during fitting?

```python
from sklearn.cluster import KMeans
import numpy as np

np.random.seed(0)
matrix = np.random.randn(1000, 100)

for i in range(100):
    k1 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k1.fit(matrix)
    k2 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k2.fit(matrix)
    print('min diff: %s, max diff: %s' % (str(np.min(k1.cluster_centers_ - k2.cluster_centers_)), str(np.max(k1.cluster_centers_ - k2.cluster_centers_))))
```

### Steps/Code to Reproduce

```python
from sklearn.cluster import KMeans
import numpy as np

np.random.seed(0)
matrix = np.random.randn(1000, 100)

for i in range(100):
    k1 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k1.fit(matrix)
    k2 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k2.fit(matrix)
    print('min diff: %s, max diff: %s' % (str(np.min(k1.cluster_centers_ - k2.cluster_centers_)), str(np.max(k1.cluster_centers_ - k2.cluster_centers_))))
```

### Expected Results

min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
...

### Actual Results

min diff: -1.1102230246251565e-16, max diff: 2.220446049250313e-16
min diff: -2.220446049250313e-16, max diff: 1.1102230246251565e-16
min diff: -1.1102230246251565e-16, max diff: 1.1102230246251565e-16
min diff: 0.0, max diff: 0.0
min diff: -1.1102230246251565e-16, max diff: 1.1102230246251565e-16
min diff: -1.1102230246251565e-16, max diff: 2.220446049250313e-16
,,,

### Versions

```shell
System:
    python: 3.8.10 (tags/v3.8.10:3d8993a, May  3 2021, 11:48:03) [MSC v.1928 64 bit (AMD64)]
executable: C:\Python38\python.exe
   machine: Windows-10-10.0.19045-SP0
Python dependencies:
      sklearn: 1.3.1
          pip: 21.1.1
   setuptools: 56.0.0
        numpy: 1.24.2
        scipy: 1.8.0
       Cython: 0.29.28
       pandas: 1.3.4
   matplotlib: 3.4.3
       joblib: 1.3.2
threadpoolctl: 3.1.0
Built with OpenMP: True
threadpoolctl info:
       user_api: openmp
   internal_api: openmp
         prefix: vcomp
       filepath: C:\Python38\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None
    num_threads: 16
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Python38\Lib\site-packages\numpy\.libs\libopenblas64__v0.3.21-gcc_10_3_0.dll
        version: 0.3.21
threading_layer: pthreads
   architecture: SkylakeX
    num_threads: 16
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Python38\Lib\site-packages\scipy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
        version: 0.3.17
threading_layer: pthreads
   architecture: Prescott
    num_threads: 16
```

COMMENT:
Are you considering to add this behaviour to the documentation? 

The current documentation seems to indicate that the random_state is the only parameter required to set to get deterministic results. In my case, this caused another custom algorithm that uses the cluster centers to give indeterministic result. Using threadpoolctl  works very well to mitigate this issue, however it took considerable time to pin the problem to this issue which would be much easier if there would be a disclaimer of some sort in the documentation.

Labels: Bug