[CUDA] Cholesky via cuSOLVER

[sashko-zakharchuk]

Proposed changes

First op from the CUDA linalg gap discussed in #1392 (and #1026); inverse would follow.

• Cholesky::eval_gpu in the CUDA backend, backed by cuSOLVER: cusolverDnXpotrf per
  matrix, switching to cusolverDnSpotrfBatched for batches of matrices up to n = 256
  (crossover measured on sm_120). Handles are cached per device the same way as the cuBLAS
  and cuDNN ones (cusolver_utils.{h,cpp}).
• A small kernel zeroes the untouched triangle after potrf, matching the CPU op's output
  exactly.
• info is allocated but never read back: reading it costs a sync, and the CPU op also
  ignores a positive info, so neither path reports a non positive definite input.
• linalg::cholesky now accepts a GPU stream when the CUDA backend is available. Metal
  still raises at graph construction with the same message as before.
• Wheel packaging: nvidia-cusolver added to install_requires, the auditwheel excludes,
  and the MLX_LOAD_CUDA_LIBS_FROM_PYTHON rpaths (the cu12 cusolver wheel resolves its
  cusparse/nvJitLink deps through its own rpath, so no further pins are needed).
• Windows: the CI toolkit install gains the cusolver subpackages, and the delay-load helper
  learns to resolve cusolver, registering the cusparse/nvjitlink wheel dirs alongside it. I
  have no Windows machine, so that path is only compile tested.
• On CUDA builds the tests add GPU checks against the CPU result on positive definite
  inputs: a single 3x3, 16 8x8 through the batched path, two 512x512 through the loop, plus
  empty and non contiguous inputs.

float64 stays CPU-only: GPU streams reject float64 at array construction, so the GPU path
only ever sees float32. Non contiguous inputs go through the copy that already runs before
the factorization, so the kernels always get dense row major matrices.

Benchmarks

RTX 5050 (sm_120), float32, against the CPU path on the same machine (Threadripper PRO
5975WX):

           shape    cpu ms    gpu ms  speedup
           16x16     0.195     0.044     4.5x
           64x64     0.035     0.080     0.4x
         128x128     0.190     0.132     1.4x
         256x256     0.436     0.218     2.0x
         512x512     1.360     0.384     3.5x
       1024x1024     6.849     0.820     8.4x
       2048x2048    13.494     1.919     7.0x
       4096x4096    44.960     4.645     9.7x
        64x16x16     0.240     0.062     3.9x
        64x64x64     0.521     0.125     4.2x
       256x32x32     0.685     0.085     8.0x
      16x256x256     4.608     0.369    12.5x

A single 64x64 is the one shape measured where the CPU is still faster. The same sweep on an
RTX PRO 6000 (GB202) lands within noise of these numbers, and the batched/loop threshold held
on both cards.

Beyond the updated unit tests, a 60-case differential run against the CPU implementation
(sizes 1 to 257, three batch shapes, both triangles, non contiguous input, empty, non
positive definite) matches everywhere at float32 tolerances.

Two behavior notes from stress testing:

• For positive semi definite input (an exact zero eigenvalue) the undefined region differs:
  LAPACK leaves finite garbage past the rank boundary, cuSOLVER usually writes NaN from that
  row on, and whether it does varies by version. The valid leading block agrees to about
  1e-5. Worth knowing because test_cholesky's matrix is singular (sqrtA there has rank
  2): on it cuSOLVER 12.6 writes NaN into the upper factor while 12.9 and 13.0 do not, so
  the new GPU checks use positive definite inputs instead.
• Two python threads running cholesky on separate mx.new_stream streams intermittently
  poison stream capture (cudaStreamEndCapture ... previous error during capture, roughly
  half of runs). Serializing our captures behind a mutex does not change the rate, and the
  same two-thread pattern with matmul does not fail at all, so I do not think it is the
  cholesky call itself. It does not happen single threaded, with threads sharing a stream,
  or with MLX_USE_CUDA_GRAPHS=0. I can open a separate issue with the repro.

Checklist

• I have read the CONTRIBUTING document
• I have run pre-commit run --all-files to format my code / installed pre-commit prior to committing changes
• I have added tests that prove my fix is effective or that my feature works
• I have updated the necessary documentation (if needed)

[sashko-zakharchuk]

Force-pushed: rebased on main, and fixed the cuda-12.6 failure.

The GPU assertions I had added to test_cholesky ran against that test's existing matrix,
sqrtA.T @ sqrtA / 81, which is singular (sqrtA is rank 2), so the factor past the rank
boundary is undefined. cuSOLVER 12.6 writes NaN there where 12.9 and 13.0 return finite
values, which is why only that one job went red and only on the upper=True assertion.

The existing assertions are back to untouched upstream code on the CPU stream, and the GPU
checks now use positive definite inputs sized to hit both cuSOLVER paths.

[HaoXuAI]

Nice work — I'd independently implemented the same op before this landed (my PR was closed as a duplicate, correctly). Two things from my testing that might be useful, and one thing yours does better than mine did.

The n <= 256 batched threshold may be hardware-dependent. On an L40S (sm_89, datacenter) potrfBatched stays substantially ahead of the loop well past n=256:

shape	batched	serialised loop
256 × 32²	0.089 ms	3.856 ms
64 × 128²	0.223 ms	4.919 ms
16 × 512²	0.721 ms	5.304 ms

So 16 x 512² would be ~7x slower routed through the loop on this GPU. Your crossover was measured on sm_120, which has a very different SM count and bandwidth — might be worth either raising the threshold or gating it on device properties.

The fill mode is asymmetric in cuSOLVER. upper=False (the MLX default) maps to CUBLAS_FILL_MODE_UPPER, which is consistently slower than FILL_MODE_LOWER at every size I measured — 512²: 2.04x vs 1.42x, 2048²: 1.96x vs 1.09x, 4096²: 1.97x vs 0.90x, all relative to torch.linalg.cholesky on the same GPU. Computing the fast mode and transposing is mathematically equivalent (cholesky(A, upper=True).T == cholesky(A, upper=False), verified numerically) but only pays off at large n: +13% at 4096², -8% at 512². Probably not worth doing, but worth knowing the default path is the slower one.

I also have a PyTorch comparison benchmark if that's useful — single matrices land at ~2x torch.linalg.cholesky, batched at 1.0-1.7x. Happy to open that separately or hand it over.

One note for my own benefit: launching the pointer-fill kernel inside the capture context so stream order handles the ordering is neater than what I did (allocating the pointer array as an mlx array just to get a graph dependency edge). Stealing that.