add spr! to LinearAlgebra.BLAS (#42830)

Author: bjarthur

NEWS.md

@@ -89,6 +89,7 @@ Standard library changes
 #### Package Manager
 
 #### LinearAlgebra
+* The BLAS submodule now supports the level-2 BLAS subroutine `spr!` ([#42830]).
 
 #### Markdown
 

stdlib/LinearAlgebra/src/blas.jl

@@ -33,6 +33,7 @@ export
     sbmv!,
     sbmv,
     spmv!,
+    spr!,
     symv!,
     symv,
     trsv!,
@@ -1145,6 +1146,71 @@ Return the updated `y`.
 """
 spmv!
 
+### spr!, (SP) symmetric packed matrix-vector operation defined as A := alpha*x*x' + A
+for (fname, elty) in ((:dspr_, :Float64),
+                      (:sspr_, :Float32))
+    @eval begin
+        function spr!(uplo::AbstractChar,
+                      n::Integer,
+                      α::$elty,
+                      x::Union{Ptr{$elty}, AbstractArray{$elty}},
+                      incx::Integer,
+                      AP::Union{Ptr{$elty}, AbstractArray{$elty}})
+
+            ccall((@blasfunc($fname), libblastrampoline), Cvoid,
+                  (Ref{UInt8},     # uplo,
+                   Ref{BlasInt},   # n,
+                   Ref{$elty},     # α,
+                   Ptr{$elty},     # x,
+                   Ref{BlasInt},   # incx,
+                   Ptr{$elty},     # AP,
+                   Clong),         # length of uplo
+                  uplo,
+                  n,
+                  α,
+                  x,
+                  incx,
+                  AP,
+                  1)
+            return AP
+        end
+    end
+end
+
+function spr!(uplo::AbstractChar,
+              α::Real, x::Union{DenseArray{T}, AbstractVector{T}},
+              AP::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasReal}
+    require_one_based_indexing(AP, x)
+    N = length(x)
+    if 2*length(AP) < N*(N + 1)
+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    return spr!(uplo, N, convert(T, α), x, stride(x, 1), AP)
+end
+
+"""
+    spr!(uplo, α, x, AP)
+
+Update matrix `A` as `α*A*x*x'`, where `A` is a symmetric matrix provided
+in packed format `AP` and `x` is a vector.
+
+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 input `α` must be real.
+
+The array inputs `x` and `AP` must all be of `Float32` or `Float64` type.
+Return the updated `AP`.
+"""
+spr!
+
 ### hbmv, (HB) Hermitian banded matrix-vector multiplication
 for (fname, elty) in ((:zhbmv_,:ComplexF64),
                       (:chbmv_,:ComplexF32))

stdlib/LinearAlgebra/test/blas.jl

@@ -314,6 +314,49 @@ Random.seed!(100)
             end
         end
 
+        # spr!
+        if elty in (Float32, Float64)
+            @testset "spr! $elty" begin
+                α = rand(elty)
+                M = rand(elty, n, n)
+                AL = Symmetric(M, :L)
+                AU = Symmetric(M, :U)
+                x = rand(elty, n)
+
+                function pack(A, uplo)
+                    AP = elty[]
+                    for j in 1:n
+                        for i in (uplo==:L ? (j:n) : (1:j))
+                            push!(AP, A[i,j])
+                        end
+                    end
+                    return AP
+                end
+
+                ALP_result_julia_lower = pack(α*x*x' + AL, :L)
+                ALP_result_blas_lower = pack(AL, :L)
+                BLAS.spr!('L', α, x, ALP_result_blas_lower)
+                @test ALP_result_julia_lower ≈ ALP_result_blas_lower
+                ALP_result_blas_lower = append!(pack(AL, :L), ones(elty, 10))
+                BLAS.spr!('L', α, x, ALP_result_blas_lower)
+                @test ALP_result_julia_lower ≈ ALP_result_blas_lower[1:end-10]
+                ALP_result_blas_lower = reshape(pack(AL, :L), 1, length(ALP_result_julia_lower), 1)
+                BLAS.spr!('L', α, x, ALP_result_blas_lower)
+                @test ALP_result_julia_lower ≈ vec(ALP_result_blas_lower)
+
+                AUP_result_julia_upper = pack(α*x*x' + AU, :U)
+                AUP_result_blas_upper = pack(AU, :U)
+                BLAS.spr!('U', α, x, AUP_result_blas_upper)
+                @test AUP_result_julia_upper ≈ AUP_result_blas_upper
+                AUP_result_blas_upper = append!(pack(AU, :U), ones(elty, 10))
+                BLAS.spr!('U', α, x, AUP_result_blas_upper)
+                @test AUP_result_julia_upper ≈ AUP_result_blas_upper[1:end-10]
+                AUP_result_blas_upper = reshape(pack(AU, :U), 1, length(AUP_result_julia_upper), 1)
+                BLAS.spr!('U', α, x, AUP_result_blas_upper)
+                @test AUP_result_julia_upper ≈ vec(AUP_result_blas_upper)
+            end
+        end
+
         #trsm
         A = triu(rand(elty,n,n))
         B = rand(elty,(n,n))