MAINT: Implement lstsq as a gufunc

[eric-wieser]

A summary of what this does

• Includes the missing s and c varaiants of gelsd (not used by lstsq, but there if people need them) (the very large diffs are simply what happens when regenerating lapack_lite)
• Writes the C wrappers that call the lapack functions and allocate/deallocate their workspaces
• Changes the implementation of lstsq to fall back on one of two gufuncs (that share a loop)

What this does not try to do yet, since there are decisions to be made for these:

• Enable broadcasting stsacks of matrices (ENH: broadcast lstsq #8720)
• Enable arrays with a 0 core dimension

Like #9976 (merged), this needs lapack 3.2.2.

This was tested for exact equality by placing asserts within lstsq comparing old and new values, and rerunning the test suite.

[mhvk]

I very much like this! And it all looks great. Ideally, we allow lstsq now to take single precision inputs, but that is probably better done in a follow-up PR as well (not sure how best to handle it; new argument to tell whether to allow single precision?)

Relatedly, since we expose the single-precision routines, we would ideally also test those, but I think that is also fine to leave to a next PR.

[mhvk]

Is this necessary given that the inputs are coerced already?

[eric-wieser]

The inputs don't look coerced already to me

[mhvk]

Ah, yes, above one just gets the types for later conversion. Should have noticed that.

[eric-wieser]

Ideally, we allow lstsq now to take single precision inputs

I think we're stuck with not doing this, in the same way that svd is

[eric-wieser]

An earlier version exposed twice as many ufuncs, lstsq_m, and lstsq1_m where the latter was a specialization for just vectors. There's precedent for this in solve, but it seemed clearer just to insert extra dimensions within the python code.

[eric-wieser]

Once this goes in, I have another commit that promotes resids to a proper ufunc return, which saves some memory, but at the expense of clarity.

[mhvk]

I think in principle we could add a dtype argument where the default is float or complex but can be explicitly set to whatever one wants. Or simply a allow_single_precision or something like it. Possibly, a linalg wide context manager could do the trick too. Anyway, for later...

[eric-wieser]

Possibly, a linalg wide context manager could do the trick too

Of course, then we run into #9444.

I'd be in favor of allowing single-precision via dtype kwarg, with a default of dtype='promote' - passing dtype = None would override this to do the normal ufunc thing. Lets leave this till another time though.

[eric-wieser]

Rebased, and slightly squashed. Good to go, @mhvk?

[mhvk]

Yes, this looks good, modulo that there should at least be new issues reminding us to introduce a dtype option to get single-precision and tests for the single-precision routines (which can be addressed in the same PR of course).

[eric-wieser]

Would you consider #9516 to be that issue?

[mhvk]

I think a new one is better (more focused); see #10888 . Will now merge this one.