From c24945019c0e5188dcff8043d57bcb0913725be7 Mon Sep 17 00:00:00 2001
From: Ben Arthur <arthurb@hhmi.org>
Date: Tue, 18 Jan 2022 17:18:08 -0500
Subject: [PATCH] add option for read-only off-diagonal elements

---
 src/PackedArrays.jl | 28 +++++++++++++++++-----------
 test/runtests.jl    | 11 +++++++++--
 2 files changed, 26 insertions(+), 13 deletions(-)

diff --git a/src/PackedArrays.jl b/src/PackedArrays.jl
index cbf7f30..1ad0024 100644
--- a/src/PackedArrays.jl
+++ b/src/PackedArrays.jl
@@ -7,15 +7,15 @@ import LinearAlgebra: mul!, BLAS, BlasFloat, generic_matvecmul!, MulAddMul
 
 export SymmetricPacked
 
-struct SymmetricPacked{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
+struct SymmetricPacked{T,S<:AbstractMatrix{T},V} <: AbstractMatrix{T}
     tri::Vector{T}
     n::Int
     uplo::Char
 
-    function SymmetricPacked{T,S}(tri, n, uplo) where {T,S<:AbstractMatrix{T}}
+    function SymmetricPacked{T,S,V}(tri, n, uplo) where {T,S<:AbstractMatrix{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
 
@@ -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
@@ -63,20 +64,20 @@ 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,AbstractMatrix{T}}(pack(A, uplo), n, char_uplo(uplo))
+    SymmetricPacked{T,AbstractMatrix{T},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
 
 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)
 
@@ -97,6 +98,11 @@ end
     return r
 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
+
 function setindex!(A::SymmetricPacked, v, i::Int, j::Int)
     @boundscheck checkbounds(A, i, j)
     if A.uplo=='U'
diff --git a/test/runtests.jl b/test/runtests.jl
index 63af985..5fe6f18 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,7 @@ using PackedArrays, Test, LinearAlgebra
 A = collect(reshape(1:9.0,3,3))
 
 @testset "upper triangle" begin
-    APU = SymmetricPacked(A, :U)
+    APU = SymmetricPacked(A, :U, Val(:RW))
 
     @test APU[1,1] == A[1,1]
     @test APU[1,2] == A[1,2]
@@ -20,7 +20,7 @@ A = collect(reshape(1:9.0,3,3))
 end
 
 @testset "lower triangle" begin
-    APL = SymmetricPacked(A, :L)
+    APL = SymmetricPacked(A, :L, Val(:RW))
 
     @test APL[1,1] == A[1,1]
     @test APL[1,2] == A[2,1]
@@ -36,6 +36,13 @@ end
     @test APL[2,1] == 0
 end
 
+@testset "read-only" begin
+    APL = SymmetricPacked(A, :L)
+    @test_throws ArgumentError APL[1,3]=3
+    APL = SymmetricPacked(A, :L, Val(:RO))
+    @test_throws ArgumentError APL[1,3]=3
+end
+
 @testset "mul!" begin
     for uplo in [:U, :L]
         y = Float64[1,2,3]