Ridge with cholesky solver returns NaN

Describe the bug

When using Ridge with a big number of features and the Cholesky solver, the model gets full of NaNs.
Playing with the number of features we can see that the imprecision grows until it becomes NaN.

The bug is architecture-dependent.
The provided snippet works perfectly in an Intel(R) Xeon(R) CPU E3-1585L v5 @ 3.00GHz.
It fails when running in an Intel(R) Xeon(R) Gold 6140 CPU @ 2.30GHz.

A workaround was found by blindly exporting OPENBLAS_CORETYPE="Broadwell".

Steps/Code to Reproduce

import numpy as np
from sklearn.linear_model import Ridge

# we need something big to reproduce:
n_inputs = 1000
n_targets = 1
n_samples = 20000

rng = np.random.default_rng(seed=42)

b = rng.random((n_targets,))
w = rng.random((n_targets, n_inputs)) / n_inputs

X = rng.random((n_samples, n_inputs))
y = rng.random((n_samples, n_targets))

for i in range(n_samples):
    y[i] = (w @ X[i]) + b

# "lsqr" does not have the same problem as "cholesky"
model = Ridge(alpha=0.0, tol=1e-6, solver="cholesky")
model.fit(X, y)

np.testing.assert_almost_equal(model.coef_, w)
np.testing.assert_almost_equal(model.intercept_, b)

error = np.abs(model.predict(X) - y).sum()
np.testing.assert_almost_equal(error, 0.0, decimal=4)

assert np.isfinite(model.coef_).all()
assert np.isfinite(model.intercept_).all()

print("done!")

Expected Results

done!

Actual Results

Numpy assertion errors as model.coef_ is full of NaN.

Versions

The software is the same on both machines:

System:
    python: 3.8.9 (default, Apr  2 2021, 11:20:07)  [GCC 10.3.0]
executable: ***/.venv/bin/python
   machine: Linux-3.10.0-1062.9.1.el7.x86_64-x86_64-with-glibc2.2.5

Python dependencies:
          pip: 20.2.3
   setuptools: 54.2.0.post0
      sklearn: 0.24.1
        numpy: 1.20.2
        scipy: 1.6.1
       Cython: 0.29.22
       pandas: 1.2.3
   matplotlib: 3.4.1
       joblib: 1.0.1
threadpoolctl: 2.1.0

Built with OpenMP: True

Start by running the provided Ridge reproduction on the affected and working CPU configurations, then inspect the Ridge cholesky solver path and the OpenBLAS behavior. Done means the reproduction passes with finite coefficients and intercepts without exporting OPENBLAS_CORETYPE.