diff --git a/stdlib/LinearAlgebra/src/blas.jl b/stdlib/LinearAlgebra/src/blas.jl index 9124161cfb3ad..052c227fa1861 100644 --- a/stdlib/LinearAlgebra/src/blas.jl +++ b/stdlib/LinearAlgebra/src/blas.jl @@ -870,6 +870,7 @@ for (fname, elty) in ((:zhpmv_, :ComplexF64), β, y, incy) + return y end end end @@ -882,7 +883,7 @@ function hpmv!(uplo::AbstractChar, 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) + if 2*length(AP) < N*(N + 1) throw(DimensionMismatch("Packed Hermitian matrix A has size smaller than length(x) = $(N).")) end return hpmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1)) @@ -891,14 +892,22 @@ 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 -. +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`. """ diff --git a/stdlib/LinearAlgebra/test/blas.jl b/stdlib/LinearAlgebra/test/blas.jl index 4a21fbb3a5a3d..ac7592db7d3fc 100644 --- a/stdlib/LinearAlgebra/test/blas.jl +++ b/stdlib/LinearAlgebra/test/blas.jl @@ -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_julia≈y_result_blas_lower + @test y_result_julia_lower ≈ y_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_julia≈y_result_blas_upper + @test y_result_julia_upper ≈ y_result_blas_upper end end