[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.