Add hpmv! to BLAS in stdlib/LinearAlgebra (#34211)

Author: iolaszlo

NEWS.md

@@ -31,6 +31,7 @@ Standard library changes
 
 #### LinearAlgebra
 
+* The BLAS submodule now supports the level-2 BLAS subroutine `hpmv!` ([#34211]).
 
 #### Markdown
 

stdlib/LinearAlgebra/src/blas.jl

@@ -28,6 +28,7 @@ export
     gemv,
     hemv!,
     hemv,
+    hpmv!,
     sbmv!,
     sbmv,
     symv!,
@@ -823,6 +824,82 @@ Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
 """
 hemv(ul, A, x)
 
+### hpmv!, (HP) Hermitian packed matrix-vector operation defined as y := alpha*A*x + beta*y.
+for (fname, elty) in ((:zhpmv_, :ComplexF64),
+                      (:chpmv_, :ComplexF32))
+    @eval begin
+        # SUBROUTINE ZHPMV(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 hpmv!(uplo::AbstractChar,
+                       n::BlasInt,
+                       α::$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, output
+                   Ref{BlasInt}),  # incy
+                  uplo,
+                  n,
+                  α,
+                  AP,
+                  x,
+                  incx,
+                  β,
+                  y,
+                  incy)
+        end
+    end
+end
+
+function hpmv!(uplo::AbstractChar,
+               α::Number, AP::Union{DenseArray{T}, AbstractVector{T}}, x::Union{DenseArray{T}, AbstractVector{T}},
+               β::Number, y::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasComplex}
+    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 length(AP) < Int64(N*(N+1)/2)
+        throw(DimensionMismatch("Packed Hermitian matrix A has size smaller than length(x) =  $(N)."))
+    end
+    GC.@preserve x y AP hpmv!(uplo, N, convert(T, α), AP, pointer(x), BlasInt(stride(x, 1)), convert(T, β), pointer(y), BlasInt(stride(y, 1)))
+    y
+end
+
+"""
+    hpmv!(uplo, α, AP, x, β, y)
+
+Update vector `y` as `α*AP*x + β*y` where `AP` is a packed Hermitian matrix.
+The storage layout for `AP` is described in the reference BLAS module, level-2 BLAS at
+<http://www.netlib.org/lapack/explore-html/>.
+
+The scalar inputs `α` and `β` shall be numbers.
+
+The array inputs `x`, `y` and `AP` must be complex one-dimensional julia arrays of the
+same type that is either `ComplexF32` or `ComplexF64`.
+
+Return the updated `y`.
+"""
+hpmv!
+
 ### sbmv, (SB) symmetric banded matrix-vector multiplication
 for (fname, elty) in ((:dsbmv_,:Float64),
                       (:ssbmv_,:Float32))

stdlib/LinearAlgebra/test/blas.jl

@@ -8,6 +8,7 @@ using LinearAlgebra: BlasReal, BlasComplex
 Random.seed!(100)
 ## BLAS tests - testing the interface code to BLAS routines
 @testset for elty in [Float32, Float64, ComplexF32, ComplexF64]
+
     @testset "syr2k!" begin
         U = randn(5,2)
         V = randn(5,2)
@@ -200,6 +201,46 @@ Random.seed!(100)
             @test_throws DimensionMismatch BLAS.trmm('R','U','N','N',one(elty),triu(Cnn),Cnm)
         end
 
+        # hpmv!
+        if elty in (ComplexF32, ComplexF64)
+            @testset "hpmv!" begin
+                # Both matrix dimensions n coincide, as we have Hermitian matrices.
+                # Define the inputs and outputs of hpmv!, y = α*A*x+β*y
+                α = rand(elty)
+                M = rand(elty, n, n)
+                A = (M+M')/elty(2.0)
+                x = rand(elty, n)
+                β = rand(elty)
+                y = rand(elty, n)
+
+                y_result_julia = α*A*x+β*y
+
+                # Create lower triangular packing of A
+                AP = typeof(A[1,1])[]
+                for j in 1:n
+                    for i in j:n
+                        push!(AP, A[i,j])
+                    end
+                end
+
+                y_result_blas_lower = copy(y)
+                BLAS.hpmv!('L', α, AP, x, β, y_result_blas_lower)
+                @test y_result_julia≈y_result_blas_lower
+
+                # Create upper triangular packing of A
+                AP = typeof(A[1,1])[]
+                for j in 1:n
+                    for i in 1:j
+                        push!(AP, A[i,j])
+                    end
+                end
+
+                y_result_blas_upper = copy(y)
+                BLAS.hpmv!('U', α, AP, x, β, y_result_blas_upper)
+                @test y_result_julia≈y_result_blas_upper
+            end
+        end
+
         #trsm
         A = triu(rand(elty,n,n))
         B = rand(elty,(n,n))