support numpy/2 and cupy/14 linalg.solve shape requirements

[sbailey]

Numpy/2.x and cupy/14.x changed the requirements on the allowable input shapes to linalg.solve(M, y), leading to a redrock.zscan.solve_matrices crash when using cupy/14.x . We had not caught this with previous numpy/2.x testing because the CPU path calls solve_matrices differently than the GPU path and thus didn't hit the API change. Unit tests fail when running on GPUs at NERSC with cupy/14.x (DESI "test-main" environment); this PR fixes that by reshaping to have an extra dimension=1 and then re-reshaping.

Example

srun -N 1 -n 64 -c 2 --gpu-bind=map_gpu:3,2,1,0 --cpu-bind=cores \
    rrdesi_mpi --gpu --max-gpuprocs 4 \
    -i /global/cfs/cdirs/desi/users/sjbailey/spectro/redux/test-26.2/tiles/cumulative/42508/20220612/coadd-0-42508-thru20220612.fits \
    -o $SCRATCH/desi/temp/redrock-0-42508-thru20220612.fits

calls solve_matrices(M,y) with M.shape=(1446, 10, 10) and y.shape=(1446, 10) which is then passed to cp.linalg.solve(M,y). Previously this was allowable, but now needs y.shape=(1446,10,1). I made the reshaping change internal to solve_matrices so that its API remains unchanged.

This update is required for Matterhorn. I have verified that with the current desi main environment (numpy/1.22.4, cupy/13.1.0) this PR produces bitwise identical output as Redrock main. With the new desi test-main environment (numpy/2.3.5, cupy/14.0.1), Redrock main crashes, while this PR produces output that differs by some rounding but no substantial differences compared out our current environment.

Background:

• https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve
• Make cupy.linalg.solve compatible with numpy v2 cupy/cupy#8629

[coveralls]

coverage: 42.195% (-0.04%) from 42.233%
when pulling c89e582 on cupy_solve
into 3aa17a7 on main.

[Copilot]

Pull request overview

Fixes a crash in redrock.zscan.solve_matrices when using numpy 2.x / cupy 14.x, which tightened linalg.solve RHS shape requirements for batched solves (GPU path).

Changes:

• Normalize solve_algorithm once (upper()), simplifying downstream comparisons.
• Reshape y to (..., m, 1) for PCA solves before calling np/cp.linalg.solve, then reshape back to preserve the existing solve_matrices API.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[sbailey]

@craigwarner-ufastro reviewed and commented by Slack:

Looks like it is correct to me.  Another valid fix would be
return cp.linalg.solve(M, y[:, :, cp.newaxis]).squeeze(-1) or to generalize
return cp.linalg.solve(M, y[...,cp.newaxis]).squeeze(-1)
but I don't think it matters which of these we choose - doesn't seem like it changes the timing at all so may as well leave the version you already have for the PR.  But I did test that it matches cp.linalg.solve(M, y) in an example and all are == (on cupy 13.1).

.squeeze is a good trick to know, but would be dangerous in the case of processing a file with 1 target — I think that would get squeezed out as well. So I'm going to take that Slack comment as PR review and approval and merge this.