Implement wrapped-SMatrix + UniformScaling (#1332)

## Context

`UniformScaling` operations on structured wrappers backed by
`StaticMatrix` were falling back to generic `LinearAlgebra` methods.
That path uses mutable copies, which caused type leakage from
`SMatrix`-backed wrappers to `MMatrix`-backed wrappers. In some cases
(for example `I - Symmetric(SMatrix(...))`), it could also fail due to
mutation of immutable storage.

I don't think there is an easy way to solve this problem in general due
to how these operations are defined in base, so instead I've just
written some appropriate wrappers.

## What this changes

This PR adds `UniformScaling` specialisations for structured wrappers
over `StaticMatrix` in `src/linalg.jl`:
  - `Symmetric`
  - `Hermitian`
  - `UpperTriangular`
  - `LowerTriangular`
  - `UnitUpperTriangular` (result type `UpperTriangular`)
  - `UnitLowerTriangular` (result type `LowerTriangular`)

I only implement methods for `A + J` and `J - A` for brevity, with
`uniformscaling.jl:L156` handling the reverse argument orders.

I include a special case for the some of Hermitian + complex
UniformScaling which removes the wrapper to be type stable, consistent
with base LinearAlgebra.jl.

## Tests

I've extended `test/linalg.jl` (`Interaction with UniformScaling`) to
cover all of the cases mentioned.

Author: THargreaves

Project.toml

@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.16"
+version = "1.9.17"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/linalg.jl

@@ -40,6 +40,50 @@ import Base: +, -, *, /, \
     end
 end
 
+@inline _plus_uniform_parent(A, λ) = _plus_uniform(Size(A), parent(A), λ)
+@inline _minus_uniform_parent(A, λ) = _plus_uniform(Size(A), -parent(A), λ)
+@inline _sym_uplo(A) = A.uplo == 'U' ? :U : :L
+
+# Self-adjoint wrappers over StaticMatrix with UniformScaling
+@inline +(A::Symmetric{T,<:StaticMatrix{N,M,T}}, J::UniformScaling) where {N,M,T} =
+    Symmetric(_plus_uniform_parent(A, J.λ), _sym_uplo(A))
+@inline -(J::UniformScaling, A::Symmetric{T,<:StaticMatrix{N,M,T}}) where {N,M,T} =
+    Symmetric(_minus_uniform_parent(A, J.λ), _sym_uplo(A))
+
+@inline +(A::Hermitian{T,<:StaticMatrix{N,M,T}}, J::UniformScaling) where {N,M,T} =
+    Hermitian(_plus_uniform_parent(A, J.λ), _sym_uplo(A))
+@inline -(J::UniformScaling, A::Hermitian{T,<:StaticMatrix{N,M,T}}) where {N,M,T} =
+    Hermitian(_minus_uniform_parent(A, J.λ), _sym_uplo(A))
+
+# Lose Hermitian wrapper when adding complex UniformScaling, in-line with Base behavior
+@inline function +(A::Hermitian{T,<:StaticMatrix{N,M,T}}, J::UniformScaling{<:Complex}) where {N,M,T}
+    TS = Base.promote_op(+, eltype(A), typeof(J))
+    _plus_uniform(Size(A), similar_type(A, TS)(A), J.λ)
+end
+@inline function -(J::UniformScaling{<:Complex}, A::Hermitian{T,<:StaticMatrix{N,M,T}}) where {N,M,T}
+    TS = Base.promote_op(+, eltype(A), typeof(J))
+    _plus_uniform(Size(A), -similar_type(A, TS)(A), J.λ)
+end
+
+# Triangular wrappers over StaticMatrix with UniformScaling
+for TWR in (:UpperTriangular, :LowerTriangular)
+    @eval begin
+        @inline +(A::$TWR{T,<:StaticMatrix{N,M,T}}, J::UniformScaling) where {N,M,T} = $TWR(_plus_uniform_parent(A, J.λ))
+        @inline -(J::UniformScaling, A::$TWR{T,<:StaticMatrix{N,M,T}}) where {N,M,T} = $TWR(_minus_uniform_parent(A, J.λ))
+    end
+end
+
+# Unit triangular wrappers over StaticMatrix with UniformScaling
+@inline _plus_uniform_materialized(A, λ) = _plus_uniform(Size(A), similar_type(A, promote_type(eltype(A), typeof(λ)))(A), λ)
+@inline _minus_uniform_materialized(A, λ) = _plus_uniform(Size(A), -similar_type(A, promote_type(eltype(A), typeof(λ)))(A), λ)
+
+for (TWR1, TWR2) in ((:UnitUpperTriangular, :UpperTriangular), (:UnitLowerTriangular, :LowerTriangular))
+    @eval begin
+        @inline +(A::$TWR1{T,<:StaticMatrix{N,M,T}}, J::UniformScaling) where {N,M,T} = $TWR2(_plus_uniform_materialized(A, J.λ))
+        @inline -(J::UniformScaling, A::$TWR1{T,<:StaticMatrix{N,M,T}}) where {N,M,T} = $TWR2(_minus_uniform_materialized(A, J.λ))
+    end
+end
+
 @inline *(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ * b
 @inline *(a::StaticMatrix, b::UniformScaling) = a * b.λ
 @inline \(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ \ b

test/linalg.jl

@@ -107,6 +107,55 @@ end
         @test @inferred(I * @SMatrix([0 1; 2 3])) === @SMatrix [0 1; 2 3]
         @test @inferred(@SMatrix([0 1; 2 3]) / I) === @SMatrix [0.0 1.0; 2.0 3.0]
         @test @inferred(I \ @SMatrix([0 1; 2 3])) === @SMatrix [0.0 1.0; 2.0 3.0]
+
+        A = @SMatrix [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
+        A_dynamic = Matrix(A)
+
+        for WR in (Symmetric, Hermitian)
+            B = WR(A, :L)
+            B_dynamic = WR(A_dynamic, :L)
+
+            @test @inferred(B + I) isa WR{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(I + B) isa WR{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(B - I) isa WR{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(I - B) isa WR{Float64,SMatrix{3,3,Float64,9}}
+            @test (B + I).uplo == 'L'
+            @test (I - B).uplo == 'L'
+
+            @test B + I == B_dynamic + I
+            @test I + B == I + B_dynamic
+            @test B - I == B_dynamic - I
+            @test I - B == I - B_dynamic
+        end
+
+        B = Hermitian(complex.(A), :U)
+        B_dynamic = Hermitian(complex.(A_dynamic), :U)
+
+        # Hermitian wrapper should be dropped iff the UniformScaling is complex, in-line with Base behavior
+        @test @inferred(B + I) isa Hermitian{ComplexF64,SMatrix{3,3,ComplexF64,9}}
+        @test @inferred(I - B) isa Hermitian{ComplexF64,SMatrix{3,3,ComplexF64,9}}
+        @test @inferred(B + (1 + 2im)I) isa SMatrix{3,3,ComplexF64,9}
+        @test @inferred((1 + 2im)I - B) isa SMatrix{3,3,ComplexF64,9}
+        @test B + (1 + 2im)I == B_dynamic + (1 + 2im)I
+        @test (1 + 2im)I - B == (1 + 2im)I - B_dynamic
+
+        for (WR1, WR2) in ((UpperTriangular, UpperTriangular),
+            (LowerTriangular, LowerTriangular),
+            (UnitUpperTriangular, UpperTriangular),
+            (UnitLowerTriangular, LowerTriangular))
+            B = WR1(A)
+            B_dynamic = WR1(A_dynamic)
+
+            @test @inferred(B + I) isa WR2{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(I + B) isa WR2{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(B - I) isa WR2{Float64,SMatrix{3,3,Float64,9}}
+            @test @inferred(I - B) isa WR2{Float64,SMatrix{3,3,Float64,9}}
+
+            @test B + I == B_dynamic + I
+            @test I + B == I + B_dynamic
+            @test B - I == B_dynamic - I
+            @test I - B == I - B_dynamic
+        end
     end
 
     @testset "Constructors from UniformScaling" begin