Skip to content

Commit

Permalink
use singleton types
Browse files Browse the repository at this point in the history
  • Loading branch information
dkarrasch committed May 20, 2021
1 parent b8d5251 commit 4e8487b
Show file tree
Hide file tree
Showing 14 changed files with 83 additions and 100 deletions.
6 changes: 3 additions & 3 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,9 @@ Standard library changes
* The shape of an `UpperHessenberg` matrix is preserved under certain arithmetic operations, e.g. when multiplying or dividing by an `UpperTriangular` matrix. ([#40039])
* `cis(A)` now supports matrix arguments ([#40194]).
* `dot` now supports `UniformScaling` with `AbstractMatrix` ([#40250]).
* `qr[!]` and `lu[!]` now support `Symbol` values as their optional `pivot` argument:
defaults are `qr(A, :none)` (vs. `qr(A, :colnorm)` for pivoting) and `lu(A, :rowmax)`
(vs. `lu(A, :none)` without pivoting); the former `Val{true/false}`-based calls are deprecated. ([#40623])
* `qr[!]` and `lu[!]` now support `PivotingStrategy` values as their optional `pivot` argument:
defaults are `qr(A, NoPivot())` (vs. `qr(A, ColNorm())` for pivoting) and `lu(A, RowMax())`
(vs. `lu(A, NoPivot())` without pivoting); the former `Val{true/false}`-based calls are deprecated. ([#40623])

#### Markdown

Expand Down
7 changes: 7 additions & 0 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,19 +35,22 @@ export
BunchKaufman,
Cholesky,
CholeskyPivoted,
ColNorm,
Eigen,
GeneralizedEigen,
GeneralizedSVD,
GeneralizedSchur,
Hessenberg,
LU,
LDLt,
NoPivot,
QR,
QRPivoted,
LQ,
Schur,
SVD,
Hermitian,
RowMax,
Symmetric,
LowerTriangular,
UpperTriangular,
Expand Down Expand Up @@ -164,6 +167,10 @@ abstract type Algorithm end
struct DivideAndConquer <: Algorithm end
struct QRIteration <: Algorithm end

abstract type PivotingStrategy end
struct NoPivot <: PivotingStrategy end
struct RowMax <: PivotingStrategy end
struct ColNorm <: PivotingStrategy end

# Check that stride of matrix/vector is 1
# Writing like this to avoid splatting penalty when called with multiple arguments,
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1371,7 +1371,7 @@ function factorize(A::StridedMatrix{T}) where T
end
return lu(A)
end
qr(A, :colnorm)
qr(A, ColNorm())
end
factorize(A::Adjoint) = adjoint(factorize(parent(A)))
factorize(A::Transpose) = transpose(factorize(parent(A)))
Expand Down
11 changes: 3 additions & 8 deletions stdlib/LinearAlgebra/src/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,9 @@ size(F::Adjoint{<:Any,<:Factorization}) = reverse(size(parent(F)))
size(F::Transpose{<:Any,<:Factorization}) = reverse(size(parent(F)))

checkpositivedefinite(info) = info == 0 || throw(PosDefException(info))
function checknonsingular(info, pivoted = :rowmax)
if info != 0
pivoted === :rowmax && throw(SingularException(info))
pivoted === :none && throw(ZeroPivotException(info))
else
return nothing
end
end
checknonsingular(info, ::RowMax) = info == 0 || throw(SingularException(info))
checknonsingular(info, ::NoPivot) = info == 0 || throw(ZeroPivotException(info))
checknonsingular(info) = checknonsingular(info, RowMax())

"""
issuccess(F::Factorization)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1141,7 +1141,7 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat)
end
return lu(A) \ B
end
return qr(A, :colnorm) \ B
return qr(A, ColNorm()) \ B
end

(\)(a::AbstractVector, b::AbstractArray) = pinv(a) * b
Expand Down
63 changes: 30 additions & 33 deletions stdlib/LinearAlgebra/src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -76,28 +76,26 @@ adjoint(F::LU) = Adjoint(F)
transpose(F::LU) = Transpose(F)

# StridedMatrix
function lu!(A::StridedMatrix{T}, pivot::Symbol = :rowmax; check::Bool = true) where {T<:BlasFloat}
if pivot === :none
return generic_lufact!(A, pivot; check = check)
elseif pivot === :rowmax
lpt = LAPACK.getrf!(A)
check && checknonsingular(lpt[3])
return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
else
throw(ArgumentError("only `:rowmax` and `:none` are supported as `pivot` argument but you supplied `$pivot`"))
end
lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMax(); check=check)
function lu!(A::StridedMatrix{T}, ::RowMax; 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])
end
function lu!(A::HermOrSym, pivot::Symbol = :rowmax; check::Bool = true)
function lu!(A::StridedMatrix{<:BlasFloat}, pivot::NoPivot; check::Bool = true)
return generic_lufact!(A, pivot; check = check)
end
function lu!(A::HermOrSym, pivot::PivotingStrategy = RowMax(); check::Bool = true)
copytri!(A.data, A.uplo, isa(A, Hermitian))
lu!(A.data, pivot; check = check)
end
# for backward compatibility
# TODO: remove towards Julia v2
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, :rowmax; check=check)
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{false}; check::Bool = true) lu!(A, :none; check=check)
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, RowMax(); check=check)
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{false}; check::Bool = true) lu!(A, NoPivot(); check=check)

"""
lu!(A, pivot = :rowmax; check = true) -> LU
lu!(A, pivot = RowMax(); check = true) -> LU
`lu!` is the same as [`lu`](@ref), but saves space by overwriting the
input `A`, instead of creating a copy. An [`InexactError`](@ref)
Expand Down Expand Up @@ -133,26 +131,26 @@ Stacktrace:
[...]
```
"""
lu!(A::StridedMatrix, pivot::Symbol = :rowmax; check::Bool = true) =
lu!(A::StridedMatrix, pivot::PivotingStrategy = RowMax(); check::Bool = true) =
generic_lufact!(A, pivot; check = check)
function generic_lufact!(A::StridedMatrix{T}, pivot = :rowmax; check::Bool = true) where T
# Check arguments
if pivot !== :rowmax && pivot !== :none
throw(ArgumentError("only `rowmax` and `none` are supported as `pivot` argument but you supplied `$pivot`"))
end

function generic_lufact!(A::StridedMatrix{T}, pivot::PivotingStrategy = RowMax();
check::Bool = true) where {T}
# Extract values
m, n = size(A)
minmn = min(m,n)

if pivot !== RowMax() && pivot !== NoPivot()
throw(ArgumentError("only `RowMax()` and `NoPivot()` are supported as `pivot` argument but you supplied `$pivot`"))
end

# Initialize variables
info = 0
ipiv = Vector{BlasInt}(undef, minmn)
@inbounds begin
for k = 1:minmn
# find index max
kp = k
if pivot === :rowmax && k < m
if pivot === RowMax() && k < m
amax = abs(A[k, k])
for i = k+1:m
absi = abs(A[i,k])
Expand Down Expand Up @@ -213,7 +211,7 @@ end

# for all other types we must promote to a type which is stable under division
"""
lu(A, pivot = :rowmax; check = true) -> F::LU
lu(A, pivot = RowMax(); check = true) -> F::LU
Compute the LU factorization of `A`.
Expand All @@ -224,7 +222,7 @@ validity (via [`issuccess`](@ref)) lies with the user.
In most cases, if `A` is a subtype `S` of `AbstractMatrix{T}` with an element
type `T` supporting `+`, `-`, `*` and `/`, the return type is `LU{T,S{T}}`. If
pivoting is chosen (default) the element type should also support [`abs`](@ref) and
[`<`](@ref). Pivoting can be turned off by passing `pivot = :none`.
[`<`](@ref). Pivoting can be turned off by passing `pivot = NoPivot()`.
The individual components of the factorization `F` can be accessed via [`getproperty`](@ref):
Expand Down Expand Up @@ -280,13 +278,13 @@ julia> l == F.L && u == F.U && p == F.p
true
```
"""
function lu(A::AbstractMatrix{T}, pivot::Symbol = :rowmax; check::Bool = true) where {T}
function lu(A::AbstractMatrix{T}, pivot::PivotingStrategy = RowMax(); check::Bool = true) where {T}
S = lutype(T)
lu!(copy_oftype(A, S), pivot; check = check)
end
# TODO: remove for Julia v2.0
@deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, :rowmax; check=check)
@deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, :none; check=check)
@deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, RowMax(); check=check)
@deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, NoPivot(); check=check)


lu(S::LU) = S
Expand Down Expand Up @@ -497,13 +495,12 @@ inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A))
# Tridiagonal

# See dgttrf.f
function lu!(A::Tridiagonal{T,V}, pivot::Symbol = :rowmax; check::Bool = true) where {T,V}
function lu!(A::Tridiagonal{T,V}, pivot::PivotingStrategy = RowMax(); check::Bool = true) where {T,V}
# Extract values
n = size(A, 1)

# Check arguments
if pivot !== :rowmax && pivot !== :none
throw(ArgumentError("only `:row` and `:none` are supported as `pivot` argument but you supplied `$pivot`"))
if pivot !== RowMax() && pivot !== NoPivot()
throw(ArgumentError("only `RowMax()` and `NoPivot()` are supported as `pivot` argument but you supplied `$pivot`"))
end

# Initialize variables
Expand All @@ -523,7 +520,7 @@ function lu!(A::Tridiagonal{T,V}, pivot::Symbol = :rowmax; check::Bool = true) w
end
for i = 1:n-2
# pivot or not?
if pivot === :none || abs(d[i]) >= abs(dl[i])
if pivot === NoPivot() || abs(d[i]) >= abs(dl[i])
# No interchange
if d[i] != 0
fact = dl[i]/d[i]
Expand All @@ -546,7 +543,7 @@ function lu!(A::Tridiagonal{T,V}, pivot::Symbol = :rowmax; check::Bool = true) w
end
if n > 1
i = n-1
if pivot === :none || abs(d[i]) >= abs(dl[i])
if pivot === NoPivot() || abs(d[i]) >= abs(dl[i])
if d[i] != 0
fact = dl[i]/d[i]
dl[i] = fact
Expand Down
42 changes: 15 additions & 27 deletions stdlib/LinearAlgebra/src/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -246,20 +246,14 @@ function qrfactPivotedUnblocked!(A::AbstractMatrix)
end

# LAPACK version
Base.@aggressive_constprop function qr!(A::StridedMatrix{<:BlasFloat}, pivot::Symbol = :none; blocksize=36)
if pivot === :none
return QRCompactWY(LAPACK.geqrt!(A, min(min(size(A)...), blocksize))...)
elseif pivot === :colnorm
return QRPivoted(LAPACK.geqp3!(A)...)
else
throw(ArgumentError("only `:colnorm` and `:none` are supported as `pivot` argument but you supplied `$pivot`"))
end
end
qr!(A::StridedMatrix{<:BlasFloat}, ::NoPivot; blocksize=36) =
QRCompactWY(LAPACK.geqrt!(A, min(min(size(A)...), blocksize))...)
qr!(A::StridedMatrix{<:BlasFloat}, ::ColNorm) = QRPivoted(LAPACK.geqp3!(A)...)

# Generic fallbacks

"""
qr!(A, pivot = :none; blocksize)
qr!(A, pivot = NoPivot(); blocksize)
`qr!` is the same as [`qr`](@ref) when `A` is a subtype of [`StridedMatrix`](@ref),
but saves space by overwriting the input `A`, instead of creating a copy.
Expand Down Expand Up @@ -298,23 +292,17 @@ Stacktrace:
[...]
```
"""
Base.@aggressive_constprop function qr!(A::AbstractMatrix, pivot::Symbol = :none)
if pivot === :none
return qrfactUnblocked!(A)
elseif pivot === :colnorm
return qrfactPivotedUnblocked!(A)
else
throw(ArgumentError("only `:colnorm` and `:none` are supported as `pivot` argument but you supplied `$pivot`"))
end
end
qr!(A::AbstractMatrix, ::NoPivot) = qrfactUnblocked!(A)
qr!(A::AbstractMatrix, ::ColNorm) = qrfactPivotedUnblocked!(A)
qr!(A::AbstractMatrix) = qr!(A, NoPivot())
# TODO: Remove in Julia v2.0
@deprecate qr!(A::AbstractMatrix, ::Val{true}) qr!(A, :colnorm)
@deprecate qr!(A::AbstractMatrix, ::Val{false}) qr!(A, :none)
@deprecate qr!(A::AbstractMatrix, ::Val{true}) qr!(A, ColNorm())
@deprecate qr!(A::AbstractMatrix, ::Val{false}) qr!(A, NoPivot())

_qreltype(::Type{T}) where T = typeof(zero(T)/sqrt(abs2(one(T))))

"""
qr(A, pivot = :none; blocksize) -> F
qr(A, pivot = NoPivot(); blocksize) -> F
Compute the QR factorization of the matrix `A`: an orthogonal (or unitary if `A` is
complex-valued) matrix `Q`, and an upper triangular matrix `R` such that
Expand All @@ -325,7 +313,7 @@ A = Q R
The returned object `F` stores the factorization in a packed format:
- if `pivot == :colnorm` then `F` is a [`QRPivoted`](@ref) object,
- if `pivot == ColNorm()` then `F` is a [`QRPivoted`](@ref) object,
- otherwise if the element type of `A` is a BLAS type ([`Float32`](@ref), [`Float64`](@ref),
`ComplexF32` or `ComplexF64`), then `F` is a [`QRCompactWY`](@ref) object,
Expand Down Expand Up @@ -355,7 +343,7 @@ and `F.Q*A` are supported. A `Q` matrix can be converted into a regular matrix w
orthogonal matrix.
The block size for QR decomposition can be specified by keyword argument
`blocksize :: Integer` when `pivot == :none` and `A isa StridedMatrix{<:BlasFloat}`.
`blocksize :: Integer` when `pivot == NoPivot()` and `A isa StridedMatrix{<:BlasFloat}`.
It is ignored when `blocksize > minimum(size(A))`. See [`QRCompactWY`](@ref).
!!! compat "Julia 1.4"
Expand Down Expand Up @@ -391,15 +379,15 @@ true
elementary reflectors, so that the `Q` and `R` matrices can be stored
compactly rather as two separate dense matrices.
"""
Base.@aggressive_constprop function qr(A::AbstractMatrix{T}, arg...; kwargs...) where T
function qr(A::AbstractMatrix{T}, arg...; kwargs...) where T
require_one_based_indexing(A)
AA = similar(A, _qreltype(T), size(A))
copyto!(AA, A)
return qr!(AA, arg...; kwargs...)
end
# TODO: remove in Julia v2.0
@deprecate qr(A::AbstractMatrix, ::Val{false}; kwargs...) qr(A, :none; kwargs...)
@deprecate qr(A::AbstractMatrix, ::Val{true}; kwargs...) qr(A, :colnorm; kwargs...)
@deprecate qr(A::AbstractMatrix, ::Val{false}; kwargs...) qr(A, NoPivot(); kwargs...)
@deprecate qr(A::AbstractMatrix, ::Val{true}; kwargs...) qr(A, ColNorm(); kwargs...)

qr(x::Number) = qr(fill(x,1,1))
function qr(v::AbstractVector)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/test/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -565,7 +565,7 @@ end
D = Diagonal(randn(5))
Q = qr(randn(5, 5)).Q
@test D * Q' == Array(D) * Q'
Q = qr(randn(5, 5), :colnorm).Q
Q = qr(randn(5, 5), ColNorm()).Q
@test_throws ArgumentError lmul!(Q, D)
end

Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/test/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -382,13 +382,13 @@ LinearAlgebra.Transpose(a::ModInt{n}) where {n} = transpose(a)
A = [ModInt{2}(1) ModInt{2}(0); ModInt{2}(1) ModInt{2}(1)]
b = [ModInt{2}(1), ModInt{2}(0)]

@test A*(lu(A, :none)\b) == b
@test A*(lu(A, NoPivot())\b) == b

# Needed for pivoting:
Base.abs(a::ModInt{n}) where {n} = a
Base.:<(a::ModInt{n}, b::ModInt{n}) where {n} = a.k < b.k

@test A*(lu(A, :rowmax)\b) == b
@test A*(lu(A, RowMax())\b) == b
end

@testset "Issue 18742" begin
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/test/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
lqa = lq(a)
x = lqa\b
l,q = lqa.L, lqa.Q
qra = qr(a, :colnorm)
qra = qr(a, ColNorm())
@testset "Basic ops" begin
@test size(lqa,1) == size(a,1)
@test size(lqa,3) == 1
Expand Down
16 changes: 8 additions & 8 deletions stdlib/LinearAlgebra/test/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ dimg = randn(n)/2
end
κd = cond(Array(d),1)
@testset "Tridiagonal LU" begin
lud = @inferred lu(d)
lud = @inferred lu(d)
@test LinearAlgebra.issuccess(lud)
@test @inferred(lu(lud)) == lud
@test_throws ErrorException lud.Z
Expand Down Expand Up @@ -229,12 +229,12 @@ end
@test_throws SingularException lu!(copy(A); check = true)
@test !issuccess(lu(A; check = false))
@test !issuccess(lu!(copy(A); check = false))
@test_throws ZeroPivotException lu(A, :none)
@test_throws ZeroPivotException lu!(copy(A), :none)
@test_throws ZeroPivotException lu(A, :none; check = true)
@test_throws ZeroPivotException lu!(copy(A), :none; check = true)
@test !issuccess(lu(A, :none; check = false))
@test !issuccess(lu!(copy(A), :none; check = false))
@test_throws ZeroPivotException lu(A, NoPivot())
@test_throws ZeroPivotException lu!(copy(A), NoPivot())
@test_throws ZeroPivotException lu(A, NoPivot(); check = true)
@test_throws ZeroPivotException lu!(copy(A), NoPivot(); check = true)
@test !issuccess(lu(A, NoPivot(); check = false))
@test !issuccess(lu!(copy(A), NoPivot(); check = false))
F = lu(A; check = false)
@test sprint((io, x) -> show(io, "text/plain", x), F) ==
"Failed factorization of type $(typeof(F))"
Expand Down Expand Up @@ -320,7 +320,7 @@ include("trickyarithmetic.jl")
@testset "lu with type whose sum is another type" begin
A = TrickyArithmetic.A[1 2; 3 4]
ElT = TrickyArithmetic.D{TrickyArithmetic.C,TrickyArithmetic.C}
B = lu(A, :none)
B = lu(A, NoPivot())
@test B isa LinearAlgebra.LU{ElT,Matrix{ElT}}
end

Expand Down
Loading

0 comments on commit 4e8487b

Please sign in to comment.