Skip to content

Commit

Permalink
Generalize storage of factorizations.
Browse files Browse the repository at this point in the history
Introduce additional type parameters to remove hard-coded arrays from structures.
  • Loading branch information
maleadt committed Feb 11, 2022
1 parent 4409a1b commit 6e35e9f
Show file tree
Hide file tree
Showing 11 changed files with 191 additions and 153 deletions.
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 ⋅
Expand Down
25 changes: 14 additions & 11 deletions stdlib/LinearAlgebra/src/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
24 changes: 13 additions & 11 deletions stdlib/LinearAlgebra/src/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
33 changes: 18 additions & 15 deletions stdlib/LinearAlgebra/src/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
44 changes: 22 additions & 22 deletions stdlib/LinearAlgebra/src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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))
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 6e35e9f

Please sign in to comment.