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

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:https://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))