Use LAPACK for T\B where T is triangular (#96)

* Use LAPACK for T\B where T is triangular

* v0.8.7

* Update triangular.jl

* Update test_ldiv.jl

Author: dlfivefifty

Project.toml

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

src/triangular.jl

@@ -215,25 +215,40 @@ for UNIT in ('U', 'N')
         @eval @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,<:AbstractColumnMajor},
                                             <:AbstractStridedLayout}) =
             BLAS.trsv!($UPLO, 'N', $UNIT, triangulardata(M.A), M.B)
+        @eval @inline materialize!(M::BlasMatLdivMat{<:TriangularLayout{$UPLO,$UNIT,<:AbstractColumnMajor},
+                                            <:AbstractStridedLayout}) =
+            LAPACK.trtrs!($UPLO, 'N', $UNIT, triangulardata(M.A), M.B)
     end
 
     @eval begin
         @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractRowMajor},
                                                         <:AbstractStridedLayout}) =
             BLAS.trsv!('L', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)
+        @inline materialize!(M::BlasMatLdivMat{<:TriangularLayout{'U',$UNIT,<:AbstractRowMajor},
+                                                        <:AbstractStridedLayout}) =
+            LAPACK.trtrs!('L', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)
 
         @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractRowMajor},
                                                         <:AbstractStridedLayout}) =
             BLAS.trsv!('U', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)
+        @inline materialize!(M::BlasMatLdivMat{<:TriangularLayout{'L',$UNIT,<:AbstractRowMajor},
+                                                        <:AbstractStridedLayout}) =
+            LAPACK.trtrs!('U', 'T', $UNIT, transpose(triangulardata(M.A)), M.B)
 
 
         @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                                         <:AbstractStridedLayout}) =
             BLAS.trsv!('L', 'C', $UNIT, triangulardata(M.A)', M.B)
+        @inline materialize!(M::BlasMatLdivMat{<:TriangularLayout{'U',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
+                                                        <:AbstractStridedLayout}) =
+            LAPACK.trtrs!('L', 'C', $UNIT, triangulardata(M.A)', M.B)            
 
         @inline materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                                         <:AbstractStridedLayout}) =
             BLAS.trsv!('U', 'C', $UNIT, triangulardata(M.A)', M.B)
+        @inline materialize!(M::BlasMatLdivMat{<:TriangularLayout{'L',$UNIT,<:ConjLayout{<:AbstractRowMajor}},
+                                                        <:AbstractStridedLayout}) =
+            LAPACK.trtrs!('U', 'C', $UNIT, triangulardata(M.A)', M.B)            
     end
 end
 

test/test_ldiv.jl

@@ -86,11 +86,24 @@ import ArrayLayouts: ApplyBroadcastStyle, QRCompactWYQLayout, QRCompactWYLayout,
 
     @testset "Triangular \\ matrix" begin
         A = randn(5,5)
-        b = randn(5,5)
-        M =  Ldiv(UpperTriangular(A), b)
+        B = randn(5,5)
+        M =  Ldiv(UpperTriangular(A), B)
         @test Base.Broadcast.broadcastable(M) === M
-        @test UpperTriangular(A) \ b ≈ copyto!(similar(b) , Ldiv(UpperTriangular(A), b)) ≈
-            (b .= Ldiv(UpperTriangular(A), b))
+        @test UpperTriangular(A) \ B ≈ copyto!(similar(B) , Ldiv(UpperTriangular(A), B)) ≈ (B .= Ldiv(UpperTriangular(A), B))
+        @test UpperTriangular(A') \ B ≈ copyto!(similar(B) , Ldiv(UpperTriangular(A'), B)) ≈ (B .= Ldiv(UpperTriangular(A'), B))
+        @test LowerTriangular(A') \ B ≈ copyto!(similar(B) , Ldiv(LowerTriangular(A'), B)) ≈ (B .= Ldiv(LowerTriangular(A'), B))
+        @test UnitUpperTriangular(A) \ B ≈ copyto!(similar(B) , Ldiv(UnitUpperTriangular(A), B)) ≈ (B .= Ldiv(UnitUpperTriangular(A), B))
+        @test UnitUpperTriangular(A') \ B ≈ copyto!(similar(B) , Ldiv(UnitUpperTriangular(A'), B)) ≈ (B .= Ldiv(UnitUpperTriangular(A'), B))
+        @test UnitLowerTriangular(A') \ B ≈ copyto!(similar(B) , Ldiv(UnitLowerTriangular(A'), B)) ≈ (B .= Ldiv(UnitLowerTriangular(A'), B))        
+
+        C = randn(5,5) + im * randn(5,5)
+        D = randn(5,5) + im * randn(5,5)
+        @test UpperTriangular(C) \ D ≈ copyto!(similar(D) , Ldiv(UpperTriangular(C), D)) ≈ (D .= Ldiv(UpperTriangular(C), D))
+        @test UpperTriangular(C') \ D ≈ copyto!(similar(D) , Ldiv(UpperTriangular(C'), D)) ≈ (D .= Ldiv(UpperTriangular(C'), D))
+        @test LowerTriangular(C') \ D ≈ copyto!(similar(D) , Ldiv(LowerTriangular(C'), D)) ≈ (D .= Ldiv(LowerTriangular(C'), D))
+        @test UnitUpperTriangular(C) \ D ≈ copyto!(similar(D) , Ldiv(UnitUpperTriangular(C), D)) ≈ (D .= Ldiv(UnitUpperTriangular(C), D))
+        @test UnitUpperTriangular(C') \ D ≈ copyto!(similar(D) , Ldiv(UnitUpperTriangular(C'), D)) ≈ (D .= Ldiv(UnitUpperTriangular(C'), D))
+        @test UnitLowerTriangular(C') \ D ≈ copyto!(similar(D) , Ldiv(UnitLowerTriangular(C'), D)) ≈ (D .= Ldiv(UnitLowerTriangular(C'), D))
     end
 
     @testset "Int" begin