[CUSOLVER] Interface larft!

Author: amontoison

lib/cusolver/base.jl

@@ -43,3 +43,23 @@ function Base.convert(::Type{cusolverEigRange_t}, range::Char)
         throw(ArgumentError("Unknown eigenvalue solver range $range."))
     end
 end
+
+function Base.convert(::Type{cusolverStorevMode_t}, storev::Char)
+    if storev == 'C'
+        CUBLAS_STOREV_COLUMNWISE
+    elseif storev == 'R'
+        CUBLAS_STOREV_ROWWISE
+    else
+        throw(ArgumentError("Unknown storage mode $storev."))
+    end
+end
+
+function Base.convert(::Type{cusolverDirectMode_t}, direct::Char)
+    if direct == 'F'
+        CUBLAS_DIRECT_FORWARD
+    elseif direct == 'B'
+        CUBLAS_DIRECT_BACKWARD
+    else
+        throw(ArgumentError("Unknown direction mode $direct."))
+    end
+end

lib/cusolver/dense_generic.jl

@@ -192,6 +192,33 @@ function trtri!(uplo::Char, diag::Char, A::StridedCuMatrix{T}) where {T <: BlasF
     A
 end
 
+# Xlarft!
+function larft!(direct::Char, v::StridedCuMatrix{T}, tau::StridedCuVector{T}, t::StridedCuMatrix{T}) where {T <: BlasFloat}
+    n, k = size(v)
+    ktau = length(tau)
+    mt, nt = size(t)
+    (n < k) && throw(ArgumentError("The number of elementary reflectors ($k) must be lower or equal to the order of block reflector H ($n)."))
+    (ktau != k) && throw(ArgumentError("The length of tau ($ktau) is not equal to the number of elementary reflectors ($k)."))
+    (mt != k || nt != k) && throw(ArgumentError("The size of the triangular factor of the block reflector is ($mt, $nt) and must be ($k, $k)."))
+    ldv = max(1, stride(v, 2))
+    ldt = max(1, stride(t, 2))
+    params = CuSolverParameters()
+
+    function bufferSize()
+        out_cpu = Ref{Csize_t}(0)
+        out_gpu = Ref{Csize_t}(0)
+        cusolverDnXlarft_bufferSize(dense_handle(), params, direct, 'C', n, k, T,
+                                    v, ldv, T, tau, T, t, ldt, T, out_gpu, out_cpu)
+        out_gpu[], out_cpu[]
+    end
+    with_workspaces(bufferSize()...) do buffer_gpu, buffer_cpu
+        cusolverDnXlarft(dense_handle(), params, direct, 'C', n, k, T, v, ldv, T, tau, T, t,
+                         ldt, T, buffer_gpu, sizeof(buffer_gpu), buffer_cpu, sizeof(buffer_cpu))
+    end
+
+    t
+end
+
 # Xgesvd
 function Xgesvd!(jobu::Char, jobvt::Char, A::StridedCuMatrix{T}) where {T <: BlasFloat}
     m, n = size(A)

lib/cusolver/libcusolver.jl

@@ -4890,8 +4890,8 @@ end
                                                dataTypeTau::cudaDataType,
                                                d_tau::CuPtr{Cvoid}, dataTypeT::cudaDataType,
                                                d_T::CuPtr{Cvoid}, ldt::Int64,
-                                               computeType::CuPtr{Cvoid},
-                                               bufferOnDevice::Ptr{Cvoid},
+                                               computeType::cudaDataType,
+                                               bufferOnDevice::CuPtr{Cvoid},
                                                workspaceInBytesOnDevice::Csize_t,
                                                bufferOnHost::Ptr{Cvoid},
                                                workspaceInBytesOnHost::Csize_t)::cusolverStatus_t

res/wrap/cusolver.toml

@@ -2338,4 +2338,4 @@ needs_context = false
 8 = "CuPtr{Cvoid}"
 11 = "CuPtr{Cvoid}"
 13 = "CuPtr{Cvoid}"
-15 = "CuPtr{Cvoid}"
+16 = "CuPtr{Cvoid}"

test/libraries/cusolver/dense_generic.jl

@@ -5,7 +5,33 @@ m = 15
 n = 10
 p = 5
 
-@testset "cusolver -- generic API -- $elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
+@testset "cusolver -- generic API -- $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
+    @testset "larft!" begin
+        @testset "direct = $direct" for direct in ('F', 'B')
+            direct == 'B' && continue
+            A = rand(elty,m,n)
+            t = rand(elty,n,n)
+
+            dA = CuMatrix(A)
+            dA, dτ = CUSOLVER.geqrf!(dA)
+            hI = Matrix{elty}(I, m, m)
+            dI = CuArray(hI)
+            dH = CUSOLVER.ormqr!('L', 'N', dA, dτ, copy(dI))
+
+            v = Array(dA)
+            for j = 1:n
+                v[j,j] = one(elty)
+                for i = 1:j-1
+                    v[i,j] = zero(elty)
+                end
+            end
+            dv = CuArray(v)
+            dt = CuMatrix(t)
+            dt = CUSOLVER.larft!(direct, dv, dτ, dt)
+            @test dI - dv * dt * dv' ≈ dH
+        end
+    end
+
     @testset "sytrs!" begin
         for uplo in ('L', 'U')
             A = rand(elty,n,n)