Turn broken rdiv tests to proper tests, optimize ldiv by diagonal (#113)

Author: dkarrasch

Project.toml

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

README.md

@@ -1,11 +1,15 @@
 # ArrayLayouts.jl
+
 A Julia package for describing array layouts and more general fast linear algebra
 
 [![Build Status](https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/workflows/CI/badge.svg)](https://github.com/JuliaMatrices/ArrayLayouts.jl/actions)
 [![codecov](https://codecov.io/gh/JuliaLinearAlgebra/ArrayLayouts.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaMatrices/ArrayLayouts.jl)
 
-This package implements a trait-based framework for describing array layouts such as column major, row major, etc. that can be dispatched 
-to appropriate BLAS or optimised Julia linear algebra routines. This supports a much wider class of matrix types than Julia's in-built `StridedArray`. Here is an example:
+This package implements a trait-based framework for describing array layouts such as column
+major, row major, etc. that can be dispatched to appropriate BLAS or optimised Julia linear
+algebra routines. This supports a much wider class of matrix types than Julia's in-built
+`StridedArray`. Here is an example:
+
 ```julia
 julia> using ArrayLayouts
 
@@ -22,7 +26,9 @@ julia> @time muladd!(1.0, V, x, 0.0, y); # ArrayLayouts does and is 3x faster as
 
 ## Internal design
 
-The package is based on assigning a `MemoryLayout` to every array, which is used for dispatch. For example, 
+The package is based on assigning a `MemoryLayout` to every array, which is used for
+dispatch. For example,
+
 ```julia
 julia> MemoryLayout(A) # Each column of A is column major, and columns stored in order
 DenseColumnMajor()
@@ -34,9 +40,11 @@ julia> MemoryLayout(V) # A symmetric version, whose storage is DenseColumnMajor
 SymmetricLayout{DenseColumnMajor}()
 ```
 
-This is then used by `muladd!(α, A, B, β, C)`, `ArrayLayouts.lmul!(A, B)`, and `ArrayLayouts.rmul!(A, B)` to
-lower to the correct BLAS calls via lazy objects `MulAdd(α, A, B, β, C)`, `Lmul(A, B)`, `Rmul(A, B)` which are materialized,
-in analogy to `Base.Broadcasted`. 
+This is then used by `muladd!(α, A, B, β, C)`, `ArrayLayouts.lmul!(A, B)`, and
+`ArrayLayouts.rmul!(A, B)` to lower to the correct BLAS calls via lazy objects
+`MulAdd(α, A, B, β, C)`, `Lmul(A, B)`, `Rmul(A, B)` which are materialized, in analogy to
+`Base.Broadcasted`.
 
-Note there is also a higher level function `mul(A, B)` that materializes via `Mul(A, B)`, which uses the layout of `A` and `B`
-to further reduce to either `MulAdd`, `Lmul`, and `Rmul`.
+Note there is also a higher level function `mul(A, B)` that materializes via `Mul(A, B)`,
+which uses the layout of `A` and `B` to further reduce to either `MulAdd`, `Lmul`, and
+`Rmul`.

src/diagonal.jl

@@ -35,8 +35,8 @@ function materialize!(M::Ldiv{<:DiagonalLayout})
     M.B
 end
 
-copy(M::Ldiv{<:DiagonalLayout,<:DiagonalLayout}) = diagonal(inv.(M.A.diag) .* M.B.diag)
-copy(M::Ldiv{<:DiagonalLayout}) = inv.(M.A.diag) .* M.B
+copy(M::Ldiv{<:DiagonalLayout,<:DiagonalLayout}) = diagonal(M.A.diag .\ M.B.diag)
+copy(M::Ldiv{<:DiagonalLayout}) = M.A.diag .\ M.B
 copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout}}) = inv(getindex_value(M.A.diag)) .* M.B
 copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = diagonal(inv(getindex_value(M.A.diag)) .* M.B.diag)
 

src/ldiv.jl

@@ -143,6 +143,10 @@ macro _layoutldiv(Typ)
         (\)(A::$Typ, x::AbstractMatrix) = ArrayLayouts.ldiv(A,x)
 
         (\)(x::AbstractMatrix, A::$Typ) = ArrayLayouts.ldiv(x,A)
+        (\)(x::LinearAlgebra.HermOrSym, A::$Typ) = ArrayLayouts.ldiv(x,A)
+        if VERSION < v"1.9-" # disambiguation
+            \(x::LinearAlgebra.HermOrSym{<:Any,<:StridedMatrix}, A::$Typ) = ArrayLayouts.ldiv(x,A)
+        end
         (\)(x::UpperTriangular, A::$Typ) = ArrayLayouts.ldiv(x,A)
         (\)(x::UnitUpperTriangular, A::$Typ) = ArrayLayouts.ldiv(x,A)
         (\)(x::LowerTriangular, A::$Typ) = ArrayLayouts.ldiv(x,A)

test/test_layoutarray.jl

@@ -94,6 +94,7 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test lu!(copy(A)).factors ≈ lu(A.A).factors
             b = randn(5)
             @test all(A \ b .≡ A.A \ b .≡ A.A \ MyVector(b) .≡ ldiv!(lu(A.A), copy(MyVector(b))))
+            @test all(A \ b .≡ ldiv!(lu(A), copy(MyVector(b))) .≡ ldiv!(lu(A), copy(b)))
             @test all(lu(A).L .≡ lu(A.A).L)
             @test all(lu(A).U .≡ lu(A.A).U)
             @test lu(A).p == lu(A.A).p
@@ -105,13 +106,26 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I))
             @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
             @test cholesky(S, CRowMaximum()).U ≈ cholesky(Matrix(S), CRowMaximum()).U
+            @test cholesky!(deepcopy(S), CRowMaximum()).U ≈ cholesky(Matrix(S), CRowMaximum()).U
             @test cholesky(S) \ b ≈ cholesky(Matrix(S)) \ b ≈ cholesky(Matrix(S)) \ MyVector(b)
             @test cholesky(S, CRowMaximum()) \ b ≈ cholesky(Matrix(S), CRowMaximum()) \ b
             @test cholesky(S, CRowMaximum()) \ b ≈ ldiv!(cholesky(Matrix(S), CRowMaximum()), copy(MyVector(b)))
-
-            S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I),:L)
+            @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-"
+                @test S \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S)) \ b
+                @test S \ b ≈ Symmetric(Matrix(S)) \ MyVector(b)
+            end
+
+            S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I), :L)
             @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
             @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-"
+                @test S \ b ≈ Matrix(S) \ b ≈ Symmetric(Matrix(S), :L) \ b
+                @test S \ b ≈ Symmetric(Matrix(S), :L) \ MyVector(b)
+            end
         end
         Bin = randn(5,5)
         B = MyMatrix(copy(Bin))
@@ -173,18 +187,18 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test_broken ldiv!(A, t) ≈ A\t
             @test ldiv!(A, copy(X)) ≈ A\X
             @test A\T ≈ A\T̃
-            @test_broken 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
-            @test_broken A/A ≈ I
+            VERSION >= v"1.9-" && @test A/A ≈ I
             @test A\MyVector(x) ≈ A\x
             @test A\MyMatrix(X) ≈ A\X
 
-            @test_broken A/A ≈ A.A / A.A
+            VERSION >= v"1.9-" && @test A/A ≈ A.A / A.A
         end
 
         @testset "dot" begin
@@ -392,5 +406,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
-    @test_broken U / MyMatrix(A) ≈ U / A
+    VERSION >= v"1.9-" && @test U / MyMatrix(A) ≈ U / A
 end