From cb11ac7e06b8a0e34e29984f70f459011b7f0272 Mon Sep 17 00:00:00 2001 From: Laszlo Papp Date: Sun, 19 Jan 2020 20:38:24 +0100 Subject: [PATCH] Minor updates to BLAS/hpmv! as in #34320. [#34439] Follow the syntax, docstring and testing conventions elaborated in #34320. --- stdlib/LinearAlgebra/src/blas.jl | 24 ++++++++++++++++-------- stdlib/LinearAlgebra/test/blas.jl | 23 +++++++++++++---------- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/stdlib/LinearAlgebra/src/blas.jl b/stdlib/LinearAlgebra/src/blas.jl index 08867fb24c875..21471fd3b6708 100644 --- a/stdlib/LinearAlgebra/src/blas.jl +++ b/stdlib/LinearAlgebra/src/blas.jl @@ -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}}, @@ -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 -. +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 a2ecf8416f555..8aa644314d38f 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