diff --git a/src/GenericLinearAlgebra.jl b/src/GenericLinearAlgebra.jl
index 2ab778b..c39aa32 100644
--- a/src/GenericLinearAlgebra.jl
+++ b/src/GenericLinearAlgebra.jl
@@ -11,5 +11,4 @@ include("eigenSelfAdjoint.jl")
 include("eigenGeneral.jl")
 include("tridiag.jl")
 include("svd.jl")
-include("rectfullpacked.jl")
 end
diff --git a/src/rectfullpacked.jl b/src/rectfullpacked.jl
deleted file mode 100644
index e6f17ae..0000000
--- a/src/rectfullpacked.jl
+++ /dev/null
@@ -1,174 +0,0 @@
-# Rectangular Full Packed Matrices
-
-using LinearAlgebra: BlasFloat
-
-import ..LAPACK2: trttf!
-import Base: \
-import LinearAlgebra: ldiv!
-
-struct HermitianRFP{T<:BlasFloat} <: AbstractMatrix{T}
-    data::Vector{T}
-    transr::Char
-    uplo::Char
-end
-
-function Base.size(A::HermitianRFP, 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::HermitianRFP)
-    n = size(A, 1)
-    return (n, n)
-end
-
-function Base.getindex(A::HermitianRFP, i::Integer, j::Integer)
-    n = size(A, 1)
-    n2 = n >> 1
-    if i < 1 || i > n || j < 0 || j > n
-        throw(BoundsError(A, (i, j)))
-    end
-    if A.uplo == 'L'
-        if i < j
-            return conj(A[j, i])
-        end
-
-        if j <= n2 + isodd(n)
-            if A.transr == 'N'
-                return A.data[i+iseven(n)+(j-1)*(n+iseven(n))]
-            else
-                return conj(A.data[(i-1)*(n+iseven(n))+j+iseven(n)])
-            end
-        else
-            if A.transr == 'N'
-                return conj(A.data[(i-n2-1)*(n+iseven(n))+j-n2-isodd(n)])
-            else
-                return A.data[i-n2-isodd(n)+(j-n2-1)*(n+iseven(n))]
-            end
-        end
-    else
-        if i > j
-            return conj(A[j, i])
-        end
-
-        if j > n2
-            if A.transr == 'N'
-                return A.data[i+(j-n2-1)*(n+iseven(n))]
-            else
-                return conj(A.data[(i-n2-1)*(n+iseven(n))+j])
-            end
-        else
-            if A.transr == 'N'
-                return conj(A.data[(i-1)*(n+iseven(n))+j+n2+1])
-            else
-                return A.data[i-n2-isodd(n)+(j-n2-1)*(n+iseven(n))]
-            end
-        end
-    end
-end
-
-function Ac_mul_A_RFP(A::Matrix{T}, uplo = :U) where {T<:BlasFloat}
-    n = size(A, 2)
-    if uplo == :U
-        C = LAPACK2.sfrk!(
-            'N',
-            'U',
-            T <: Complex ? 'C' : 'T',
-            1.0,
-            A,
-            0.0,
-            Vector{T}(undef, (n * (n + 1)) >> 1),
-        )
-        return HermitianRFP(C, 'N', 'U')
-    elseif uplo == :L
-        C = LAPACK2.sfrk!(
-            'N',
-            'L',
-            T <: Complex ? 'C' : 'T',
-            1.0,
-            A,
-            0.0,
-            Vector{T}(undef, (n * (n + 1)) >> 1),
-        )
-        return HermitianRFP(C, 'N', 'L')
-    else
-        throw(ArgumentError("uplo must be either :L or :U"))
-    end
-end
-
-Base.copy(A::HermitianRFP) = HermitianRFP(copy(A.data), A.transr, A.uplo)
-
-struct TriangularRFP{T<:BlasFloat} <: AbstractMatrix{T}
-    data::Vector{T}
-    transr::Char
-    uplo::Char
-end
-
-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
-
-Base.copy(A::TriangularRFP) = TriangularRFP(copy(A.data), A.transr, A.uplo)
-
-function Base.Array(A::TriangularRFP)
-    C = LAPACK2.tfttr!(A.transr, A.uplo, A.data)
-    if A.uplo == 'U'
-        return triu!(C)
-    else
-        return tril!(C)
-    end
-end
-
-LinearAlgebra.inv!(A::TriangularRFP) =
-    TriangularRFP(LAPACK2.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} =
-    LAPACK2.tfsm!(A.transr, 'L', A.uplo, 'N', 'N', one(T), A.data, B)
-(\)(A::TriangularRFP, B::StridedVecOrMat) = ldiv!(A, copy(B))
-
-struct CholeskyRFP{T<:BlasFloat} <: Factorization{T}
-    data::Vector{T}
-    transr::Char
-    uplo::Char
-end
-
-LinearAlgebra.cholesky!(A::HermitianRFP{T}) where {T<:BlasFloat} =
-    CholeskyRFP(LAPACK2.pftrf!(A.transr, A.uplo, copy(A.data)), A.transr, A.uplo)
-LinearAlgebra.cholesky(A::HermitianRFP{T}) where {T<:BlasFloat} = cholesky!(copy(A))
-LinearAlgebra.factorize(A::HermitianRFP) = cholesky(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) = cholesky(A) \ B
-
-LinearAlgebra.inv!(A::CholeskyRFP) =
-    HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo)
-LinearAlgebra.inv(A::CholeskyRFP) = LinearAlgebra.inv!(copy(A))
-LinearAlgebra.inv(A::HermitianRFP) = LinearAlgebra.inv!(cholesky(A))
diff --git a/test/rectfullpacked.jl b/test/rectfullpacked.jl
deleted file mode 100644
index 45a0336..0000000
--- a/test/rectfullpacked.jl
+++ /dev/null
@@ -1,82 +0,0 @@
-using Test
-import GenericLinearAlgebra
-import GenericLinearAlgebra: Ac_mul_A_RFP, TriangularRFP
-
-@testset "Rectangular Full Pack Format" begin
-
-    @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 = 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(cholesky(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 = lu(rand(elty, n, n)).U
-        A = uplo == :U ? A : copy(A')
-        A_RFP = TriangularRFP(A, uplo)
-        o = ones(elty, n)
-
-        @test_broken A ≈ A_RFP
-        @test A ≈ Array(A_RFP)
-        @test A \ o ≈ A_RFP \ o
-        @test inv(A) ≈ Array(inv(A_RFP))
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index ecc35c7..4de5615 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,6 +8,5 @@ include("eigenselfadjoint.jl")
 include("eigengeneral.jl")
 include("tridiag.jl")
 include("svd.jl")
-include("rectfullpacked.jl")
 include("lapack.jl")
 # end