Now that we have access to BLAS functions cython bindings from scipy, we have access to level 3 functions which weren't present in vendored CBLAS. It would be interesting to check if they could be used in some places instead of a loop of lower level functions.

For example, in cd_fast I saw that:

# XtA = np.dot(X.T, R) - beta * w
for i in range(n_features):
    XtA[i] = (_dot(n_samples, &X[0, i], 1, &R[0], 1)
              - beta * w[i])

loop of dot (level 1) which could be replaced by a gemv (level 2)