Skip to content

Commit 6a9a926

Browse files
authored
Blas hpmv minor updates (#34439)
1 parent e803d78 commit 6a9a926

File tree

2 files changed

+29
-17
lines changed

2 files changed

+29
-17
lines changed

stdlib/LinearAlgebra/src/blas.jl

Lines changed: 16 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -870,6 +870,7 @@ for (fname, elty) in ((:zhpmv_, :ComplexF64),
870870
β,
871871
y,
872872
incy)
873+
return y
873874
end
874875
end
875876
end
@@ -882,7 +883,7 @@ function hpmv!(uplo::AbstractChar,
882883
if N != length(y)
883884
throw(DimensionMismatch("x has length $(N), but y has length $(length(y))"))
884885
end
885-
if length(AP) < Int64(N*(N+1)/2)
886+
if 2*length(AP) < N*(N + 1)
886887
throw(DimensionMismatch("Packed Hermitian matrix A has size smaller than length(x) = $(N)."))
887888
end
888889
return hpmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1))
@@ -891,14 +892,22 @@ end
891892
"""
892893
hpmv!(uplo, α, AP, x, β, y)
893894
894-
Update vector `y` as `α*AP*x + β*y` where `AP` is a packed Hermitian matrix.
895-
The storage layout for `AP` is described in the reference BLAS module, level-2 BLAS at
896-
<http://www.netlib.org/lapack/explore-html/>.
895+
Update vector `y` as `α*A*x + β*y`, where `A` is a Hermitian matrix provided
896+
in packed format `AP`.
897+
898+
With `uplo = 'U'`, the array AP must contain the upper triangular part of the
899+
Hermitian matrix packed sequentially, column by column, so that `AP[1]`
900+
contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[1, 2]` and `A[2, 2]`
901+
respectively, and so on.
902+
903+
With `uplo = 'L'`, the array AP must contain the lower triangular part of the
904+
Hermitian matrix packed sequentially, column by column, so that `AP[1]`
905+
contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[2, 1]` and `A[3, 1]`
906+
respectively, and so on.
897907
898-
The scalar inputs `α` and `β` shall be numbers.
908+
The scalar inputs `α` and `β` must be complex or real numbers.
899909
900-
The array inputs `x`, `y` and `AP` must be complex one-dimensional julia arrays of the
901-
same type that is either `ComplexF32` or `ComplexF64`.
910+
The array inputs `x`, `y` and `AP` must all be of `ComplexF32` or `ComplexF64` type.
902911
903912
Return the updated `y`.
904913
"""

stdlib/LinearAlgebra/test/blas.jl

Lines changed: 13 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -208,36 +208,39 @@ Random.seed!(100)
208208
# Define the inputs and outputs of hpmv!, y = α*A*x+β*y
209209
α = rand(elty)
210210
M = rand(elty, n, n)
211-
A = (M+M')/elty(2.0)
211+
AL = Hermitian(M, :L)
212+
AU = Hermitian(M, :U)
212213
x = rand(elty, n)
213214
β = rand(elty)
214215
y = rand(elty, n)
215216

216-
y_result_julia = α*A*x+β*y
217+
y_result_julia_lower = α*AL*x + β*y
217218

218-
# Create lower triangular packing of A
219-
AP = typeof(A[1,1])[]
219+
# Create lower triangular packing of AL
220+
AP = typeof(AL[1,1])[]
220221
for j in 1:n
221222
for i in j:n
222-
push!(AP, A[i,j])
223+
push!(AP, AL[i,j])
223224
end
224225
end
225226

226227
y_result_blas_lower = copy(y)
227228
BLAS.hpmv!('L', α, AP, x, β, y_result_blas_lower)
228-
@test y_result_juliay_result_blas_lower
229+
@test y_result_julia_lower y_result_blas_lower
229230

230-
# Create upper triangular packing of A
231-
AP = typeof(A[1,1])[]
231+
y_result_julia_upper = α*AU*x + β*y
232+
233+
# Create upper triangular packing of AU
234+
AP = typeof(AU[1,1])[]
232235
for j in 1:n
233236
for i in 1:j
234-
push!(AP, A[i,j])
237+
push!(AP, AU[i,j])
235238
end
236239
end
237240

238241
y_result_blas_upper = copy(y)
239242
BLAS.hpmv!('U', α, AP, x, β, y_result_blas_upper)
240-
@test y_result_juliay_result_blas_upper
243+
@test y_result_julia_upper y_result_blas_upper
241244
end
242245
end
243246

0 commit comments

Comments
 (0)