Out-of-place triu/tril for Symmetric in each branch

src/adjtrans.jl

@@ -332,6 +332,10 @@ wrapperop(::Transpose) = transpose
 _wrapperop(x) = wrapperop(x)
 _wrapperop(::Adjoint{<:Real}) = transpose
 
+# equivalent to wrapperop, but returns the type of the wrapper
+wrappertype(::Adjoint) = Adjoint
+wrappertype(::Transpose) = Transpose
+
 # the following fallbacks can be removed if Adjoint/Transpose are restricted to AbstractVecOrMat
 size(A::AdjOrTrans) = reverse(size(A.parent))
 axes(A::AdjOrTrans) = reverse(axes(A.parent))
@@ -391,7 +395,7 @@ similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, axes(A))
 similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims)
 
 # AbstractMatrix{T} constructor for adjtrans vector: preserve wrapped type
-AbstractMatrix{T}(A::AdjOrTransAbsVec) where {T} = wrapperop(A)(AbstractVector{T}(A.parent))
+AbstractMatrix{T}(A::AdjOrTransAbsVec) where {T} = wrappertype(A)(AbstractVector{T}(A.parent))
 
 # sundry basic definitions
 parent(A::AdjOrTrans) = A.parent

src/symmetric.jl

@@ -474,49 +474,49 @@ Base.conj!(A::HermOrSym) = typeof(A)(parentof_applytri(conj!, A), A.uplo)
 # tril/triu
 function tril(A::Hermitian, k::Integer=0)
     if A.uplo == 'U' && k <= 0
-        return tril!(copy(A.data'),k)
+        return tril_maybe_inplace(copy(A.data'),k)
     elseif A.uplo == 'U' && k > 0
-        return tril!(copy(A.data'),-1) + tril!(triu(A.data),k)
+        return tril_maybe_inplace(copy(A.data'),-1) + tril_maybe_inplace(triu(A.data),k)
     elseif A.uplo == 'L' && k <= 0
         return tril(A.data,k)
     else
-        return tril(A.data,-1) + tril!(triu!(copy(A.data')),k)
+        return tril(A.data,-1) + tril_maybe_inplace(triu_maybe_inplace(copy(A.data')),k)
     end
 end
 
 function tril(A::Symmetric, k::Integer=0)
     if A.uplo == 'U' && k <= 0
-        return tril!(copy(transpose(A.data)),k)
+        return tril_maybe_inplace(copy(transpose(A.data)),k)
     elseif A.uplo == 'U' && k > 0
-        return tril!(copy(transpose(A.data)),-1) + tril!(triu(A.data),k)
+        return tril_maybe_inplace(copy(transpose(A.data)),-1) + tril_maybe_inplace(triu(A.data),k)
     elseif A.uplo == 'L' && k <= 0
         return tril(A.data,k)
     else
-        return tril(A.data,-1) + tril!(triu!(copy(transpose(A.data))),k)
+        return tril(A.data,-1) + tril_maybe_inplace(triu_maybe_inplace(copy(transpose(A.data))),k)
     end
 end
 
 function triu(A::Hermitian, k::Integer=0)
     if A.uplo == 'U' && k >= 0
         return triu(A.data,k)
     elseif A.uplo == 'U' && k < 0
-        return triu(A.data,1) + triu!(tril!(copy(A.data')),k)
+        return triu(A.data,1) + triu_maybe_inplace(tril_maybe_inplace(copy(A.data')),k)
     elseif A.uplo == 'L' && k >= 0
-        return triu!(copy(A.data'),k)
+        return triu_maybe_inplace(copy(A.data'),k)
     else
-        return triu!(copy(A.data'),1) + triu!(tril(A.data),k)
+        return triu_maybe_inplace(copy(A.data'),1) + triu_maybe_inplace(tril(A.data),k)
     end
 end
 
 function triu(A::Symmetric, k::Integer=0)
     if A.uplo == 'U' && k >= 0
         return triu(A.data,k)
     elseif A.uplo == 'U' && k < 0
-        return triu(A.data,1) + triu!(tril!(copy(transpose(A.data))),k)
+        return triu(A.data,1) + triu_maybe_inplace(tril_maybe_inplace(copy(transpose(A.data))),k)
     elseif A.uplo == 'L' && k >= 0
-        return triu!(copy(transpose(A.data)),k)
+        return triu_maybe_inplace(copy(transpose(A.data)),k)
     else
-        return triu!(copy(transpose(A.data)),1) + triu!(tril(A.data),k)
+        return triu_maybe_inplace(copy(transpose(A.data)),1) + triu_maybe_inplace(tril(A.data),k)
     end
 end
 

src/triangular.jl

@@ -484,6 +484,11 @@ function tril!(A::UnitLowerTriangular, k::Integer=0)
     return tril!(LowerTriangular(A.data), k)
 end
 
+tril_maybe_inplace(A, k::Integer=0) = tril(A, k)
+triu_maybe_inplace(A, k::Integer=0) = triu(A, k)
+tril_maybe_inplace(A::StridedMatrix, k::Integer=0) = tril!(A, k)
+triu_maybe_inplace(A::StridedMatrix, k::Integer=0) = triu!(A, k)
+
 adjoint(A::LowerTriangular) = UpperTriangular(adjoint(A.data))
 adjoint(A::UpperTriangular) = LowerTriangular(adjoint(A.data))
 adjoint(A::UnitLowerTriangular) = UnitUpperTriangular(adjoint(A.data))

test/symmetric.jl

@@ -1206,4 +1206,18 @@ end
     end
 end
 
+@testset "triu/tril with immutable arrays" begin
+    A = ImmutableArray([1 2; 3 4])
+    for T in (Symmetric, Hermitian), uplo in (:U, :L)
+        H = T(A, uplo)
+        MH = Matrix(H)
+        @test triu(H,-1) == triu(MH,-1)
+        @test triu(H) == triu(MH)
+        @test triu(H,1) == triu(MH,1)
+        @test tril(H,1) == tril(MH,1)
+        @test tril(H) == tril(MH)
+        @test tril(H,-1) == tril(MH,-1)
+    end
+end
+
 end # module TestSymmetric

test/testhelpers/ImmutableArrays.jl

@@ -28,4 +28,7 @@ AbstractArray{T,N}(A::ImmutableArray{S,N}) where {S,T,N} = ImmutableArray(Abstra
 Base.copy(A::ImmutableArray) = ImmutableArray(copy(A.data))
 Base.zero(A::ImmutableArray) = ImmutableArray(zero(A.data))
 
+Base.adjoint(A::ImmutableArray) = ImmutableArray(adjoint(A.data))
+Base.transpose(A::ImmutableArray) = ImmutableArray(transpose(A.data))
+
 end