Overload copymutable_oftype

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "1.0.13"
+version = "1.1.0"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/ArrayLayouts.jl

@@ -270,6 +270,12 @@ copyto!(dest::AbstractMatrix, src::SubArray{<:Any,2,<:AdjOrTrans{<:Any,<:LayoutA
 # ambiguity from sparsematrix.jl
 copyto!(dest::LayoutMatrix, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
 copyto!(dest::SubArray{<:Any,2,<:LayoutMatrix}, src::SparseArrays.AbstractSparseMatrixCSC) = _copyto!(dest, src)
+if isdefined(LinearAlgebra, :copymutable_oftype)
+    LinearAlgebra.copymutable_oftype(A::Union{LayoutArray,Symmetric{<:Any,<:LayoutMatrix},Hermitian{<:Any,<:LayoutMatrix},
+                                                                UpperOrLowerTriangular{<:Any,<:LayoutMatrix},
+                                                                AdjOrTrans{<:Any,<:LayoutMatrix}}, ::Type{T}) where T = copymutable_oftype_layout(MemoryLayout(A), A, T)
+end
+copymutable_oftype_layout(_, A, ::Type{S}) where S = copyto!(similar(A, S), A)
 
 # avoid bad copy in Base
 Base.map(::typeof(copy), D::Diagonal{<:LayoutArray}) = Diagonal(map(copy, D.diag))

src/ldiv.jl

@@ -23,6 +23,8 @@ for Typ in (:Ldiv, :Rdiv)
     end
 end
 
+similar(A::Rdiv{<:DualLayout}, ::Type{T}, (ax1,ax2)) where T = dualadjoint(A.A)(similar(Array{T}, (ax2,)))
+
 @inline _ldivaxes(::Tuple{}, ::Tuple{}) = ()
 @inline _ldivaxes(::Tuple{}, Bax::Tuple) = Bax
 @inline _ldivaxes(::Tuple{<:Any}, ::Tuple{<:Any}) = ()
@@ -80,8 +82,18 @@ __ldiv!(_, F, B) = LinearAlgebra.ldiv!(F, B)
 
 @inline _ldiv!(dest, A, B; kwds...) = ldiv!(dest, factorize(A), B; kwds...)
 @inline _ldiv!(dest, A::Factorization, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
-@inline _ldiv!(dest, A::Transpose{<:Any,<:Factorization}, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
-@inline _ldiv!(dest, A::Adjoint{<:Any,<:Factorization}, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
+
+if VERSION ≥ v"1.10-"
+    using LinearAlgebra: TransposeFactorization, AdjointFactorization
+else
+    const TransposeFactorization = Transpose
+    const AdjointFactorization = Adjoint
+
+end
+@inline _ldiv!(dest, A::TransposeFactorization{<:Any,<:Factorization}, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
+@inline _ldiv!(dest, A::AdjointFactorization{<:Any,<:Factorization}, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
+
+
 
 @inline ldiv(A, B; kwds...) = materialize(Ldiv(A,B); kwds...)
 @inline rdiv(A, B; kwds...) = materialize(Rdiv(A,B); kwds...)
@@ -94,7 +106,10 @@ __ldiv!(_, F, B) = LinearAlgebra.ldiv!(F, B)
 
 @inline materialize!(M::Ldiv) = _ldiv!(M.A, M.B)
 @inline materialize!(M::Rdiv) = ldiv!(M.B', M.A')'
-@inline copyto!(dest::AbstractArray, M::Rdiv; kwds...) = copyto!(dest', Ldiv(M.B', M.A'); kwds...)'
+@inline function copyto!(dest::AbstractArray, M::Rdiv; kwds...)
+    adj = dualadjoint(dest)
+    adj(copyto!(adj(dest), Ldiv(adj(M.B), adj(M.A)); kwds...))
+end
 @inline copyto!(dest::AbstractArray, M::Ldiv; kwds...) = _ldiv!(dest, M.A, copy(M.B); kwds...)
 
 const MatLdivVec{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractVector{V}}
@@ -139,6 +154,8 @@ macro _layoutldiv(Typ)
 
         LinearAlgebra.ldiv!(A::Bidiagonal, B::$Typ; kwds...) = ArrayLayouts.ldiv!(A,B; kwds...)
 
+        LinearAlgebra.rdiv!(A::AbstractMatrix, B::$Typ; kwds...) = ArrayLayouts.rdiv!(A,B; kwds...)
+
 
         (\)(A::$Typ, x::AbstractVector; kwds...) = ArrayLayouts.ldiv(A,x; kwds...)
         (\)(A::$Typ, x::AbstractMatrix; kwds...) = ArrayLayouts.ldiv(A,x; kwds...)
@@ -171,6 +188,8 @@ macro _layoutldiv(Typ)
         (/)(A::$Typ, D::Diagonal; kwds...) = ArrayLayouts.rdiv(A,D; kwds...)
 
         (/)(x::$Typ, A::$Typ; kwds...) = ArrayLayouts.rdiv(x,A; kwds...)
+        (/)(D::Adjoint{<:Any,<:AbstractVector}, A::$Typ; kwds...) = ArrayLayouts.rdiv(D,A; kwds...)
+        (/)(D::Transpose{<:Any,<:AbstractVector}, A::$Typ; kwds...) = ArrayLayouts.rdiv(D,A; kwds...)
     end
     if Typ ≠ :LayoutVector
         ret = quote

test/test_layoutarray.jl

@@ -64,8 +64,14 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test A[kr,jr] == A.A[kr,jr]
         end
         b = randn(5)
+        B = randn(5,5)
         for Tri in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
             @test ldiv!(Tri(A), copy(b)) ≈ ldiv!(Tri(A.A), copy(b)) ≈ Tri(A.A) \ MyVector(b)
+            @test ldiv!(Tri(A), copy(B)) ≈ ldiv!(Tri(A.A), copy(B)) ≈ Tri(A.A) \ MyMatrix(B)
+            if VERSION ≥ v"1.9"
+                @test rdiv!(copy(b)', Tri(A)) ≈ rdiv!(copy(b)', Tri(A.A)) ≈ MyVector(b)' / Tri(A.A)
+                @test rdiv!(copy(B), Tri(A)) ≈ rdiv!(copy(B), Tri(A.A)) ≈ B / Tri(A.A)
+            end
             @test lmul!(Tri(A), copy(b)) ≈ lmul!(Tri(A.A), copy(b)) ≈ Tri(A.A) * MyVector(b)
         end
 
@@ -112,7 +118,7 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test cholesky(S, CRowMaximum()) \ b ≈ ldiv!(cholesky(Matrix(S), CRowMaximum()), copy(MyVector(b)))
             @test cholesky(S) \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S)) \ b
             @test cholesky(S) \ b ≈ Symmetric(Matrix(S)) \ MyVector(b)
-            if VERSION >= v"1.9-"
+            if VERSION >= v"1.9"
                 @test S \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S)) \ b
                 @test S \ b ≈ Symmetric(Matrix(S)) \ MyVector(b)
             end
@@ -122,19 +128,19 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test cholesky(S,CRowMaximum()).U ≈ cholesky(Matrix(S),CRowMaximum()).U
             @test cholesky(S) \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S), :L) \ b
             @test cholesky(S) \ b ≈ Symmetric(Matrix(S), :L) \ MyVector(b)
-            if VERSION >= v"1.9-"
+            if VERSION >= v"1.9"
                 @test S \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S), :L) \ b
                 @test S \ b ≈ Symmetric(Matrix(S), :L) \ MyVector(b)
             end
 
             @testset "ldiv!" begin
                 c = MyVector(randn(5))
-                if VERSION < v"1.9-"
+                if VERSION < v"1.9"
                     @test_broken ldiv!(lu(A), MyVector(copy(c))) ≈ A \ c
                 else
                     @test ldiv!(lu(A), MyVector(copy(c))) ≈ A \ c
                 end
-                if VERSION < v"1.9-" || VERSION >= v"1.10-"
+                if VERSION < v"1.9" || VERSION >= v"1.10-"
                     @test_throws ErrorException ldiv!(qr(A), MyVector(copy(c)))
                 else
                     @test_throws MethodError ldiv!(qr(A), MyVector(copy(c)))
@@ -221,18 +227,24 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test_broken ldiv!(A, t) ≈ A\t
             @test ldiv!(A, copy(X)) ≈ A\X
             @test A\T ≈ A\T̃
-            VERSION >= v"1.9-" && @test A/T ≈ A/T̃
+            VERSION >= v"1.9" && @test A/T ≈ A/T̃
             @test_broken ldiv!(A, T) ≈ A\T
             @test B\A ≈ B\Matrix(A)
             @test D \ A ≈ D \ Matrix(A)
             @test transpose(B)\A ≈ transpose(B)\Matrix(A) ≈ Transpose(B)\A ≈ Adjoint(B)\A
             @test B'\A ≈ B'\Matrix(A)
             @test A\A ≈ I
-            VERSION >= v"1.9-" && @test A/A ≈ I
+            VERSION >= v"1.9" && @test A/A ≈ I
             @test A\MyVector(x) ≈ A\x
             @test A\MyMatrix(X) ≈ A\X
 
-            VERSION >= v"1.9-" && @test A/A ≈ A.A / A.A
+            if VERSION >= v"1.9"
+                @test A/A ≈ A.A / A.A
+                @test x' / A ≈ x' / A.A
+                @test transpose(x) / A ≈ transpose(x) / A.A 
+                @test transpose(x) / A isa Transpose
+                @test x' / A isa Adjoint
+            end
         end
 
         @testset "dot" begin
@@ -455,6 +467,13 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
         @test UpperTriangular(A) * UnitUpperTriangular(A') ≈ UpperTriangular(A.A) * UnitUpperTriangular(A.A')
         @test UpperTriangular(A') * UnitUpperTriangular(A') ≈ UpperTriangular(A.A') * UnitUpperTriangular(A.A')
     end
+
+    if isdefined(LinearAlgebra, :copymutable_oftype)
+        @testset "copymutable_oftype" begin
+            A = MyMatrix(randn(3,3))
+            @test LinearAlgebra.copymutable_oftype(A, BigFloat) == A
+        end
+    end
 end
 
 struct MyUpperTriangular{T} <: AbstractMatrix{T}
@@ -498,5 +517,5 @@ triangulardata(A::MyUpperTriangular) = triangulardata(A.A)
     @test_skip lmul!(U,view(copy(B),collect(1:5),1:5)) ≈ U * B
 
     @test MyMatrix(A) / U ≈ A / U
-    VERSION >= v"1.9-" && @test U / MyMatrix(A) ≈ U / A
+    VERSION >= v"1.9" && @test U / MyMatrix(A) ≈ U / A
 end