Skip to content

Commit

Permalink
add spr! to LinearAlgebra.BLAS (#42830)
Browse files Browse the repository at this point in the history
  • Loading branch information
bjarthur authored Nov 22, 2021
1 parent a21cc80 commit ae336ab
Show file tree
Hide file tree
Showing 3 changed files with 110 additions and 0 deletions.
1 change: 1 addition & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ Standard library changes
#### Package Manager

#### LinearAlgebra
* The BLAS submodule now supports the level-2 BLAS subroutine `spr!` ([#42830]).

#### Markdown

Expand Down
66 changes: 66 additions & 0 deletions stdlib/LinearAlgebra/src/blas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ export
sbmv!,
sbmv,
spmv!,
spr!,
symv!,
symv,
trsv!,
Expand Down Expand Up @@ -1145,6 +1146,71 @@ Return the updated `y`.
"""
spmv!

### spr!, (SP) symmetric packed matrix-vector operation defined as A := alpha*x*x' + A
for (fname, elty) in ((:dspr_, :Float64),
(:sspr_, :Float32))
@eval begin
function spr!(uplo::AbstractChar,
n::Integer,
α::$elty,
x::Union{Ptr{$elty}, AbstractArray{$elty}},
incx::Integer,
AP::Union{Ptr{$elty}, AbstractArray{$elty}})

ccall((@blasfunc($fname), libblastrampoline), Cvoid,
(Ref{UInt8}, # uplo,
Ref{BlasInt}, # n,
Ref{$elty}, # α,
Ptr{$elty}, # x,
Ref{BlasInt}, # incx,
Ptr{$elty}, # AP,
Clong), # length of uplo
uplo,
n,
α,
x,
incx,
AP,
1)
return AP
end
end
end

function spr!(uplo::AbstractChar,
α::Real, x::Union{DenseArray{T}, AbstractVector{T}},
AP::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasReal}
require_one_based_indexing(AP, x)
N = length(x)
if 2*length(AP) < N*(N + 1)
throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
end
return spr!(uplo, N, convert(T, α), x, stride(x, 1), AP)
end

"""
spr!(uplo, α, x, AP)
Update matrix `A` as `α*A*x*x'`, where `A` is a symmetric matrix provided
in packed format `AP` and `x` is a vector.
With `uplo = 'U'`, the array AP must contain the upper triangular part of the
symmetric 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
symmetric 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 input `α` must be real.
The array inputs `x` and `AP` must all be of `Float32` or `Float64` type.
Return the updated `AP`.
"""
spr!

### hbmv, (HB) Hermitian banded matrix-vector multiplication
for (fname, elty) in ((:zhbmv_,:ComplexF64),
(:chbmv_,:ComplexF32))
Expand Down
43 changes: 43 additions & 0 deletions stdlib/LinearAlgebra/test/blas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,49 @@ Random.seed!(100)
end
end

# spr!
if elty in (Float32, Float64)
@testset "spr! $elty" begin
α = rand(elty)
M = rand(elty, n, n)
AL = Symmetric(M, :L)
AU = Symmetric(M, :U)
x = rand(elty, n)

function pack(A, uplo)
AP = elty[]
for j in 1:n
for i in (uplo==:L ? (j:n) : (1:j))
push!(AP, A[i,j])
end
end
return AP
end

ALP_result_julia_lower = pack*x*x' + AL, :L)
ALP_result_blas_lower = pack(AL, :L)
BLAS.spr!('L', α, x, ALP_result_blas_lower)
@test ALP_result_julia_lower ALP_result_blas_lower
ALP_result_blas_lower = append!(pack(AL, :L), ones(elty, 10))
BLAS.spr!('L', α, x, ALP_result_blas_lower)
@test ALP_result_julia_lower ALP_result_blas_lower[1:end-10]
ALP_result_blas_lower = reshape(pack(AL, :L), 1, length(ALP_result_julia_lower), 1)
BLAS.spr!('L', α, x, ALP_result_blas_lower)
@test ALP_result_julia_lower vec(ALP_result_blas_lower)

AUP_result_julia_upper = pack*x*x' + AU, :U)
AUP_result_blas_upper = pack(AU, :U)
BLAS.spr!('U', α, x, AUP_result_blas_upper)
@test AUP_result_julia_upper AUP_result_blas_upper
AUP_result_blas_upper = append!(pack(AU, :U), ones(elty, 10))
BLAS.spr!('U', α, x, AUP_result_blas_upper)
@test AUP_result_julia_upper AUP_result_blas_upper[1:end-10]
AUP_result_blas_upper = reshape(pack(AU, :U), 1, length(AUP_result_julia_upper), 1)
BLAS.spr!('U', α, x, AUP_result_blas_upper)
@test AUP_result_julia_upper vec(AUP_result_blas_upper)
end
end

#trsm
A = triu(rand(elty,n,n))
B = rand(elty,(n,n))
Expand Down

0 comments on commit ae336ab

Please sign in to comment.