Deprecate schur(...) in favor of schurfact(...) and factorization des…

…tructuring.

NEWS.md

@@ -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]).
 

stdlib/LinearAlgebra/docs/src/index.md

@@ -345,7 +345,6 @@ LinearAlgebra.hessfact
 LinearAlgebra.hessfact!
 LinearAlgebra.schurfact
 LinearAlgebra.schurfact!
-LinearAlgebra.schur
 LinearAlgebra.ordschur
 LinearAlgebra.ordschur!
 LinearAlgebra.svdfact

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -131,7 +131,6 @@ export
     lqfact,
     rank,
     rdiv!,
-    schur,
     schurfact!,
     schurfact,
     svd,

stdlib/LinearAlgebra/src/dense.jl

@@ -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
@@ -605,7 +605,7 @@ matrix function is returned whenever possible.
 If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is
 used, if `A` is triangular an improved version of the inverse scaling and squaring method is
 employed (see [^AH12] and [^AHR13]). For general matrices, the complex Schur form
-([`schur`](@ref)) is computed and the triangular algorithm is used on the
+([`schurfact`](@ref)) is computed and the triangular algorithm is used on the
 triangular factor.
 
 [^AH12]: Awad H. Al-Mohy and Nicholas J. Higham, "Improved inverse  scaling and squaring algorithms for the matrix logarithm", SIAM Journal on Scientific Computing, 34(4), 2012, C153-C169. [doi:10.1137/110852553](https://doi.org/10.1137/110852553)
@@ -662,7 +662,7 @@ that is the unique matrix ``X`` with eigenvalues having positive real part such
 
 If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is
 used to compute the square root. Otherwise, the square root is determined by means of the
-Björck-Hammarling method [^BH83], which computes the complex Schur form ([`schur`](@ref))
+Björck-Hammarling method [^BH83], which computes the complex Schur form ([`schurfact`](@ref))
 and then the complex square root of the triangular factor.
 
 [^BH83]:
@@ -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)
@@ -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)

stdlib/LinearAlgebra/src/deprecated.jl

@@ -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

stdlib/LinearAlgebra/src/schur.jl

@@ -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
 
@@ -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
@@ -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
 
@@ -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
 
@@ -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}
@@ -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:")

stdlib/LinearAlgebra/test/lapack.jl

@@ -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]
@@ -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]

stdlib/LinearAlgebra/test/schur.jl

@@ -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)
@@ -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)