Make inv and solves work for RFP

Author: andreasnoack

src/lapack.jl

@@ -2,7 +2,7 @@ module LAPACK2
 
     using Base.LinAlg: BlasInt, chkstride1, LAPACKException
     using Base.LinAlg.BLAS: @blasfunc
-    using Base.LinAlg.LAPACK: chkdiag, chkuplo
+    using Base.LinAlg.LAPACK: chkdiag, chkside, chkuplo
 
     # LAPACK wrappers
     import Base.BLAS.@blasfunc
@@ -604,12 +604,19 @@ for (f, elty) in ((:dtfsm_, :Float64),
                    (:ctfsm_, :Complex64))
 
     @eval begin
-        function pftrs!(transr::Char, side::Char, uplo::Char, trans::Char, diag::Char, alpha::Real, A::StridedVector{$elty}, B::StridedMatrix{$elty})
+        function tfsm!(transr::Char,
+                       side::Char,
+                       uplo::Char,
+                       trans::Char,
+                       diag::Char,
+                       alpha::$elty,
+                       A::StridedVector{$elty},
+                       B::StridedVecOrMat{$elty})
             chkuplo(uplo)
             chkside(side)
             chkdiag(diag)
             chkstride1(B)
-            m, n = size(B)
+            m, n = size(B, 1), size(B, 2)
             if round(Int, div(sqrt(8length(A)), 2)) != m
                 throw(DimensionMismatch("First dimension of B must equal $(round(Int, div(sqrt(8length(A)), 2))), got $m"))
             end
@@ -623,7 +630,7 @@ for (f, elty) in ((:dtfsm_, :Float64),
                 &diag, &m, &n, &alpha,
                 A, B, &ldb)
 
-            B
+            return B
         end
     end
 end

src/rectfullpacked.jl

@@ -128,6 +128,10 @@ end
 Base.LinAlg.inv!(A::TriangularRFP) = TriangularRFP(LAPACK2.tftri!(A.transr, A.uplo, 'N', A.data), A.transr, A.uplo)
 Base.LinAlg.inv(A::TriangularRFP)  = Base.LinAlg.inv!(copy(A))
 
+A_ldiv_B!(A::TriangularRFP{T}, B::StridedVecOrMat{T}) where T =
+    LAPACK2.tfsm!(A.transr, 'L', A.uplo, 'N', 'N', one(T), A.data, B)
+(\)(A::TriangularRFP, B::StridedVecOrMat) = A_ldiv_B!(A, copy(B))
+
 struct CholeskyRFP{T<:BlasFloat} <: Factorization{T}
     data::Vector{T}
     transr::Char
@@ -141,8 +145,8 @@ Base.LinAlg.factorize(A::HermitianRFP) = cholfact(A)
 Base.copy(F::CholeskyRFP{T}) where T = CholeskyRFP{T}(copy(F.data), F.transr, F.uplo)
 
 # Solve
-\(A::CholeskyRFP, B::StridedVecOrMat) = LAPACK2.pftrs!(A.transr, A.uplo, A.data, copy(B))
-\(A::HermitianRFP, B::StridedVecOrMat) = cholfact(A)\B
+(\)(A::CholeskyRFP, B::StridedVecOrMat) = LAPACK2.pftrs!(A.transr, A.uplo, A.data, copy(B))
+(\)(A::HermitianRFP, B::StridedVecOrMat) = cholfact(A)\B
 
 Base.LinAlg.inv!(A::CholeskyRFP) = HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo)
 Base.LinAlg.inv(A::CholeskyRFP)  = Base.LinAlg.inv!(copy(A))

test/rectfullpacked.jl

@@ -60,12 +60,14 @@ import LinearAlgebra: Ac_mul_A_RFP, TriangularRFP
                                                                            n in (6, 7),
                                                                         uplo in (:L, :U)
 
-        A     = triu(rand(elty, n, n))
-        A_RFP = TriangularRFP(A)
+        A     = lufact(rand(elty, n, n))[:U]
+        A     = uplo == :U ? A : A'
+        A_RFP = TriangularRFP(A, uplo)
         o     = ones(elty, n)
 
         @test_broken A ≈ A_RFP
         @test        A ≈ full(A_RFP)
+        @test      A\o ≈ A_RFP\o
         @test   inv(A) ≈ full(inv(A_RFP))
     end
 end