Improve test coverage for RFP code and fix the broken parts

Author: andreasnoack

src/rectfullpacked.jl

@@ -2,8 +2,9 @@
 
 import Base: \
 import Base.LinAlg: BlasFloat
+import ..LinearAlgebra.LAPACK2: trttf!
 
-immutable HermitianRFP{T<:BlasFloat} <: AbstractMatrix{T}
+struct HermitianRFP{T<:BlasFloat} <: AbstractMatrix{T}
     data::Vector{T}
     transr::Char
     uplo::Char
@@ -83,14 +84,37 @@ end
 
 Base.copy(A::HermitianRFP) = HermitianRFP(copy(A.data), A.transr, A.uplo)
 
-immutable TriangularRFP{T<:BlasFloat} <: AbstractMatrix{T}
+struct TriangularRFP{T<:BlasFloat} <: AbstractMatrix{T}
     data::Vector{T}
     transr::Char
     uplo::Char
 end
-TriangularRFP(A::Matrix) = TriangularRFP(trttf!('N', 'U', A), 'N', 'U')
 
-function full(A::TriangularRFP)
+function TriangularRFP(A::StridedMatrix, uplo::Symbol = :U)
+    if uplo == :U
+        return TriangularRFP(trttf!('N', 'U', A), 'N', 'U')
+    elseif uplo == :L
+        return TriangularRFP(trttf!('N', 'L', A), 'N', 'L')
+    else
+        throw(ArgumentError("uplo must be either :U or :L but was :$uplo"))
+    end
+end
+
+function Base.size(A::TriangularRFP, i::Integer)
+    if i == 1 || i == 2
+        return (isqrt(8*length(A.data) + 1) - 1) >> 1
+    elseif i > 2
+        return 1
+    else
+        return size(A.data, i)
+    end
+end
+function Base.size(A::TriangularRFP)
+    n = size(A, 1)
+    return (n, n)
+end
+
+function Base.full(A::TriangularRFP)
     C = LAPACK2.tfttr!(A.transr, A.uplo, A.data)
     if A.uplo == 'U'
         return triu!(C)
@@ -99,7 +123,7 @@ function full(A::TriangularRFP)
     end
 end
 
-type CholeskyRFP{T<:BlasFloat} <: Factorization{T}
+struct CholeskyRFP{T<:BlasFloat} <: Factorization{T}
     data::Vector{T}
     transr::Char
     uplo::Char
@@ -109,10 +133,12 @@ Base.LinAlg.cholfact!{T<:BlasFloat}(A::HermitianRFP{T}) = CholeskyRFP(LAPACK2.pf
 Base.LinAlg.cholfact{T<:BlasFloat}(A::HermitianRFP{T}) = cholfact!(copy(A))
 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
 
 inv!(A::CholeskyRFP) = HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo)
-Base.LinAlg.inv(A::CholeskyRFP)  = inv(copy(A))
+Base.LinAlg.inv(A::CholeskyRFP)  = inv!(copy(A))
 Base.LinAlg.inv(A::HermitianRFP) = inv!(cholfact(A))

test/rectfullpacked.jl

@@ -1,18 +1,70 @@
 using Base.Test
-using LinearAlgebra
+import LinearAlgebra
+import LinearAlgebra: Ac_mul_A_RFP, TriangularRFP
 
 @testset "Rectuangular Full Pack Format" begin
 
-    @testset "Element type: $elty. Problem size: $n" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}),
-        n in (6, 7)
+    @testset "Core generic functionality" for n in (6, 7),
+                                           uplo in (:U, :L)
+
+        A = rand(10, n)
+
+        @testset "Hermitian" begin
+            AcA_RFP = Ac_mul_A_RFP(A, uplo)
+
+            @test size(AcA_RFP, 1) == n
+            @test size(AcA_RFP, 2) == n
+            @test size(AcA_RFP, 3) == 1
+            @test_throws BoundsError AcA_RFP[0, 1]
+            @test_throws BoundsError AcA_RFP[1, 0]
+            @test_throws BoundsError AcA_RFP[n + 1, 1]
+            @test_throws BoundsError AcA_RFP[1, n + 1]
+
+            @test AcA_RFP[2, 1]             == AcA_RFP[1, 2]
+            @test AcA_RFP[end - 2, end - 1] == AcA_RFP[end - 1, end - 2]
+        end
+
+        @testset "Triangular" begin
+            Atr_RFP = TriangularRFP(triu(A'A), uplo)
+            @test size(Atr_RFP, 1) == n
+            @test size(Atr_RFP, 2) == n
+            @test size(Atr_RFP, 3) == 1
+
+            # Indexing not implemented yet for Atr_RFP
+            # @test_throws BoundsError Atr_RFP[0, 1]
+            # @test_throws BoundsError Atr_RFP[1, 0]
+            # @test_throws BoundsError Atr_RFP[n + 1, 1]
+            # @test_throws BoundsError Atr_RFP[1, n + 1]
+
+            @test_broken Atr_RFP[2, 1]             == Atr_RFP[1, 2]
+            @test_broken Atr_RFP[end - 2, end - 1] == Atr_RFP[end - 1, end - 2]
+        end
+    end
+
+    @testset "Hermitian with element type: $elty. Problem size: $n" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}),
+                                                                           n in (6, 7),
+                                                                        uplo in (:L, :U)
 
         A       = rand(elty, 10, n)
         AcA     = A'A
-        AcA_RFP = LinearAlgebra.Ac_mul_A_RFP(A)
+        AcA_RFP = Ac_mul_A_RFP(A, uplo)
         o       = ones(elty, n)
 
         @test AcA   ≈ A'A
         @test AcA\o ≈ AcA_RFP\o
         @test inv(AcA) ≈ inv(AcA_RFP)
+        @test inv(cholfact(AcA)) ≈ inv(factorize(AcA_RFP))
+    end
+
+    @testset "Hermitian with element type: $elty. Problem size: $n" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}),
+                                                                           n in (6, 7),
+                                                                        uplo in (:L, :U)
+
+        A     = triu(rand(elty, n, n))
+        A_RFP = TriangularRFP(A)
+        o     = ones(elty, n)
+
+        @test_broken A ≈ A_RFP
+        @test        A ≈ full(A_RFP)
     end
 end