Skip to content

Commit

Permalink
add option for writable off-diagonal elements
Browse files Browse the repository at this point in the history
  • Loading branch information
@bjarthur
bjarthur committed Feb 6, 2022
1 parent 553c0f1 commit f3d0e50
Show file tree
Hide file tree
Showing 3 changed files with 100 additions and 15 deletions.
44 changes: 44 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,3 +101,47 @@ julia> mul!(copy(y), AP, x, α, β) # or use the equivalent julian interface
0.9592679333345454
0.27783932357851576
```

`SymmetricPacked` differs in one important aspect to `LinearAlgebra.Symmetric`
in that it optionally permits off diagonal elements to be set. Continuing
from the example above:

```
julia> AP[1,1]=5
5
julia> AP
5×5 SymmetricPacked{Float64, Matrix{Float64}, Val{:RO}()}:
5.0 0.736451 0.704095 0.842738 0.207618
0.736451 0.63549 0.861638 0.834968 0.0660772
0.704095 0.861638 0.131074 0.727387 0.0995039
0.842738 0.834968 0.727387 0.96513 0.616451
0.207618 0.0660772 0.0995039 0.616451 0.369353
julia> AP[1,5]=5
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
[1] setindex!(A::SymmetricPacked{Float64, Matrix{Float64}, Val{:RO}()}, v::Int64, i::Int64, j::Int64)
@ PackedArrays ~/projects/darshan/PackedArrays/src/PackedArrays.jl:114
[2] top-level scope
@ REPL[7]:1
julia> APRW = SymmetricPacked(A, :L, Val(:RW)) # default is :RO
5×5 SymmetricPacked{Float64, Matrix{Float64}, Val{:RW}()}:
0.364856 0.0797498 0.152769 0.201612 0.283603
0.0797498 0.63549 0.00271082 0.0429332 0.818754
0.152769 0.00271082 0.131074 0.682401 0.76548
0.201612 0.0429332 0.682401 0.96513 0.284176
0.283603 0.818754 0.76548 0.284176 0.369353
julia> APRW[1,5]=5
5
julia> APRW
5×5 SymmetricPacked{Float64, Matrix{Float64}, Val{:RW}()}:
0.364856 0.0797498 0.152769 0.201612 5.0
0.0797498 0.63549 0.00271082 0.0429332 0.818754
0.152769 0.00271082 0.131074 0.682401 0.76548
0.201612 0.0429332 0.682401 0.96513 0.284176
5.0 0.818754 0.76548 0.284176 0.369353
```
37 changes: 22 additions & 15 deletions src/SymmetricFormats.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,15 @@ import LinearAlgebra: mul!, BLAS, BlasFloat, generic_matvecmul!, MulAddMul

export SymmetricPacked, packedsize

struct SymmetricPacked{T,S<:AbstractVecOrMat{<:T}} <: AbstractMatrix{T}
struct SymmetricPacked{T,S<:AbstractVecOrMat{<:T},V} <: AbstractMatrix{T}
tri::Vector{T}
n::Int
uplo::Char

function SymmetricPacked{T,S}(tri, n, uplo) where {T,S<:AbstractVecOrMat{<:T}}
function SymmetricPacked{T,S,V}(tri, n, uplo) where {T,S<:AbstractVecOrMat{<:T},V}
require_one_based_indexing(tri)
uplo=='U' || uplo=='L' || throw(ArgumentError("uplo must be either 'U' (upper) or 'L' (lower)"))
new{T,S}(tri, n, uplo)
new{T,S,V}(tri, n, uplo)
end
end

Expand All @@ -33,10 +33,11 @@ function pack(A::AbstractMatrix{T}, uplo::Symbol) where {T}
end

"""
SymmetricPacked(A, uplo=:U)
SymmetricPacked(A, uplo=:U, offdiag=Val(:RO))
Construct a `Symmetric` matrix in packed form of the upper (if `uplo = :U`)
or lower (if `uplo = :L`) triangle of the matrix `A`.
or lower (if `uplo = :L`) triangle of the matrix `A`. `offdiag` specifies
whether elements not on the diagaonal can be set (if `:RW`) or not (if `:RO`).
# Examples
```jldoctest
Expand All @@ -63,13 +64,13 @@ julia> Base.summarysize(AP)
184
```
"""
function SymmetricPacked(A::AbstractMatrix{T}, uplo::Symbol=:U) where {T}
function SymmetricPacked(A::AbstractMatrix{T}, uplo::Symbol=:U, offdiag=Val(:RO)) where {T}
n = checksquare(A)
SymmetricPacked{T,typeof(A)}(pack(A, uplo), n, char_uplo(uplo))
SymmetricPacked{T,typeof(A),offdiag}(pack(A, uplo), n, char_uplo(uplo))
end

function SymmetricPacked(x::SymmetricPacked{T,S}) where{T,S}
SymmetricPacked{T,S}(T.(x.tri), x.n, x.uplo)
function SymmetricPacked(x::SymmetricPacked{T,S,V}) where{T,S,V}
SymmetricPacked{T,S,V}(T.(x.tri), x.n, x.uplo)
end

function SymmetricPacked(V::AbstractVector{T}, uplo::Symbol=:U) where {T}
Expand All @@ -80,9 +81,9 @@ end

checksquare(x::SymmetricPacked) = x.n

convert(::Type{SymmetricPacked{T,S}}, x::SymmetricPacked) where {T,S} = SymmetricPacked{T,S}(T.(x.tri), x.n, x.uplo)
convert(::Type{SymmetricPacked{T,S,V}}, x::SymmetricPacked) where {T,S,V} = SymmetricPacked{T,S}(T.(x.tri), x.n, x.uplo)

unsafe_convert(::Type{Ptr{T}}, A::SymmetricPacked{T,S}) where {T,S} = Base.unsafe_convert(Ptr{T}, A.tri)
unsafe_convert(::Type{Ptr{T}}, A::SymmetricPacked{T,S,V}) where {T,S,V} = Base.unsafe_convert(Ptr{T}, A.tri)

size(A::SymmetricPacked) = (A.n,A.n)

Expand Down Expand Up @@ -115,8 +116,7 @@ end
return r
end

function setindex!(A::SymmetricPacked, v, i::Int, j::Int)
i!=j && throw(ArgumentError("Cannot set a non-diagonal index in a symmetric matrix"))
function _setindex!(A::SymmetricPacked, v, i::Int, j::Int)
@boundscheck checkbounds(A, i, j)
if A.uplo=='U'
i,j = minmax(i,j)
Expand All @@ -128,15 +128,22 @@ function setindex!(A::SymmetricPacked, v, i::Int, j::Int)
return v
end

function setindex!(A::SymmetricPacked{T,S,Val(:RO)}, v, i::Int, j::Int) where {T,S}
i!=j && throw(ArgumentError("Cannot set a non-diagonal index in a symmetric matrix"))
_setindex!(A, v, i, j)
end

setindex!(A::SymmetricPacked{T,S,Val(:RW)}, v, i::Int, j::Int) where {T,S} = _setindex!(A, v, i, j)

function copy(A::SymmetricPacked{T,S}) where {T,S}
B = copy(A.tri)
SymmetricPacked{T,S}(B, A.n, A.uplo)
end

@inline function mul!(y::StridedVector{T},
AP::SymmetricPacked{T,<:StridedMatrix},
AP::SymmetricPacked{T,<:StridedMatrix,V},
x::StridedVector{T},
α::Number, β::Number) where {T<:BlasFloat}
α::Number, β::Number) where {T<:BlasFloat,V}
alpha, beta = promote(α, β, zero(T))
if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
BLAS.spmv!(AP.uplo, alpha, AP.tri, x, beta, y)
Expand Down
34 changes: 34 additions & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,40 @@ A = collect(reshape(1:9.0,3,3))
@test_throws ArgumentError APL[1,3]=3
end

@testset "read/write upper" begin
APU = SymmetricPacked(A, :U, Val(:RW))

@test APU[1,1] == A[1,1]
@test APU[1,2] == A[1,2]
@test APU[1,3] == A[1,3]
@test APU[2,1] == A[1,2]
@test APU[2,2] == A[2,2]
@test APU[2,3] == A[2,3]
@test APU[3,1] == A[1,3]
@test APU[3,2] == A[2,3]
@test APU[3,3] == A[3,3]

APU[1,2] = 0
@test APU[2,1] == 0
end

@testset "read/write lower" begin
APL = SymmetricPacked(A, :L, Val(:RW))

@test APL[1,1] == A[1,1]
@test APL[1,2] == A[2,1]
@test APL[1,3] == A[3,1]
@test APL[2,1] == A[2,1]
@test APL[2,2] == A[2,2]
@test APL[2,3] == A[3,2]
@test APL[3,1] == A[3,1]
@test APL[3,2] == A[3,2]
@test APL[3,3] == A[3,3]

APL[1,2] = 0
@test APL[2,1] == 0
end

@testset "mul!" begin
for uplo in [:U, :L]
y = Float64[1,2,3]
Expand Down

0 comments on commit f3d0e50

Please sign in to comment.