Support sparse opnorm (#1466)

Author: Roger-luo

lib/cusparse/CUSPARSE.jl

@@ -27,6 +27,7 @@ include("libcusparse_deprecated.jl")
 include("array.jl")
 include("util.jl")
 include("types.jl")
+include("linalg.jl")
 
 # low-level wrappers
 include("helpers.jl")

lib/cusparse/linalg.jl

@@ -0,0 +1,36 @@
+using LinearAlgebra
+
+function sum_dim1(A::CuSparseMatrixCSR)
+    function kernel(Tnorm, out, dA)
+        idx = (blockIdx().x-1) * blockDim().x + threadIdx().x
+        idx < length(dA.rowPtr) || return
+        s = zero(Tnorm)
+        for k in dA.rowPtr[idx]:dA.rowPtr[idx+1]-1
+            s += abs(dA.nzVal[k])
+        end
+        out[idx] = s
+        return
+    end
+
+    m, n = size(A)
+    Tnorm = typeof(float(real(zero(eltype(A)))))
+    Tsum = promote_type(Float64,Tnorm)
+    rowsum = CUDA.CuArray{Tsum}(undef, m)
+    kernel_f = @cuda launch=false kernel(Tnorm, rowsum, A)
+    
+    config = launch_configuration(kernel_f.fun)
+    threads = min(n, config.threads)
+    blocks = cld(n, threads)
+    kernel_f(Tnorm, rowsum, A; threads, blocks)
+    return rowsum
+end
+
+function LinearAlgebra.opnorm(A::CuSparseMatrixCSR, p::Real=2)
+    if p == Inf
+        return maximum(sum_dim1(A))
+    else
+        error("p=$p is not supported")
+    end
+end
+
+LinearAlgebra.opnorm(A::CuSparseMatrixCSC, p::Real=2) = opnorm(CuSparseMatrixCSR(A), p)

test/cusparse/linalg.jl

@@ -0,0 +1,10 @@
+using CUDA.CUSPARSE
+using LinearAlgebra, SparseArrays
+
+@testset "opnorm" for T in [Float32, Float64, ComplexF32, ComplexF64]
+    S = sprand(T, 10, 10, 0.1)
+    dS_csc = CuSparseMatrixCSC(S)
+    dS_csr = CuSparseMatrixCSR(S)
+    @test opnorm(S, Inf) ≈ opnorm(dS_csc, Inf)
+    @test opnorm(S, Inf) ≈ opnorm(dS_csr, Inf)
+end