Skip to content

Commit

Permalink
Deprecate schur(...) in favor of schurfact(...) and factorization des…
Browse files Browse the repository at this point in the history
…tructuring.
  • Loading branch information
Sacha0 committed May 13, 2018
1 parent dd01064 commit 1624676
Show file tree
Hide file tree
Showing 8 changed files with 155 additions and 66 deletions.
15 changes: 15 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -1025,6 +1025,21 @@ Deprecated or removed
(`vals, vecs = eigfact(A, B)`), or as a
`GeneralizedEigen` object (`eigf = eigfact(A, B)`) ([#26997]).

* `schur(A::AbstractMatrix)` has been deprecated in favor of `schurfact(A)`.
Whereas the former returns a tuple of arrays, the latter returns a `Schur` object.
So for a direct replacement, use `(schurfact(A)...,)`. But going forward, consider
using the direct result of `schurfact(A)` instead, either destructured into
its components (`T, Z, λ = schurfact(A)`) or as a `Schur` object
(`schurf = schurfact(A)`) ([#26997]).

* `schur(A::StridedMatrix, B::StridedMatrix)` has been deprecated in favor of
`schurfact(A, B)`. Whereas the former returns a tuple of arrays, the latter
returns a `GeneralizedSchur` object. So for a direct replacement, use
`(schurfact(A, B)...,)`. But going forward, consider using the direct result
of `schurfact(A, B)` instead, either destructured into its components
(`S, T, Q, Z, α, β = schurfact(A, B)`) or as a `GeneralizedSchur` object
(`schurf = schurfact(A, B)`) ([#26997]).

* The timing functions `tic`, `toc`, and `toq` are deprecated in favor of `@time` and `@elapsed`
([#17046]).

Expand Down
1 change: 0 additions & 1 deletion stdlib/LinearAlgebra/docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -345,7 +345,6 @@ LinearAlgebra.hessfact
LinearAlgebra.hessfact!
LinearAlgebra.schurfact
LinearAlgebra.schurfact!
LinearAlgebra.schur
LinearAlgebra.ordschur
LinearAlgebra.ordschur!
LinearAlgebra.svdfact
Expand Down
1 change: 0 additions & 1 deletion stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,6 @@ export
lqfact,
rank,
rdiv!,
schur,
schurfact!,
schurfact,
svd,
Expand Down
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -418,7 +418,7 @@ function schurpow(A::AbstractMatrix, p)
retmat = retmat * powm!(UpperTriangular(float.(A)), real(p - floor(p)))
end
else
S,Q,d = schur(complex(A))
S,Q,d = schurfact(complex(A))
# Integer part
R = S ^ floor(p)
# Real part
Expand Down Expand Up @@ -1409,8 +1409,8 @@ julia> A*X + X*B + C
```
"""
function sylvester(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) where T<:BlasFloat
RA, QA = schur(A)
RB, QB = schur(B)
RA, QA = schurfact(A)
RB, QB = schurfact(B)

D = -(adjoint(QA) * (C*QB))
Y, scale = LAPACK.trsyl!('N','N', RA, RB, D)
Expand Down Expand Up @@ -1453,7 +1453,7 @@ julia> A*X + X*A' + B
```
"""
function lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:BlasFloat}
R, Q = schur(A)
R, Q = schurfact(A)

D = -(adjoint(Q) * (C*Q))
Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
Expand Down
23 changes: 23 additions & 0 deletions stdlib/LinearAlgebra/src/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1322,3 +1322,26 @@ function _geneig(A, B)
"or as a `GeneralizedEigen` object (`eigf = eigfact(A, B)`)."), :eig)
return (eigfact(A, B)...,)
end

# deprecate schur(...) in favor of schurfact(...) and factorization destructuring
export schur
function schur(A::Union{StridedMatrix,Symmetric,Hermitian,UpperTriangular,LowerTriangular,Tridiagonal})
depwarn(string("`schur(A::AbstractMatrix)` has been deprecated in favor of ",
"`schurfact(A)`. Whereas the former returns a tuple of arrays, the ",
"latter returns a `Schur` object. So for a direct replacement, ",
"use `(schurfact(A)...,)`. But going forward, consider using ",
"the direct result of `schurfact(A)` instead, either destructured ",
"into its components (`T, Z, λ = schurfact(A)`) or as a `Schur` ",
"object (`schurf = schurfact(A)`)."), :schur)
return (schurfact(A)...,)
end
function schur(A::StridedMatrix, B::StridedMatrix)
depwarn(string("`schur(A::StridedMatrix, B::StridedMatrix)` has been ",
"deprecated in favor of `schurfact(A, B)`. Whereas the former returns ",
"a tuple of arrays, the latter returns a `GeneralizedSchur` object. ",
"So for a direct replacement, use `(schurfact(A, B)...,)`. But going ",
"forward, consider using the direct result of `schurfact(A, B)` instead, ",
"either destructured into its components (`S, T, Q, Z, α, β = schurfact(A, B)`) ",
" or as a `GeneralizedSchur` object (`schurf = schurfact(A, B)`)."), :schur)
return (schurfact(A, B)...,)
end
155 changes: 104 additions & 51 deletions stdlib/LinearAlgebra/src/schur.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,14 @@ struct Schur{Ty,S<:AbstractMatrix} <: Factorization{Ty}
end
Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) where {Ty} = Schur{Ty, typeof(T)}(T, Z, values)

# iteration for destructuring into factors
Base.start(::Schur) = Val(:T)
Base.next(F::Schur, ::Val{:T}) = (F.T, Val(:Z))
Base.next(F::Schur, ::Val{:Z}) = (F.Z, Val(:values))
Base.next(F::Schur, ::Val{:values}) = (F.values, Val(:done))
Base.done(F::Schur, ::Val{:done}) = true
Base.done(F::Schur, ::Any) = false

"""
schurfact!(A::StridedMatrix) -> F::Schur
Expand Down Expand Up @@ -78,10 +86,34 @@ julia> F.vectors * F.Schur * F.vectors'
2×2 Array{Float64,2}:
5.0 7.0
-2.0 -4.0
julia> T, Z, λ = schurfact(A); # destructuring via iteration
julia> T
2×2 Array{Float64,2}:
3.0 9.0
0.0 -2.0
julia> Z
2×2 Array{Float64,2}:
0.961524 0.274721
-0.274721 0.961524
julia> λ
2-element Array{Float64,1}:
3.0
-2.0
```
"""
schurfact(A::StridedMatrix{<:BlasFloat}) = schurfact!(copy(A))
schurfact(A::StridedMatrix{T}) where T = schurfact!(copy_oftype(A, eigtype(T)))
schurfact(A::Tridiagonal) = _schurfact_structured(A)
schurfact(A::Symmetric) = _schurfact_annotated(A)
schurfact(A::Hermitian) = _schurfact_annotated(A)
schurfact(A::UpperTriangular) = _schurfact_annotated(A)
schurfact(A::LowerTriangular) = _schurfact_annotated(A)
_schurfact_annotated(A) = schurfact(copyto!(similar(parent(A)), A))
_schurfact_structured(A) = schurfact(Matrix(A))

function getproperty(F::Schur, d::Symbol)
if d == :Schur
Expand All @@ -106,47 +138,6 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Schur)
show(io, mime, F.values)
end

"""
schur(A::StridedMatrix) -> T::Matrix, Z::Matrix, λ::Vector
Computes the Schur factorization of the matrix `A`. The methods return the (quasi)
triangular Schur factor `T` and the orthogonal/unitary Schur vectors `Z` such that
`A = Z * T * Z'`. The eigenvalues of `A` are returned in the vector `λ`.
See [`schurfact`](@ref).
# Examples
```jldoctest
julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
5.0 7.0
-2.0 -4.0
julia> T, Z, lambda = schur(A)
([3.0 9.0; 0.0 -2.0], [0.961524 0.274721; -0.274721 0.961524], [3.0, -2.0])
julia> Z * Z'
2×2 Array{Float64,2}:
1.0 0.0
0.0 1.0
julia> Z * T * Z'
2×2 Array{Float64,2}:
5.0 7.0
-2.0 -4.0
```
"""
function schur(A::StridedMatrix)
SchurF = schurfact(A)
SchurF.T, SchurF.Z, SchurF.values
end
schur(A::Symmetric) = schur(copyto!(similar(parent(A)), A))
schur(A::Hermitian) = schur(copyto!(similar(parent(A)), A))
schur(A::UpperTriangular) = schur(copyto!(similar(parent(A)), A))
schur(A::LowerTriangular) = schur(copyto!(similar(parent(A)), A))
schur(A::Tridiagonal) = schur(Matrix(A))


"""
ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur
Expand Down Expand Up @@ -209,6 +200,17 @@ function GeneralizedSchur(S::AbstractMatrix{Ty}, T::AbstractMatrix{Ty}, alpha::V
GeneralizedSchur{Ty, typeof(S)}(S, T, alpha, beta, Q, Z)
end

# iteration for destructuring into factors
Base.start(::GeneralizedSchur) = Val(:S)
Base.next(F::GeneralizedSchur, ::Val{:S}) = (F.S, Val(:T))
Base.next(F::GeneralizedSchur, ::Val{:T}) = (F.T, Val(:Q))
Base.next(F::GeneralizedSchur, ::Val{:Q}) = (F.Q, Val(:Z))
Base.next(F::GeneralizedSchur, ::Val{:Z}) = (F.Z, Val())
Base.next(F::GeneralizedSchur, ::Val{:α}) = (F.α, Val())
Base.next(F::GeneralizedSchur, ::Val{:β}) = (F.β, Val(:done))
Base.done(F::GeneralizedSchur, ::Val{:done}) = true
Base.done(F::GeneralizedSchur, ::Any) = false

"""
schurfact!(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur
Expand All @@ -226,6 +228,67 @@ and `F.T`, the left unitary/orthogonal Schur vectors can be obtained with `F.lef
`F.Q` and the right unitary/orthogonal Schur vectors can be obtained with `F.right` or
`F.Z` such that `A=F.left*F.S*F.right'` and `B=F.left*F.T*F.right'`. The
generalized eigenvalues of `A` and `B` can be obtained with `F.α./F.β`.
# Examples
```jldoctest
julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
5.0 7.0
-2.0 -4.0
julia> B = [6. 8.; -3. -5.]
2×2 Array{Float64,2}:
6.0 8.0
-3.0 -5.0
julia> F = schurfact(A, B)
GeneralizedSchur{Float64,Array{Float64,2}}
S factor:
2×2 Array{Float64,2}:
9.65794 0.323267
-0.504339 0.604369
T factor:
2×2 Array{Float64,2}:
11.5642 0.0
0.0 0.518842
Q factor:
2×2 Array{Float64,2}:
-0.864443 0.50273
0.50273 0.864443
Z factor:
2×2 Array{Float64,2}:
-0.578929 0.815378
-0.815378 -0.578929
α:
2-element Array{Complex{Float64},1}:
2.0000000000000044 + 2.738953629966988e-8im
2.9999999999999907 - 4.10843044495046e-8im
β:
2-element Array{Float64,1}:
2.0000000000000058
2.9999999999999925
julia> F.Q * F.S * F.Z'
2×2 Array{Float64,2}:
5.0 7.0
-2.0 -4.0
julia> F.Q * F.T * F.Z'
2×2 Array{Float64,2}:
6.0 8.0
-3.0 -5.0
julia> S, T, Q, Z, α, β = schurfact(A, B); # destructuring via iteration
julia> Q * S * Z' ≈ A
true
julia> Q * T * Z' ≈ B
true
julia> α. / β ≈ eigvals(A, B)
true
```
"""
schurfact(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:BlasFloat} = schurfact!(copy(A),copy(B))
function schurfact(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
Expand Down Expand Up @@ -300,16 +363,6 @@ end
Base.propertynames(F::GeneralizedSchur) =
(:values, :left, :right, fieldnames(typeof(F))...)

"""
schur(A::StridedMatrix, B::StridedMatrix) -> S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, α::Vector, β::Vector
See [`schurfact`](@ref).
"""
function schur(A::StridedMatrix, B::StridedMatrix)
SchurF = schurfact(A, B)
SchurF.S, SchurF.T, SchurF.Q, SchurF.Z, SchurF.α, SchurF.β
end

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::GeneralizedSchur)
println(io, summary(F))
println(io, "S factor:")
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/test/lapack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -590,7 +590,7 @@ end
for job in ('N', 'E', 'V', 'B')
for c in ('V', 'N')
A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
T,Q,d = schur(A)
T,Q,d = schurfact(A)
s, sep = LinearAlgebra.LAPACK.trsen!(job,c,Array{LinearAlgebra.BlasInt}([0,1,0,0]),T,Q)[4:5]
@test d[1] T[2,2]
@test d[2] T[1,1]
Expand All @@ -616,7 +616,7 @@ end
@testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
for c in ('V', 'N')
A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
T,Q,d = schur(A)
T,Q,d = schurfact(A)
LinearAlgebra.LAPACK.trexc!(c,LinearAlgebra.BlasInt(1),LinearAlgebra.BlasInt(2),T,Q)
@test d[1] T[2,2]
@test d[2] T[1,1]
Expand Down
14 changes: 7 additions & 7 deletions stdlib/LinearAlgebra/test/schur.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,15 @@ aimg = randn(n,n)/2
@test convert(Array, f) a
@test_throws ErrorException f.A

sch, vecs, vals = schur(UpperTriangular(triu(a)))
sch, vecs, vals = schurfact(UpperTriangular(triu(a)))
@test vecs*sch*vecs' triu(a)
sch, vecs, vals = schur(LowerTriangular(tril(a)))
sch, vecs, vals = schurfact(LowerTriangular(tril(a)))
@test vecs*sch*vecs' tril(a)
sch, vecs, vals = schur(Hermitian(asym))
sch, vecs, vals = schurfact(Hermitian(asym))
@test vecs*sch*vecs' asym
sch, vecs, vals = schur(Symmetric(a + transpose(a)))
sch, vecs, vals = schurfact(Symmetric(a + transpose(a)))
@test vecs*sch*vecs' a + transpose(a)
sch, vecs, vals = schur(Tridiagonal(a + transpose(a)))
sch, vecs, vals = schurfact(Tridiagonal(a + transpose(a)))
@test vecs*sch*vecs' Tridiagonal(a + transpose(a))

tstring = sprint((t, s) -> show(t, "text/plain", s), f.T)
Expand Down Expand Up @@ -111,15 +111,15 @@ aimg = randn(n,n)/2
@test S.Z SchurNew.Z
@test S.alpha SchurNew.alpha
@test S.beta SchurNew.beta
sS,sT,sQ,sZ = schur(a1_sf,a2_sf)
sS,sT,sQ,sZ = schurfact(a1_sf,a2_sf)
@test NS.Q sQ
@test NS.T sT
@test NS.S sS
@test NS.Z sZ
end
end
@testset "0x0 matrix" for A in (zeros(eltya, 0, 0), view(rand(eltya, 2, 2), 1:0, 1:0))
T, Z, λ = LinearAlgebra.schur(A)
T, Z, λ = LinearAlgebra.schurfact(A)
@test T == A
@test Z == A
@test λ == zeros(0)
Expand Down

0 comments on commit 1624676

Please sign in to comment.