Merge pull request #24163 from Sacha0/similarsym

make similar(A<:HermOrSym, opts...) consistently preserve A's storage type and uplo flag

base/linalg/symmetric.jl

@@ -137,16 +137,19 @@ function setindex!(A::Hermitian, v, i::Integer, j::Integer)
     end
 end
 
-similar(A::Symmetric, ::Type{T}) where {T} = Symmetric(similar(A.data, T))
-# Hermitian version can be simplified when check for imaginary part of
-# diagonal in Hermitian has been removed
+# For A<:Union{Symmetric,Hermitian}, similar(A[, neweltype]) should yield a matrix with the same
+# symmetry type, uplo flag, and underlying storage type as A. The following methods cover these cases.
+similar(A::Symmetric, ::Type{T}) where {T} = Symmetric(similar(parent(A), T), ifelse(A.uplo == 'U', :U, :L))
+# If the the Hermitian constructor's check ascertaining that the wrapped matrix's
+# diagonal is strictly real is removed, the following method can be simplified.
 function similar(A::Hermitian, ::Type{T}) where T
-    B = similar(A.data, T)
-    for i = 1:size(A,1)
-        B[i,i] = 0
-    end
-    return Hermitian(B)
+    B = similar(parent(A), T)
+    for i in 1:size(B, 1) B[i, i] = 0 end
+    return Hermitian(B, ifelse(A.uplo == 'U', :U, :L))
 end
+# On the other hand, similar(A, [neweltype,] shape...) should yield a matrix of the underlying
+# storage type of A (not wrapped in a symmetry type). The following method covers these cases.
+similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = similar(parent(A), T, dims)
 
 # Conversion
 convert(::Type{Matrix}, A::Symmetric) = copytri!(convert(Matrix, copy(A.data)), A.uplo)

test/linalg/symmetric.jl

@@ -444,3 +444,17 @@ end
         @test norm(inv(Hermitian(H))*(H*ones(8)) .- 1) ≈ 0 atol = 1e-5
     end
 end
+
+@testset "similar should preserve underlying storage type and uplo flag" begin
+    m, n = 4, 3
+    sparsemat = sprand(m, m, 0.5)
+    for SymType in (Symmetric, Hermitian)
+        symsparsemat = SymType(sparsemat)
+        @test isa(similar(symsparsemat), typeof(symsparsemat))
+        @test similar(symsparsemat).uplo == symsparsemat.uplo
+        @test isa(similar(symsparsemat, Float32), SymType{Float32,<:SparseMatrixCSC{Float32}})
+        @test similar(symsparsemat, Float32).uplo == symsparsemat.uplo
+        @test isa(similar(symsparsemat, (n, n)), typeof(sparsemat))
+        @test isa(similar(symsparsemat, Float32, (n, n)), SparseMatrixCSC{Float32})
+    end
+end