Generalize storage of factorizations. (JuliaLang#42594)

Introduce additional type parameters to remove hard-coded arrays from structures.

stdlib/LinearAlgebra/docs/src/index.md

@@ -60,7 +60,7 @@ julia> A = [1.5 2 -4; 3 -1 -6; -10 2.3 4]
  -10.0   2.3   4.0
 
 julia> factorize(A)
-LU{Float64, Matrix{Float64}}
+LU{Float64, Matrix{Float64}, Vector{Int64}}
 L factor:
 3×3 Matrix{Float64}:
   1.0    0.0       0.0
@@ -84,7 +84,7 @@ julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
  -4.0  -3.0   5.0
 
 julia> factorize(B)
-BunchKaufman{Float64, Matrix{Float64}}
+BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
 D factor:
 3×3 Tridiagonal{Float64, Vector{Float64}}:
  -1.64286   0.0   ⋅

stdlib/LinearAlgebra/src/bunchkaufman.jl

@@ -28,7 +28,7 @@ julia> A = [1 2; 2 3]
  2  3
 
 julia> S = bunchkaufman(A) # A gets wrapped internally by Symmetric(A)
-BunchKaufman{Float64, Matrix{Float64}}
+BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
 D factor:
 2×2 Tridiagonal{Float64, Vector{Float64}}:
  -0.333333  0.0
@@ -48,7 +48,7 @@ julia> d == S.D && u == S.U && p == S.p
 true
 
 julia> S = bunchkaufman(Symmetric(A, :L))
-BunchKaufman{Float64, Matrix{Float64}}
+BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
 D factor:
 2×2 Tridiagonal{Float64, Vector{Float64}}:
  3.0   0.0
@@ -63,22 +63,25 @@ permutation:
  1
 ```
 """
-struct BunchKaufman{T,S<:AbstractMatrix} <: Factorization{T}
+struct BunchKaufman{T,S<:AbstractMatrix,P<:AbstractVector{<:Integer}} <: Factorization{T}
     LD::S
-    ipiv::Vector{BlasInt}
+    ipiv::P
     uplo::Char
     symmetric::Bool
     rook::Bool
     info::BlasInt
 
-    function BunchKaufman{T,S}(LD, ipiv, uplo, symmetric, rook, info) where {T,S<:AbstractMatrix}
+    function BunchKaufman{T,S,P}(LD, ipiv, uplo, symmetric, rook, info) where {T,S<:AbstractMatrix,P<:AbstractVector}
         require_one_based_indexing(LD)
-        new(LD, ipiv, uplo, symmetric, rook, info)
+        new{T,S,P}(LD, ipiv, uplo, symmetric, rook, info)
     end
 end
-BunchKaufman(A::AbstractMatrix{T}, ipiv::Vector{BlasInt}, uplo::AbstractChar, symmetric::Bool,
-             rook::Bool, info::BlasInt) where {T} =
-        BunchKaufman{T,typeof(A)}(A, ipiv, uplo, symmetric, rook, info)
+BunchKaufman(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, uplo::AbstractChar,
+             symmetric::Bool, rook::Bool, info::BlasInt) where {T} =
+        BunchKaufman{T,typeof(A),typeof(ipiv)}(A, ipiv, uplo, symmetric, rook, info)
+# backwards-compatible constructors (remove with Julia 2.0)
+@deprecate(BunchKaufman(LD, ipiv, uplo, symmetric, rook, info) where {T,S},
+           BunchKaufman{T,S,typeof(ipiv)}(LD, ipiv, uplo, symmetric, rook, info))
 
 # iteration for destructuring into components
 Base.iterate(S::BunchKaufman) = (S.D, Val(:UL))
@@ -148,7 +151,7 @@ julia> A = [1 2; 2 3]
  2  3
 
 julia> S = bunchkaufman(A) # A gets wrapped internally by Symmetric(A)
-BunchKaufman{Float64, Matrix{Float64}}
+BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
 D factor:
 2×2 Tridiagonal{Float64, Vector{Float64}}:
  -0.333333  0.0
@@ -173,7 +176,7 @@ julia> S.U*S.D*S.U' - S.P*A*S.P'
  0.0  0.0
 
 julia> S = bunchkaufman(Symmetric(A, :L))
-BunchKaufman{Float64, Matrix{Float64}}
+BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
 D factor:
 2×2 Tridiagonal{Float64, Vector{Float64}}:
  3.0   0.0

stdlib/LinearAlgebra/src/cholesky.jl

@@ -127,7 +127,7 @@ julia> X = [1.0, 2.0, 3.0, 4.0];
 julia> A = X * X';
 
 julia> C = cholesky(A, RowMaximum(), check = false)
-CholeskyPivoted{Float64, Matrix{Float64}}
+CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}
 U factor with rank 1:
 4×4 UpperTriangular{Float64, Matrix{Float64}}:
  4.0  2.0  3.0  1.0
@@ -150,23 +150,25 @@ julia> l == C.L && u == C.U
 true
 ```
 """
-struct CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T}
+struct CholeskyPivoted{T,S<:AbstractMatrix,P<:AbstractVector{<:Integer}} <: Factorization{T}
     factors::S
     uplo::Char
-    piv::Vector{BlasInt}
+    piv::P
     rank::BlasInt
     tol::Real
     info::BlasInt
 
-    function CholeskyPivoted{T,S}(factors, uplo, piv, rank, tol, info) where {T,S<:AbstractMatrix}
+    function CholeskyPivoted{T,S,P}(factors, uplo, piv, rank, tol, info) where {T,S<:AbstractMatrix,P<:AbstractVector}
         require_one_based_indexing(factors)
-        new(factors, uplo, piv, rank, tol, info)
+        new{T,S,P}(factors, uplo, piv, rank, tol, info)
     end
 end
-function CholeskyPivoted(A::AbstractMatrix{T}, uplo::AbstractChar, piv::Vector{<:Integer},
-                            rank::Integer, tol::Real, info::Integer) where T
-    CholeskyPivoted{T,typeof(A)}(A, uplo, piv, rank, tol, info)
-end
+CholeskyPivoted(A::AbstractMatrix{T}, uplo::AbstractChar, piv::AbstractVector{<:Integer},
+                rank::Integer, tol::Real, info::Integer) where T =
+    CholeskyPivoted{T,typeof(A),typeof(piv)}(A, uplo, piv, rank, tol, info)
+# backwards-compatible constructors (remove with Julia 2.0)
+@deprecate(CholeskyPivoted{T,S}(factors, uplo, piv, rank, tol, info) where {T,S<:AbstractMatrix},
+           CholeskyPivoted{T,S,typeof(piv)}(factors, uplo, piv, rank, tol, info))
 
 
 # iteration for destructuring into components
@@ -306,7 +308,7 @@ end
 function cholesky!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix},
                    ::RowMaximum; tol = 0.0, check::Bool = true)
     AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
-    C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
+    C = CholeskyPivoted{eltype(AA),typeof(AA),typeof(piv)}(AA, A.uplo, piv, rank, tol, info)
     check && chkfullrank(C)
     return C
 end
@@ -438,7 +440,7 @@ julia> X = [1.0, 2.0, 3.0, 4.0];
 julia> A = X * X';
 
 julia> C = cholesky(A, RowMaximum(), check = false)
-CholeskyPivoted{Float64, Matrix{Float64}}
+CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}
 U factor with rank 1:
 4×4 UpperTriangular{Float64, Matrix{Float64}}:
  4.0  2.0  3.0  1.0

stdlib/LinearAlgebra/src/lq.jl

@@ -22,13 +22,13 @@ julia> A = [5. 7.; -2. -4.]
  -2.0  -4.0
 
 julia> S = lq(A)
-LQ{Float64, Matrix{Float64}}
+LQ{Float64, Matrix{Float64}, Vector{Float64}}
 L factor:
 2×2 Matrix{Float64}:
  -8.60233   0.0
   4.41741  -0.697486
 Q factor:
-2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}}:
+2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}:
  -0.581238  -0.813733
  -0.813733   0.581238
 
@@ -43,28 +43,31 @@ julia> l == S.L &&  q == S.Q
 true
 ```
 """
-struct LQ{T,S<:AbstractMatrix{T}} <: Factorization{T}
+struct LQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: Factorization{T}
     factors::S
-    τ::Vector{T}
+    τ::C
 
-    function LQ{T,S}(factors, τ) where {T,S<:AbstractMatrix{T}}
+    function LQ{T,S,C}(factors, τ) where {T,S<:AbstractMatrix{T},C<:AbstractVector{T}}
         require_one_based_indexing(factors)
-        new{T,S}(factors, τ)
+        new{T,S,C}(factors, τ)
     end
 end
-LQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = LQ{T,typeof(factors)}(factors, τ)
-function LQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T}
-    LQ(convert(AbstractMatrix{T}, factors), convert(Vector{T}, τ))
-end
+LQ(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} =
+    LQ{T,typeof(factors),typeof(τ)}(factors, τ)
+LQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} =
+    LQ(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ))
+# backwards-compatible constructors (remove with Julia 2.0)
+@deprecate(LQ{T,S}(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T,S},
+           LQ{T,S,typeof(τ)}(factors, τ))
 
 # iteration for destructuring into components
 Base.iterate(S::LQ) = (S.L, Val(:Q))
 Base.iterate(S::LQ, ::Val{:Q}) = (S.Q, Val(:done))
 Base.iterate(S::LQ, ::Val{:done}) = nothing
 
-struct LQPackedQ{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
+struct LQPackedQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: AbstractMatrix{T}
     factors::S
-    τ::Vector{T}
+    τ::C
 end
 
 
@@ -96,13 +99,13 @@ julia> A = [5. 7.; -2. -4.]
  -2.0  -4.0
 
 julia> S = lq(A)
-LQ{Float64, Matrix{Float64}}
+LQ{Float64, Matrix{Float64}, Vector{Float64}}
 L factor:
 2×2 Matrix{Float64}:
  -8.60233   0.0
   4.41741  -0.697486
 Q factor:
-2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}}:
+2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}:
  -0.581238  -0.813733
  -0.813733   0.581238
 
@@ -134,7 +137,7 @@ Array(A::LQ) = Matrix(A)
 
 adjoint(A::LQ) = Adjoint(A)
 Base.copy(F::Adjoint{T,<:LQ{T}}) where {T} =
-    QR{T,typeof(F.parent.factors)}(copy(adjoint(F.parent.factors)), copy(F.parent.τ))
+    QR{T,typeof(F.parent.factors),typeof(F.parent.τ)}(copy(adjoint(F.parent.factors)), copy(F.parent.τ))
 
 function getproperty(F::LQ, d::Symbol)
     m, n = size(F)

stdlib/LinearAlgebra/src/lu.jl

@@ -28,7 +28,7 @@ julia> A = [4 3; 6 3]
  6  3
 
 julia> F = lu(A)
-LU{Float64, Matrix{Float64}}
+LU{Float64, Matrix{Float64}, Vector{Int64}}
 L factor:
 2×2 Matrix{Float64}:
  1.0       0.0
@@ -47,24 +47,24 @@ julia> l == F.L && u == F.U && p == F.p
 true
 ```
 """
-struct LU{T,S<:AbstractMatrix{T}} <: Factorization{T}
+struct LU{T,S<:AbstractMatrix{T},P<:AbstractVector{<:Integer}} <: Factorization{T}
     factors::S
-    ipiv::Vector{BlasInt}
+    ipiv::P
     info::BlasInt
 
-    function LU{T,S}(factors, ipiv, info) where {T,S<:AbstractMatrix{T}}
+    function LU{T,S,P}(factors, ipiv, info) where {T, S<:AbstractMatrix{T}, P<:AbstractVector{<:Integer}}
         require_one_based_indexing(factors)
-        new{T,S}(factors, ipiv, info)
+        new{T,S,P}(factors, ipiv, info)
     end
 end
-function LU(factors::AbstractMatrix{T}, ipiv::Vector{BlasInt}, info::BlasInt) where {T}
-    LU{T,typeof(factors)}(factors, ipiv, info)
-end
-function LU{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T}
-    LU(convert(AbstractMatrix{T}, factors),
-       convert(Vector{BlasInt}, ipiv),
-       BlasInt(info))
-end
+LU(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, info::BlasInt) where {T} =
+    LU{T,typeof(factors),typeof(ipiv)}(factors, ipiv, info)
+LU{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T} =
+    LU(convert(AbstractMatrix{T}, factors), ipiv, BlasInt(info))
+# backwards-compatible constructors (remove with Julia 2.0)
+@deprecate(LU{T,S}(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
+                   info::BlasInt) where {T,S},
+           LU{T,S,typeof(ipiv)}(factors, ipiv, info))
 
 # iteration for destructuring into components
 Base.iterate(S::LU) = (S.L, Val(:U))
@@ -80,7 +80,7 @@ lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMaximum(); ch
 function lu!(A::StridedMatrix{T}, ::RowMaximum; check::Bool = true) where {T<:BlasFloat}
     lpt = LAPACK.getrf!(A)
     check && checknonsingular(lpt[3])
-    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
+    return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
 end
 function lu!(A::StridedMatrix{<:BlasFloat}, pivot::NoPivot; check::Bool = true)
     return generic_lufact!(A, pivot; check = check)
@@ -111,7 +111,7 @@ julia> A = [4. 3.; 6. 3.]
  6.0  3.0
 
 julia> F = lu!(A)
-LU{Float64, Matrix{Float64}}
+LU{Float64, Matrix{Float64}, Vector{Int64}}
 L factor:
 2×2 Matrix{Float64}:
  1.0       0.0
@@ -184,7 +184,7 @@ function generic_lufact!(A::StridedMatrix{T}, pivot::Union{RowMaximum,NoPivot} =
         end
     end
     check && checknonsingular(info, pivot)
-    return LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
+    return LU{T,typeof(A),typeof(ipiv)}(A, ipiv, convert(BlasInt, info))
 end
 
 function lutype(T::Type)
@@ -256,7 +256,7 @@ julia> A = [4 3; 6 3]
  6  3
 
 julia> F = lu(A)
-LU{Float64, Matrix{Float64}}
+LU{Float64, Matrix{Float64}, Vector{Int64}}
 L factor:
 2×2 Matrix{Float64}:
  1.0       0.0
@@ -295,13 +295,13 @@ end
 
 function LU{T}(F::LU) where T
     M = convert(AbstractMatrix{T}, F.factors)
-    LU{T,typeof(M)}(M, F.ipiv, F.info)
+    LU{T,typeof(M),typeof(F.ipiv)}(M, F.ipiv, F.info)
 end
-LU{T,S}(F::LU) where {T,S} = LU{T,S}(convert(S, F.factors), F.ipiv, F.info)
+LU{T,S,P}(F::LU) where {T,S,P} = LU{T,S,P}(convert(S, F.factors), convert(P, F.ipiv), F.info)
 Factorization{T}(F::LU{T}) where {T} = F
 Factorization{T}(F::LU) where {T} = LU{T}(F)
 
-copy(A::LU{T,S}) where {T,S} = LU{T,S}(copy(A.factors), copy(A.ipiv), A.info)
+copy(A::LU{T,S,P}) where {T,S,P} = LU{T,S,P}(copy(A.factors), copy(A.ipiv), A.info)
 
 size(A::LU)    = size(getfield(A, :factors))
 size(A::LU, i) = size(getfield(A, :factors), i)
@@ -564,7 +564,7 @@ function lu!(A::Tridiagonal{T,V}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(
     end
     B = Tridiagonal{T,V}(dl, d, du, du2)
     check && checknonsingular(info, pivot)
-    return LU{T,Tridiagonal{T,V}}(B, ipiv, convert(BlasInt, info))
+    return LU{T,Tridiagonal{T,V},typeof(ipiv)}(B, ipiv, convert(BlasInt, info))
 end
 
 factorize(A::Tridiagonal) = lu(A)
@@ -664,7 +664,7 @@ function ldiv!(transA::Transpose{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVec
 end
 
 # Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B)
-function ldiv!(adjA::Adjoint{<:Any,LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V}
+function ldiv!(adjA::Adjoint{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V}
     require_one_based_indexing(B)
     A = adjA.parent
     n = size(A,1)