Add rdiv!, ldiv! methods and tests (#19)

* Non-working version with rdiv!, ldiv!

* Change tests

* Change from character to symbol in comparison.

* Test transposed rdiv! and ldiv!, pass 'C' for cong. tr.

Author: dmbates

src/lapack.jl

@@ -206,7 +206,7 @@ for (f, elty) in (
             chkstride1(B)
             m, n = size(B, 1), size(B, 2)
             k, l = transr == 'N' ? size(A) : reverse(size(A))
-            if k - (2l ≤ k) != m
+            if k - (2l ≤ k) != (side == 'L' ? m : n)
                 throw(
                     DimensionMismatch(
                         "First dimension of B must equal $(k - (2l ≤ k)), got $m",

src/triangular.jl

@@ -56,6 +56,24 @@ LinearAlgebra.inv!(A::TriangularRFP) =
     TriangularRFP(LAPACK_RFP.tftri!(A.transr, A.uplo, 'N', A.data), A.transr, A.uplo)
 LinearAlgebra.inv(A::TriangularRFP) = LinearAlgebra.inv!(copy(A))
 
-ldiv!(A::TriangularRFP{T}, B::StridedVecOrMat{T}) where {T} =
+LinearAlgebra.ldiv!(A::TriangularRFP{T}, B::StridedVecOrMat{T}) where {T} =
     LAPACK_RFP.tfsm!(A.transr, 'L', A.uplo, 'N', 'N', one(T), A.data, B)
+function LinearAlgebra.ldiv!(
+    A::Adjoint{T, TriangularRFP{T}},
+    B::StridedVecOrMat{T}) where {T}
+    Ap = A.parent
+    tr = T <: Complex ? 'C' : 'T'
+    return LAPACK_RFP.tfsm!(Ap.transr, 'L', Ap.uplo, tr, 'N', one(T), Ap.data, B)
+end
+
+LinearAlgebra.rdiv!(A::StridedVecOrMat{T}, B::TriangularRFP{T}) where {T} =
+    LAPACK_RFP.tfsm!(B.transr, 'R', B.uplo, 'N', 'N', one(T), B.data, A)
+function LinearAlgebra.rdiv!(
+    A::StridedVecOrMat{T},
+    B::Adjoint{T, TriangularRFP{T}}) where {T}
+    Bp = B.parent
+    tr = T <: Complex ? 'C' : 'T'
+    return LAPACK_RFP.tfsm!(Bp.transr, 'R', Bp.uplo, tr, 'N', one(T), Bp.data, A)
+end
+
 (\)(A::TriangularRFP, B::StridedVecOrMat) = ldiv!(A, copy(B))

test/runtests.jl

@@ -84,19 +84,24 @@ import RectangularFullPacked: Ac_mul_A_RFP, TriangularRFP
             Complex{Float32},
             Complex{Float64},
         ),
-        n in (6, 7),
+        n in (6,7),
         uplo in (:L, :U),
         transr in (:N, elty <: Complex ? :C : :T)
 
         A = lu(rand(elty, n, n)).U
         A = uplo == :U ? A : copy(A')
         A_RFP = TriangularRFP(A, uplo; transr)
+        Atri = uplo == :U ? UpperTriangular(copy(A)) : LowerTriangular(A)
         o = ones(elty, n)
 
         @test A ≈ A_RFP
         @test A ≈ Array(A_RFP)
         @test A \ o ≈ A_RFP \ o
-        @test inv(A) ≈ Array(inv(A_RFP))
+        @test Array(inv(A)) ≈ Array(inv(A_RFP))
+        @test ldiv!(Atri, copy(o)) ≈ ldiv!(A_RFP, copy(o))
+        @test ldiv!(Atri', copy(o)) ≈ ldiv!(A_RFP', copy(o))
+        @test rdiv!(collect(o'), Atri) ≈ rdiv!(collect(o'), A_RFP)
+        @test rdiv!(collect(o'), Atri') ≈ rdiv!(collect(o'), A_RFP')
     end
 
     @testset "In-place scalar multiplication" begin