Skip to content

Commit

Permalink
Minor updates to BLAS/hpmv! as in JuliaLang#34320. [JuliaLang#34439]
Browse files Browse the repository at this point in the history
Follow the syntax, docstring and testing conventions
elaborated in JuliaLang#34320.
  • Loading branch information
iolaszlo committed Jan 19, 2020
1 parent e6828c8 commit cb11ac7
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 18 deletions.
24 changes: 16 additions & 8 deletions stdlib/LinearAlgebra/src/blas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -838,7 +838,7 @@ for (fname, elty) in ((:zhpmv_, :ComplexF64),
# * .. Array Arguments ..
# DOUBLE PRECISION A(N,N),X(N),Y(N)
function hpmv!(uplo::AbstractChar,
n::BlasInt,
n::Integer,
α::$elty,
AP::Union{Ptr{$elty}, AbstractArray{$elty}},
x::Union{Ptr{$elty}, AbstractArray{$elty}},
Expand Down Expand Up @@ -881,21 +881,29 @@ function hpmv!(uplo::AbstractChar,
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, BlasInt(N), convert(T, α), AP, pointer(x), BlasInt(stride(x, 1)), convert(T, β), pointer(y), BlasInt(stride(y, 1)))
GC.@preserve x y AP hpmv!(uplo, N, convert(T, α), AP, pointer(x), stride(x, 1), convert(T, β), pointer(y), 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/>.
Update vector `y` as `α*A*x + β*y`, where `A` is a Hermitian matrix provided
in packed format `AP`.
With `uplo = 'U'`, the array AP must contain the upper triangular part of the
Hermitian 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
Hermitian 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 `β` shall be numbers.
The scalar inputs `α` and `β` must be complex or real 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`.
The array inputs `x`, `y` and `AP` must all be of `ComplexF32` or `ComplexF64` type.
Return the updated `y`.
"""
Expand Down
23 changes: 13 additions & 10 deletions stdlib/LinearAlgebra/test/blas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -208,36 +208,39 @@ Random.seed!(100)
# Define the inputs and outputs of hpmv!, y = α*A*x+β*y
α = rand(elty)
M = rand(elty, n, n)
A = (M+M')/elty(2.0)
AL = Hermitian(M, :L)
AU = Hermitian(M, :U)
x = rand(elty, n)
β = rand(elty)
y = rand(elty, n)

y_result_julia = α*A*x+β*y
y_result_julia_lower = α*AL*x+β*y

# Create lower triangular packing of A
AP = typeof(A[1,1])[]
# Create lower triangular packing of AL
AP = typeof(AL[1,1])[]
for j in 1:n
for i in j:n
push!(AP, A[i,j])
push!(AP, AL[i,j])
end
end

y_result_blas_lower = copy(y)
BLAS.hpmv!('L', α, AP, x, β, y_result_blas_lower)
@test y_result_juliay_result_blas_lower
@test y_result_julia_lowery_result_blas_lower

# Create upper triangular packing of A
AP = typeof(A[1,1])[]
y_result_julia_upper = α*AU*x+β*y

# Create upper triangular packing of AU
AP = typeof(AU[1,1])[]
for j in 1:n
for i in 1:j
push!(AP, A[i,j])
push!(AP, AU[i,j])
end
end

y_result_blas_upper = copy(y)
BLAS.hpmv!('U', α, AP, x, β, y_result_blas_upper)
@test y_result_juliay_result_blas_upper
@test y_result_julia_uppery_result_blas_upper
end
end

Expand Down

0 comments on commit cb11ac7

Please sign in to comment.