From 76cc0344cab1a343116c6dd50bfc32c6eea95035 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 7 Jan 2022 13:08:30 +0100 Subject: [PATCH 1/4] Cleanup `copy_oftype`-like functions in factorizations --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 21 ++++----------------- stdlib/LinearAlgebra/src/lu.jl | 6 +----- stdlib/LinearAlgebra/src/schur.jl | 13 ++++--------- 3 files changed, 9 insertions(+), 31 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 8afe9905e9a30..21344c9b7c0a2 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -366,11 +366,11 @@ In general, the type of the output corresponds to that of `similar(A, T)`. There are three often used methods in LinearAlgebra to create a mutable copy of an array with a given eltype. These copies can be passed to in-place -algorithms (such as ldiv!, rdiv!, lu! and so on). Which one to use in practice +algorithms (such as `ldiv!`, `rdiv!`, `lu!` and so on). Which one to use in practice depends on what is known (or assumed) about the structure of the array in that algorithm. -See also: `copy_similar`, `copy_to_array`. +See also: `copy_similar`. """ copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T), A) @@ -380,25 +380,12 @@ copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T), A) Copy `A` to a mutable array with eltype `T` based on `similar(A, T, size(A))`. Compared to `copy_oftype`, the result can be more flexible. In general, the type -of the output corresponds to that of the three-argument method `similar(A, T, size(s))`. +of the output corresponds to that of the three-argument method `similar(A, T, size(A))`. -See also: `copy_oftype`, `copy_to_array`. +See also: `copy_oftype`. """ copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A) -""" - copy_to_array(A, T) - -Copy `A` to a regular dense `Array` with element type `T`. - -The resulting array is mutable. It can be used, for example, to pass the data of -`A` to an efficient in-place method for a matrix factorization such as `lu!`, in -cases where a more specific implementation of `lu!` (or `lu`) is not available. - -See also: `copy_oftype`, `copy_similar`. -""" -copy_to_array(A::AbstractArray, ::Type{T}) where {T} = copyto!(Array{T}(undef, size(A)...), A) - # The three copy functions above return mutable arrays with eltype T. # To only ensure a certain eltype, and if a mutable copy is not needed, it is # more efficient to use: diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 22ed1fc5a49dc..a6b0b7a3f57a5 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -76,9 +76,6 @@ adjoint(F::LU) = Adjoint(F) transpose(F::LU) = Transpose(F) # StridedMatrix -lu(A::StridedMatrix, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) = - lu!(copy_oftype(A, lutype(eltype(A))), pivot; check=check) - lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMaximum(); check=check) function lu!(A::StridedMatrix{T}, ::RowMaximum; check::Bool = true) where {T<:BlasFloat} lpt = LAPACK.getrf!(A) @@ -282,8 +279,7 @@ true ``` """ function lu(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T} - S = lutype(T) - lu!(copy_to_array(A, S), pivot; check = check) + lu!(copy_similar(A, lutype(T)), pivot; check = check) end # TODO: remove for Julia v2.0 @deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, RowMaximum(); check=check) diff --git a/stdlib/LinearAlgebra/src/schur.jl b/stdlib/LinearAlgebra/src/schur.jl index 610067fe51452..2db1bbb879c1b 100644 --- a/stdlib/LinearAlgebra/src/schur.jl +++ b/stdlib/LinearAlgebra/src/schur.jl @@ -97,7 +97,7 @@ julia> A schur!(A::StridedMatrix{<:BlasFloat}) = Schur(LinearAlgebra.LAPACK.gees!('V', A)...) """ - schur(A::StridedMatrix) -> F::Schur + schur(A::AbstractMatrix) -> F::Schur Computes the Schur factorization of the matrix `A`. The (quasi) triangular Schur factor can be obtained from the `Schur` object `F` with either `F.Schur` or `F.T` and the @@ -146,25 +146,20 @@ julia> t == F.T && z == F.Z && vals == F.values true ``` """ -schur(A::StridedMatrix{<:BlasFloat}) = schur!(copy(A)) -schur(A::StridedMatrix{T}) where T = schur!(copy_oftype(A, eigtype(T))) - -schur(A::AbstractMatrix{T}) where {T} = schur!(copy_to_array(A, eigtype(T))) +schur(A::AbstractMatrix{T}) where {T} = schur!(copy_similar(A, eigtype(T))) function schur(A::RealHermSymComplexHerm) F = eigen(A; sortby=nothing) return Schur(typeof(F.vectors)(Diagonal(F.values)), F.vectors, F.values) end function schur(A::Union{UnitUpperTriangular{T},UpperTriangular{T}}) where {T} t = eigtype(T) - Z = Matrix{t}(undef, size(A)...) - copyto!(Z, A) + Z = copy_similar(A, t) return Schur(Z, Matrix{t}(I, size(A)), convert(Vector{t}, diag(A))) end function schur(A::Union{UnitLowerTriangular{T},LowerTriangular{T}}) where {T} t = eigtype(T) # double flip the matrix A - Z = Matrix{t}(undef, size(A)...) - copyto!(Z, A) + Z = copy_similar(A, t) reverse!(reshape(Z, :)) # construct "reverse" identity n = size(A, 1) From 122f4e1ddc508c254b0dfac4f0df3f0e8cf3e86f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 7 Jan 2022 16:24:37 +0100 Subject: [PATCH 2/4] further reduce number of `lu` methods Co-authored-by: @daanhb --- stdlib/LinearAlgebra/src/lu.jl | 11 ++++------- stdlib/LinearAlgebra/src/schur.jl | 11 +++-------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index a6b0b7a3f57a5..bae26a5868aaf 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -86,9 +86,6 @@ function lu!(A::StridedMatrix{<:BlasFloat}, pivot::NoPivot; check::Bool = true) return generic_lufact!(A, pivot; check = check) end -lu(A::HermOrSym, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) = - lu!(copy_oftype(A, lutype(eltype(A))), pivot; check=check) - function lu!(A::HermOrSym, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) copytri!(A.data, A.uplo, isa(A, Hermitian)) lu!(A.data, pivot; check = check) @@ -279,12 +276,15 @@ true ``` """ function lu(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T} - lu!(copy_similar(A, lutype(T)), pivot; check = check) + lu!(_lucopy(A, lutype(T)), pivot; check = check) end # TODO: remove for Julia v2.0 @deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, RowMaximum(); check=check) @deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, NoPivot(); check=check) +_lucopy(A::AbstractMatrix, T) = copy_similar(A, T) +_lucopy(A::HermOrSym, T) = copy_oftype(A, T) +_lucopy(A::Tridiagonal, T) = copy_oftype(A, T) lu(S::LU) = S function lu(x::Number; check::Bool=true) @@ -493,9 +493,6 @@ inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A)) # Tridiagonal -lu(A::Tridiagonal{T}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where T = - lu!(copy_oftype(A, lutype(T)), pivot; check = check) - # See dgttrf.f function lu!(A::Tridiagonal{T,V}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T,V} # Extract values diff --git a/stdlib/LinearAlgebra/src/schur.jl b/stdlib/LinearAlgebra/src/schur.jl index 2db1bbb879c1b..549ea5cdf3d3a 100644 --- a/stdlib/LinearAlgebra/src/schur.jl +++ b/stdlib/LinearAlgebra/src/schur.jl @@ -97,7 +97,7 @@ julia> A schur!(A::StridedMatrix{<:BlasFloat}) = Schur(LinearAlgebra.LAPACK.gees!('V', A)...) """ - schur(A::AbstractMatrix) -> F::Schur + schur(A) -> F::Schur Computes the Schur factorization of the matrix `A`. The (quasi) triangular Schur factor can be obtained from the `Schur` object `F` with either `F.Schur` or `F.T` and the @@ -333,7 +333,7 @@ schur!(A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat} = GeneralizedSchur(LinearAlgebra.LAPACK.gges!('V', 'V', A, B)...) """ - schur(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur + schur(A, B) -> F::GeneralizedSchur Computes the Generalized Schur (or QZ) factorization of the matrices `A` and `B`. The (quasi) triangular Schur factors can be obtained from the `Schur` object `F` with `F.S` @@ -345,14 +345,9 @@ generalized eigenvalues of `A` and `B` can be obtained with `F.α./F.β`. Iterating the decomposition produces the components `F.S`, `F.T`, `F.Q`, `F.Z`, `F.α`, and `F.β`. """ -schur(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:BlasFloat} = schur!(copy(A),copy(B)) -function schur(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} - S = promote_type(eigtype(TA), TB) - return schur!(copy_oftype(A, S), copy_oftype(B, S)) -end function schur(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB} S = promote_type(eigtype(TA), TB) - return schur!(copy_oftype(A, S), copy_oftype(B, S)) + return schur!(copy_similar(A, S), copy_similar(B, S)) end """ From e26db9b7c687650611a4086590ccac6aeda2d292 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 8 Jan 2022 19:37:20 +0100 Subject: [PATCH 3/4] extend specialized lu to *diagonal --- stdlib/LinearAlgebra/src/bidiag.jl | 2 +- stdlib/LinearAlgebra/src/special.jl | 9 ++++++++ stdlib/LinearAlgebra/test/lu.jl | 36 ++++++++++++++++++----------- 3 files changed, 32 insertions(+), 15 deletions(-) diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index a4f9ac468cd34..01580c9246b09 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/stdlib/LinearAlgebra/src/bidiag.jl @@ -71,7 +71,7 @@ end #To allow Bidiagonal's where the "dv" is Vector{T} and "ev" Vector{S}, #where T and S can be promoted -function LinearAlgebra.Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Symbol) where {T,S} +function Bidiagonal(dv::Vector{T}, ev::Vector{S}, uplo::Symbol) where {T,S} TS = promote_type(T,S) return Bidiagonal{TS,Vector{TS}}(dv, ev, uplo) end diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 9674599b4e334..d884636b43720 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -50,6 +50,15 @@ Bidiagonal(A::AbstractTriangular) = isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal throw(ArgumentError("matrix cannot be represented as Bidiagonal")) +_lucopy(A::Bidiagonal, T) = copy_oftype(Tridiagonal(A), T) +_lucopy(A::Diagonal, T) = copy_oftype(Tridiagonal(A), T) +function _lucopy(A::SymTridiagonal, T) + du = copy_similar(_evview(A), T) + dl = copy.(transpose.(du)) + d = copy_similar(A.dv, T) + return Tridiagonal(dl, d, du) +end + const ConvertibleSpecialMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,AbstractTriangular} const PossibleTriangularMatrix = Union{Diagonal, Bidiagonal, AbstractTriangular} diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index cc3f1be2d1627..4ac49a9cd4e1f 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -402,20 +402,28 @@ end end end -@testset "lu(A) has a fallback for abstract matrices (#40831)" begin - # check that lu works for some structured arrays - A0 = rand(5, 5) - @test lu(Diagonal(A0)) isa LU - @test Matrix(lu(Diagonal(A0))) ≈ Diagonal(A0) - @test lu(Bidiagonal(A0, :U)) isa LU - @test Matrix(lu(Bidiagonal(A0, :U))) ≈ Bidiagonal(A0, :U) - - # lu(A) copies A and then invokes lu!, make sure that the most efficient - # implementation of lu! continues to be used - A1 = Tridiagonal(rand(2), rand(3), rand(2)) - @test lu(A1) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} - @test lu(A1, RowMaximum()) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} - @test lu(A1, RowMaximum(); check = false) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} +@testset "lu on *diagonal matrices" begin + dl = rand(3) + d = rand(4) + Bl = Bidiagonal(d, dl, :L) + Bu = Bidiagonal(d, dl, :U) + Tri = Tridiagonal(dl, d, dl) + Sym = SymTridiagonal(d, dl) + D = Diagonal(d) + b = ones(4) + B = rand(4,4) + for A in (Bl, Bu, Tri, Sym, D), pivot in (NoPivot(), RowMaximum()) + @test A\b ≈ lu(A, pivot)\b + @test B/A ≈ B/lu(A, pivot) + @test B/A ≈ B/Matrix(A) + @test Matrix(lu(A, pivot)) ≈ A + @test @inferred(lu(A)) isa LU + if A isa Union{Bidiagonal, Diagonal, Tridiagonal, SymTridiagonal} + @test lu(A) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} + @test lu(A, pivot) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} + @test lu(A, pivot; check = false) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} + end + end end end # module TestLU From 2bbde96fc3bd79b2606068b010249cecc4684525 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 13 Jan 2022 15:32:33 +0100 Subject: [PATCH 4/4] whitespace --- stdlib/LinearAlgebra/src/special.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index d884636b43720..e876e40c8065d 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -55,10 +55,10 @@ _lucopy(A::Diagonal, T) = copy_oftype(Tridiagonal(A), T) function _lucopy(A::SymTridiagonal, T) du = copy_similar(_evview(A), T) dl = copy.(transpose.(du)) - d = copy_similar(A.dv, T) + d = copy_similar(A.dv, T) return Tridiagonal(dl, d, du) end - + const ConvertibleSpecialMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,AbstractTriangular} const PossibleTriangularMatrix = Union{Diagonal, Bidiagonal, AbstractTriangular}