Wrap CUBLAS.spmv and spr (#1248)

Author: bjarthur

lib/cublas/wrappers.jl

@@ -415,6 +415,32 @@ for (fname, elty) in ((:cublasDgbmv_v2,:Float64),
     end
 end
 
+### spmv
+for (fname, elty) in ((:cublasDspmv_v2,:Float64),
+                      (:cublasSspmv_v2,:Float32))
+    @eval begin
+        function spmv!(uplo::Char,
+                       alpha::Number,
+                       AP::StridedCuVector{$elty},
+                       x::StridedCuVector{$elty},
+                       beta::Number,
+                       y::StridedCuVector{$elty})
+            n = round(Int, (sqrt(8*length(AP))-1)/2)
+            if n != length(x) || n != length(y) throw(DimensionMismatch("")) end
+            incx = stride(x,1)
+            incy = stride(y,1)
+            $fname(handle(), uplo, n, alpha, AP, x, incx, beta, y, incy)
+            y
+        end
+        function spmv(uplo::Char, alpha::Number, AP::StridedCuVector{$elty}, x::StridedCuVector{$elty})
+                spmv!(uplo, alpha, AP, x, zero($elty), similar(x))
+        end
+        function spmv(uplo::Char, AP::StridedCuVector{$elty}, x::StridedCuVector{$elty})
+            spmv(uplo, one($elty), AP, x)
+        end
+    end
+end
+
 ### symv
 for (fname, elty) in ((:cublasDsymv_v2,:Float64),
                       (:cublasSsymv_v2,:Float32),
@@ -695,6 +721,23 @@ for (fname, elty) in ((:cublasDger_v2,:Float64),
     end
 end
 
+### spr
+for (fname, elty) in ((:cublasDspr_v2,:Float64),
+                      (:cublasSspr_v2,:Float32))
+    @eval begin
+        function spr!(uplo::Char,
+                      alpha::Number,
+                      x::StridedCuVector{$elty},
+                      AP::StridedCuVector{$elty})
+            n = round(Int, (sqrt(8*length(AP))-1)/2)
+            length(x) == n || throw(DimensionMismatch("Length of vector must be the same as the matrix dimensions"))
+            incx = stride(x,1)
+            $fname(handle(), uplo, n, alpha, x, incx, AP)
+            AP
+        end
+    end
+end
+
 ### syr
 # TODO: check calls in julia b/c blas may not define syr for Z and C
 for (fname, elty) in ((:cublasDsyr_v2,:Float64),

test/cublas.jl

@@ -298,6 +298,69 @@ end
         dhA = CuArray(hA)
         x = rand(elty,m)
         dx = CuArray(x)
+
+        function pack(A, uplo)
+            AP = Vector{elty}(undef, (n*(n+1))>>1)
+            k = 1
+            for j in 1:n
+                for i in (uplo==:L ? (j:n) : (1:j))
+                    AP[k] = A[i,j]
+                    k += 1
+                end
+            end
+            return AP
+        end
+
+        if elty in ["Float32", "Float64"]
+            # pack matrices
+            sAPU = pack(sA, :U)
+            dsAPU = CuVector(sAPU)
+            sAPL = pack(sA, :L)
+            dsAPL = CuVector(sAPL)
+
+            @testset "spmv!" begin
+                # generate vectors
+                y = rand(elty,m)
+                # copy to device
+                dy = CuArray(y)
+                # execute on host
+                BLAS.spmv!('U',alpha,sAPU,x,beta,y)
+                # execute on device
+                CUBLAS.spmv!('U',alpha,dsAPU,dx,beta,dy)
+                # compare results
+                hy = Array(dy)
+                @test y ≈ hy
+                # execute on host
+                BLAS.spmv!('U',alpha,sAPL,x,beta,y)
+                # execute on device
+                CUBLAS.spmv!('U',alpha,dsAPL,dx,beta,dy)
+                # compare results
+                hy = Array(dy)
+                @test y ≈ hy
+            end
+
+            @testset "spr!" begin
+                # execute on host
+                VERSION>=v"1.8.0-DEV.1049" && BLAS.spr!('U',alpha,x,sAPU)
+                # execute on device
+                CUBLAS.spr!('U',alpha,dx,dsAPU)
+                # compare results
+                if VERSION>=v"1.8.0-DEV.1049"
+                    hsAPU = Array(dsAPU)
+                    @test sAPU ≈ hsAPU
+                end
+                # execute on host
+                VERSION>=v"1.8.0-DEV.1049" && BLAS.spr!('U',alpha,x,sAPL)
+                # execute on device
+                CUBLAS.spr!('U',alpha,dx,dsAPL)
+                # compare results
+                if VERSION>=v"1.8.0-DEV.1049"
+                    hAPL = Array(dAPL)
+                    @test sAPL ≈ hAPL
+                end
+            end
+        end
+
         @testset "symv!" begin
             # generate vectors
             y = rand(elty,m)