Add spmv! to BLAS in stdlib/LinearAlgebra (#34320)

NEWS.md

@@ -79,6 +79,7 @@ Standard library changes
 * The BLAS submodule now supports the level-2 BLAS subroutine `hpmv!` ([#34211]).
 * `normalize` now supports multidimensional arrays ([#34239])
 * `lq` factorizations can now be used to compute the minimum-norm solution to under-determined systems ([#34350]).
+* The BLAS submodule now supports the level-2 BLAS subroutine `spmv!` ([#34320]).
 
 #### Markdown
 

stdlib/LinearAlgebra/src/blas.jl

@@ -31,6 +31,7 @@ export
     hpmv!,
     sbmv!,
     sbmv,
+    spmv!,
     symv!,
     symv,
     trsv!,
@@ -977,6 +978,90 @@ Return the updated `y`.
 """
 sbmv!
 
+### spmv!, (SP) symmetric packed matrix-vector operation defined as y := alpha*A*x + beta*y.
+for (fname, elty) in ((:dspmv_, :Float64),
+                      (:sspmv_, :Float32))
+    @eval begin
+        # SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
+        # Y <- ALPHA*AP*X + BETA*Y
+        # *     .. Scalar Arguments ..
+        #       DOUBLE PRECISION ALPHA,BETA
+        #       INTEGER INCX,INCY,N
+        #       CHARACTER UPLO
+        # *     .. Array Arguments ..
+        #       DOUBLE PRECISION A(N,N),X(N),Y(N)
+        function spmv!(uplo::AbstractChar,
+                       n::Integer,
+                       α::$elty,
+                       AP::Union{Ptr{$elty}, AbstractArray{$elty}},
+                       x::Union{Ptr{$elty}, AbstractArray{$elty}},
+                       incx::Integer,
+                       β::$elty,
+                       y::Union{Ptr{$elty}, AbstractArray{$elty}},
+                       incy::Integer)
+
+            ccall((@blasfunc($fname), libblas), Cvoid,
+                  (Ref{UInt8},     # uplo,
+                   Ref{BlasInt},   # n,
+                   Ref{$elty},     # α,
+                   Ptr{$elty},     # AP,
+                   Ptr{$elty},     # x,
+                   Ref{BlasInt},   # incx,
+                   Ref{$elty},     # β,
+                   Ptr{$elty},     # y, out
+                   Ref{BlasInt}),  # incy
+                  uplo,
+                  n,
+                  α,
+                  AP,
+                  x,
+                  incx,
+                  β,
+                  y,
+                  incy)
+            return y
+        end
+    end
+end
+
+function spmv!(uplo::AbstractChar,
+               α::Real, AP::Union{DenseArray{T}, AbstractVector{T}}, x::Union{DenseArray{T}, AbstractVector{T}},
+               β::Real, y::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasReal}
+    require_one_based_indexing(AP, x, y)
+    N = length(x)
+    if N != length(y)
+        throw(DimensionMismatch("x has length $(N), but y has length $(length(y))"))
+    end
+    if 2*length(AP) < N*(N + 1)
+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    return spmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1))
+end
+
+"""
+    spmv!(uplo, α, AP, x, β, y)
+
+Update vector `y` as `α*A*x + β*y`, where `A` is a symmetric matrix provided
+in packed format `AP`.
+
+With `uplo = 'U'`, the array AP must contain the upper triangular part of the
+symmetric matrix packed sequentially, column by column, so that `AP[1]`
+contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[1, 2]` and `A[2, 2]`
+respectively, and so on.
+
+With `uplo = 'L'`, the array AP must contain the lower triangular part of the
+symmetric matrix packed sequentially, column by column, so that `AP[1]`
+contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[2, 1]` and `A[3, 1]`
+respectively, and so on.
+
+The scalar inputs `α` and `β` must be real.
+
+The array inputs `x`, `y` and `AP` must all be of `Float32` or `Float64` type.
+
+Return the updated `y`.
+"""
+spmv!
+
 ### hbmv, (HB) Hermitian banded matrix-vector multiplication
 for (fname, elty) in ((:zhbmv_,:ComplexF64),
                       (:chbmv_,:ComplexF32))

stdlib/LinearAlgebra/test/blas.jl

@@ -244,6 +244,50 @@ Random.seed!(100)
             end
         end
 
+        # spmv!
+        if elty in (Float32, Float64)
+            @testset "spmv!" begin
+                # Both matrix dimensions n coincide, as we have symmetric matrices.
+                # Define the inputs and outputs of spmv!, y = α*A*x+β*y
+                α = rand(elty)
+                M = rand(elty, n, n)
+                AL = Symmetric(M, :L)
+                AU = Symmetric(M, :U)
+                x = rand(elty, n)
+                β = rand(elty)
+                y = rand(elty, n)
+
+                y_result_julia_lower = α*AL*x + β*y
+
+                # Create lower triangular packing of AL
+                AP = typeof(M[1,1])[]
+                for j in 1:n
+                    for i in j:n
+                        push!(AP, AL[i,j])
+                    end
+                end
+
+                y_result_blas_lower = copy(y)
+                BLAS.spmv!('L', α, AP, x, β, y_result_blas_lower)
+                @test y_result_julia_lower ≈ y_result_blas_lower
+
+
+                y_result_julia_upper = α*AU*x + β*y
+
+                # Create upper triangular packing of AU
+                AP = typeof(M[1,1])[]
+                for j in 1:n
+                    for i in 1:j
+                        push!(AP, AU[i,j])
+                    end
+                end
+
+                y_result_blas_upper = copy(y)
+                BLAS.spmv!('U', α, AP, x, β, y_result_blas_upper)
+                @test y_result_julia_upper ≈ y_result_blas_upper
+            end
+        end
+
         #trsm
         A = triu(rand(elty,n,n))
         B = rand(elty,(n,n))