`Xgesvdp!` failure when only requesting singular values

[Jutho]

Describe the bug

I am trying to use the CUSOLVER Xgesvdp! method, which computes the singular value decomposition via polar decomposition and symmetric eigenvalue decomposition of the positive factor. It works when requesting both singular values and singular vectors, but not when only requesting singular values. I did not manage to decipher what is causing this bug:

To reproduce

The Minimal Working Example (MWE) for this bug:

julia> using CUDA, LinearAlgebra

julia> A = CUDA.randn(Float64, 10, 10);

julia> U, S, V, err = CUDA.CUSOLVER.Xgesvdp!('V', 0, copy(A));

julia> A ≈ U * Diagonal(S) * V' # succes!
true

julia> S, err = CUDA.CUSOLVER.Xgesvdp!('N', 1, copy(A));
ERROR: CUSOLVERError: an invalid value was used as an argument (code 3, CUSOLVER_STATUS_INVALID_VALUE)
Stacktrace:
 [1] throw_api_error(res::CUDA.CUSOLVER.cusolverStatus_t)
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/libcusolver.jl:14
 [2] check
   @ ~/.julia/packages/CUDA/oymHm/lib/cusolver/libcusolver.jl:27 [inlined]
 [3] cusolverDnXgesvdp
   @ ~/.julia/packages/GPUToolbox/cZlg7/src/ccalls.jl:33 [inlined]
 [4] (::CUDA.CUSOLVER.var"#1476#1478"{…})(buffer_gpu::CuArray{…}, buffer_cpu::Vector{…})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/dense_generic.jl:326
 [5] with_workspaces(f::CUDA.CUSOLVER.var"#1476#1478"{…}, cache_gpu::CuArray{…}, cache_cpu::Vector{…}, size_gpu::UInt64, size_cpu::UInt64)
   @ CUDA.APIUtils ~/.julia/packages/CUDA/oymHm/lib/utils/call.jl:89
 [6] Xgesvdp!(jobz::Char, econ::Int64, A::CuArray{Float64, 2, CUDA.DeviceMemory})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/dense_generic.jl:325
 [7] top-level scope
   @ REPL[17]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> S, err = CUDA.CUSOLVER.Xgesvdp!('N', 0, copy(A));
ERROR: CUSOLVERError: an invalid value was used as an argument (code 3, CUSOLVER_STATUS_INVALID_VALUE)
Stacktrace:
 [1] throw_api_error(res::CUDA.CUSOLVER.cusolverStatus_t)
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/libcusolver.jl:14
 [2] check
   @ ~/.julia/packages/CUDA/oymHm/lib/cusolver/libcusolver.jl:27 [inlined]
 [3] cusolverDnXgesvdp
   @ ~/.julia/packages/GPUToolbox/cZlg7/src/ccalls.jl:33 [inlined]
 [4] (::CUDA.CUSOLVER.var"#1476#1478"{…})(buffer_gpu::CuArray{…}, buffer_cpu::Vector{…})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/dense_generic.jl:326
 [5] with_workspaces(f::CUDA.CUSOLVER.var"#1476#1478"{…}, cache_gpu::CuArray{…}, cache_cpu::Vector{…}, size_gpu::UInt64, size_cpu::UInt64)
   @ CUDA.APIUtils ~/.julia/packages/CUDA/oymHm/lib/utils/call.jl:89
 [6] Xgesvdp!(jobz::Char, econ::Int64, A::CuArray{Float64, 2, CUDA.DeviceMemory})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/oymHm/lib/cusolver/dense_generic.jl:325
 [7] top-level scope
   @ REPL[18]:1
Some type information was truncated. Use `show(err)` to see complete types.

Expected behavior

Get a vector of singular values, and some error measure intrinsic to this method.

Version info

Details on CUDA:

julia> CUDA.versioninfo()
CUDA runtime 12.8, artifact installation
CUDA driver 12.4
NVIDIA driver 550.120.0

CUDA libraries: 
- CUBLAS: 12.8.4
- CURAND: 10.3.9
- CUFFT: 11.3.3
- CUSOLVER: 11.7.3
- CUSPARSE: 12.5.8
- CUPTI: 2025.1.1 (API 26.0.0)
- NVML: 12.0.0+550.120

Julia packages: 
- CUDA: 5.7.3
- CUDA_Driver_jll: 0.12.1+1
- CUDA_Runtime_jll: 0.16.1+0

Toolchain:
- Julia: 1.11.5
- LLVM: 16.0.6

1 device:
  0: NVIDIA TITAN V (sm_70, 11.351 GiB / 12.000 GiB available)

Details on Julia:

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 20 × Intel(R) Core(TM) i9-7900X CPU @ 3.30GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake-avx512)
Threads: 1 default, 0 interactive, 1 GC (on 20 virtual cores)

[maleadt]

CUSOLVER_STATUS_INVALID_VALUE

    Invalid parameters were passed (m,n<0 or lda<max(1,m) or ldu<max(1,m) or ldv<max(1,n) ).

So maybe ldu and ldv still have to be valid, i.e., m and n instead of 1? Pytorch also does this (for a different but similar API): https://github.com/pytorch/pytorch/blob/ad81eeb7c7c906e0cdd04a5cc8fdb9592281c317/aten/src/ATen/native/cuda/linalg/BatchLinearAlgebraLib.cpp#L442-L443

Trying that out results in illegal memory accesses. Of course, NVIDIA's example only returns everything: https://github.com/NVIDIA/CUDALibrarySamples/blob/master/cuSOLVER/Xgesvdp/cusolver_Xgesvdp_example.cu. However, looking at Pytorch again it looks like we still have to allocate U and V 🤦 https://github.com/pytorch/pytorch/blob/ad81eeb7c7c906e0cdd04a5cc8fdb9592281c317/aten/src/ATen/native/cuda/linalg/BatchLinearAlgebraLib.cpp#L445-L448

So something like this works, and probably needs to be applied to a couple of other APIs too:

diff --git a/lib/cusolver/dense_generic.jl b/lib/cusolver/dense_generic.jl
index 3dbb48874..0a5222437 100644
--- a/lib/cusolver/dense_generic.jl
+++ b/lib/cusolver/dense_generic.jl
@@ -287,28 +287,24 @@ function Xgesvdp!(jobz::Char, econ::Int, A::StridedCuMatrix{T}) where {T <: Blas
     p = min(m, n)
     R = real(T)
     econ ∈ (0, 1) || throw(ArgumentError("econ is incorrect. The values accepted are 0 and 1."))
-    U = if jobz == 'V' && econ == 1
+    U = if econ == 1
         CuMatrix{T}(undef, m, p)
-    elseif jobz == 'V' && econ == 0
+    elseif econ == 0
         CuMatrix{T}(undef, m, m)
-    elseif jobz == 'N'
-        CU_NULL
     else
         throw(ArgumentError("jobz is incorrect. The values accepted are 'V' and 'N'."))
     end
     Σ = CuVector{R}(undef, p)
-    V = if jobz == 'V' && econ == 1
+    V = if econ == 1
         CuMatrix{T}(undef, n, p)
-    elseif jobz == 'V' && econ == 0
+    elseif econ == 0
         CuMatrix{T}(undef, n, n)
-    elseif jobz == 'N'
-        CU_NULL
     else
         throw(ArgumentError("jobz is incorrect. The values accepted are 'V' and 'N'."))
     end
     lda = max(1, stride(A, 2))
-    ldu = U == CU_NULL ? 1 : max(1, stride(U, 2))
-    ldv = V == CU_NULL ? 1 : max(1, stride(V, 2))
+    ldu = U == CU_NULL ? m : max(1, stride(U, 2))
+    ldv = V == CU_NULL ? n : max(1, stride(V, 2))
     h_err_sigma = Ref{Cdouble}(0)
     params = CuSolverParameters()
     dh = dense_handle()
@@ -331,6 +327,8 @@ function Xgesvdp!(jobz::Char, econ::Int, A::StridedCuMatrix{T}) where {T <: Blas
     flag = @allowscalar dh.info[1]
     chklapackerror(flag |> BlasInt)
     if jobz == 'N'
+        unsafe_free!(U)
+        unsafe_free!(V)
         return Σ, h_err_sigma[]
     elseif jobz == 'V'
         return U, Σ, V, h_err_sigma[]

cc @kshyatt

[Jutho]

I should have read the CUSOLVER manual better. It does indeed not say that it doesn't need U and V matrices if only requesting singular values. In fact, it says U and V always need to be at least m x m and n x n, even if using the economy mode where you only want the first min(m, n) columns. So it also uses these as workspace during the computation I assume.

[amontoison]

It's probably my fault, for the normal SVD we can avoid some allocations and I did the same kind of modifications for this routine.

I can fix it the week after next one.

[amontoison]

@maleadt @Jutho
I confirm that something is wrong in CUSOLVER. If I only provide a CuMatrix with one column for U / V and jobz = 'N', it works and return the correct result (except for ComplexF64).
In theory U and V should be CU_NULL and any additional GPU memory needed should in the GPU buffer.

I opened #2763 to temporary fix the issue.
However, I will open a ticket on the website of NVIDIA.