Implement schur_full for Diagonal inputs, deprecate schur_vals

docs/src/changelog.md

@@ -33,10 +33,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Deprecated
 
+- `schur_vals` and `schur_vals!` in favour of `eig_vals` and `eig_vals!`. LAPACK's `gees` balances
+  without scaling, so the eigenvalues it returns are those of `eig_vals(A; scale = false)`, which is
+  considerably less accurate for badly scaled matrices and no faster. Note that `eig_vals` does not
+  accept the `expert` keyword argument.
+
 ### Removed
 
 ### Fixed
 
+- `schur_full` now supports `Diagonal` inputs through `DiagonalAlgorithm` ([#276](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/issues/276)).
 - LQ decompositions no longer gauge fix `Q` when `positive = false` and `L` is not computed.
 
 ### Performance

docs/src/user_interface/decompositions.md

@@ -113,11 +113,10 @@ The [Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition) tra
 It rewrites an arbitrary complex square matrix as unitarily similar to an upper triangular matrix whose diagonal elements are the eigenvalues of `A`.
 For real matrices, the same decomposition can be achieved in real arithmetic by allowing `T` to be quasi-upper triangular, i.e. triangular with blocks of size `(1, 1)` and `(2, 2)` on the diagonal.
 
-This decomposition is also useful for computing the eigenvalues of a matrix, which is exposed through the [`schur_vals`](@ref) function.
+The eigenvalues of `A` are returned alongside `T` and `Z`; use [`eig_vals`](@ref) to compute them on their own.
 
 ```@docs; canonical=false
 schur_full
-schur_vals
 ```
 
 The following algorithms are available for the Schur decomposition:

ext/MatrixAlgebraKitGenericSchurExt.jl

@@ -2,7 +2,7 @@ module MatrixAlgebraKitGenericSchurExt
 
 using MatrixAlgebraKit
 using MatrixAlgebraKit: check_input, GS, Driver
-import MatrixAlgebraKit: geev!, geevx!, gees!, eig_full!, eig_vals!, schur_full!, schur_vals!
+import MatrixAlgebraKit: geev!, geevx!, gees!, eig_full!, eig_vals!, schur_full!
 using LinearAlgebra: Diagonal, sorteig!
 using GenericSchur
 
@@ -49,9 +49,5 @@ Base.@deprecate(
     schur_full!(A, TZv, alg::GS_QRIteration),
     schur_full!(A, TZv, QRIteration(; driver = GS(), alg.kwargs...))
 )
-Base.@deprecate(
-    schur_vals!(A, vals, alg::GS_QRIteration),
-    schur_vals!(A, vals, QRIteration(; driver = GS(), alg.kwargs...))
-)
 
 end

src/implementations/schur.jl

@@ -1,10 +1,16 @@
 # Inputs
 # ------
 copy_input(::typeof(schur_full), A) = copy_input(eig_full, A)
-copy_input(::typeof(schur_vals), A) = copy_input(eig_vals, A)
 
 # check input
 function check_input(::typeof(schur_full!), A::AbstractMatrix, TZv, ::AbstractAlgorithm)
+    return _check_schur_full_input(A, TZv)
+end
+function check_input(::typeof(schur_full!), A::AbstractMatrix, TZv, ::DiagonalAlgorithm)
+    isdiag(A) || throw(DimensionMismatch("diagonal input matrix expected"))
+    return _check_schur_full_input(A, TZv)
+end
+function _check_schur_full_input(A::AbstractMatrix, TZv)
     m = LinearAlgebra.checksquare(A)
     T, Z, vals = TZv
     @assert T isa AbstractMatrix && Z isa AbstractMatrix && vals isa AbstractVector
@@ -16,13 +22,6 @@ function check_input(::typeof(schur_full!), A::AbstractMatrix, TZv, ::AbstractAl
     @check_scalar(vals, A, complex)
     return nothing
 end
-function check_input(::typeof(schur_vals!), A::AbstractMatrix, vals, ::AbstractAlgorithm)
-    m = LinearAlgebra.checksquare(A)
-    @assert vals isa AbstractVector
-    @check_size(vals, (m,))
-    @check_scalar(vals, A, complex)
-    return nothing
-end
 
 # Outputs
 # -------
@@ -32,21 +31,19 @@ function initialize_output(::typeof(schur_full!), A::AbstractMatrix, ::AbstractA
     vals = similar(A, complex(eltype(A)), n)
     return (A, Z, vals)
 end
-function initialize_output(::typeof(schur_vals!), A::AbstractMatrix, ::AbstractAlgorithm)
+# a diagonal matrix is already in Schur form, so `Z` can retain the diagonal structure
+function initialize_output(::typeof(schur_full!), A::Diagonal, ::DiagonalAlgorithm)
     n = size(A, 1) # square check will happen later
-    vals = similar(A, complex(eltype(A)), n)
-    return vals
+    return (A, similar(A), similar(A, complex(eltype(A)), n))
 end
 
 # DefaultAlgorithm intercepts
 # ---------------------------
-for f! in (:schur_full!, :schur_vals!)
-    @eval function $f!(A::AbstractMatrix, alg::DefaultAlgorithm)
-        return $f!(A, select_algorithm($f!, A, nothing; alg.kwargs...))
-    end
-    @eval function $f!(A::AbstractMatrix, out, alg::DefaultAlgorithm)
-        return $f!(A, out, select_algorithm($f!, A, nothing; alg.kwargs...))
-    end
+function schur_full!(A::AbstractMatrix, alg::DefaultAlgorithm)
+    return schur_full!(A, select_algorithm(schur_full!, A, nothing; alg.kwargs...))
+end
+function schur_full!(A::AbstractMatrix, out, alg::DefaultAlgorithm)
+    return schur_full!(A, out, select_algorithm(schur_full!, A, nothing; alg.kwargs...))
 end
 
 # ==========================
@@ -67,13 +64,9 @@ end
 # driver dispatch
 @inline schur_full_qr_iteration!(A, TZv; driver::Driver = DefaultDriver(), kwargs...) =
     schur_full_qr_iteration!(driver, A, TZv; kwargs...)
-@inline schur_vals_qr_iteration!(A, vals; driver::Driver = DefaultDriver(), kwargs...) =
-    schur_vals_qr_iteration!(driver, A, vals; kwargs...)
 
 @inline schur_full_qr_iteration!(::DefaultDriver, A, TZv; kwargs...) =
     schur_full_qr_iteration!(default_driver(QRIteration, A), A, TZv; kwargs...)
-@inline schur_vals_qr_iteration!(::DefaultDriver, A, vals; kwargs...) =
-    schur_vals_qr_iteration!(default_driver(QRIteration, A), A, vals; kwargs...)
 
 # Implementation
 function schur_full_qr_iteration!(driver::Driver, A, TZv; expert::Bool = false)
@@ -82,35 +75,34 @@ function schur_full_qr_iteration!(driver::Driver, A, TZv; expert::Bool = false)
     T === A || copy!(T, A)
     return TZv
 end
-function schur_vals_qr_iteration!(driver::Driver, A, vals; expert::Bool = false)
-    Z = similar(A, eltype(A), (size(A, 1), 0))
-    expert ? geesx!(driver, A, Z, vals) : gees!(driver, A, Z, vals)
-    return vals
-end
 
 # Top-level QRIteration dispatch
 function schur_full!(A::AbstractMatrix, TZv, alg::QRIteration)
     check_input(schur_full!, A, TZv, alg)
     schur_full_qr_iteration!(A, TZv; alg.kwargs...)
     return TZv
 end
-function schur_vals!(A::AbstractMatrix, vals, alg::QRIteration)
-    check_input(schur_vals!, A, vals, alg)
-    schur_vals_qr_iteration!(A, vals; alg.kwargs...)
-    return vals
+
+# Diagonal logic
+# --------------
+# `A` is already in Schur form, so `T = A` and `Z = I`, without any reordering
+function schur_full!(A::AbstractMatrix, TZv, alg::DiagonalAlgorithm)
+    check_input(schur_full!, A, TZv, alg)
+    T, Z, vals = TZv
+    if !has_equal_storage(A, T)
+        zero!(T)
+        diagview(T) .= diagview(A)
+    end
+    one!(Z)
+    vals .= diagview(A)
+    return TZv
 end
 
 # Deprecations
 # ------------
 for (lapack_algtype, expert_val) in ((:LAPACK_Simple, false), (:LAPACK_Expert, true))
-    @eval begin
-        Base.@deprecate(
-            schur_full!(A::AbstractMatrix, TZv, alg::$lapack_algtype),
-            schur_full!(A, TZv, QRIteration(; expert = $expert_val, alg.kwargs...))
-        )
-        Base.@deprecate(
-            schur_vals!(A::AbstractMatrix, vals, alg::$lapack_algtype),
-            schur_vals!(A, vals, QRIteration(; expert = $expert_val, alg.kwargs...))
-        )
-    end
+    @eval Base.@deprecate(
+        schur_full!(A::AbstractMatrix, TZv, alg::$lapack_algtype),
+        schur_full!(A, TZv, QRIteration(; expert = $expert_val, alg.kwargs...))
+    )
 end

src/interface/schur.jl

@@ -19,30 +19,17 @@ eigenvalues of `A`, as extracted from the (quasi-)diagonal of `T`.
 """
 @functiondef schur_full
 
-# TODO: is this useful? Is there any difference with simply `eig_vals`?
-"""
-    schur_vals(A; kwargs...) -> vals
-    schur_vals(A, alg::AbstractAlgorithm) -> vals
-    schur_vals!(A, [vals]; kwargs...) -> vals
-    schur_vals!(A, [vals], alg::AbstractAlgorithm) -> vals
-
-Compute the list of eigenvalues of `A` by computing the Schur decomposition of `A`.
-
-!!! note
-    The bang method `schur_vals!` optionally accepts the output structure and
-    possibly destroys the input matrix `A`. Always use the return value of the function
-    as it may not always be possible to use the provided `vals` as output.
-
-See also [`eig_full(!)`](@ref eig_full) and [`eig_trunc(!)`](@ref eig_trunc).
-"""
-@functiondef schur_vals
-
 # TODO: partial or truncated schur? Do we ever want or use this?
 
 # Algorithm selection
 # -------------------
-for f in (:schur_full!, :schur_vals!)
-    @eval function default_algorithm(::typeof($f), ::Type{A}; kwargs...) where {A}
-        return default_eig_algorithm(A; kwargs...)
-    end
+function default_algorithm(::typeof(schur_full!), ::Type{A}; kwargs...) where {A}
+    return default_eig_algorithm(A; kwargs...)
 end
+
+# Deprecations
+# ------------
+# `gees!` balances without scaling, so its eigenvalues are those of `eig_vals` with
+# `scale = false`, which is less accurate for badly scaled matrices and no faster
+Base.@deprecate schur_vals eig_vals false
+Base.@deprecate schur_vals! eig_vals! false

src/precompile.jl

@@ -24,7 +24,6 @@ using PrecompileTools: @compile_workload
         svd_vals(A)
 
         schur_full(A)
-        schur_vals(A)
 
         eigh_full(A)
         eigh_vals(A)

test/common/algorithms.jl

@@ -25,7 +25,7 @@ using MatrixAlgebraKit: LAPACK_SVDAlgorithm, PolarViaSVD, TruncatedAlgorithm,
     for f in (qr_full!, qr_full, qr_compact!, qr_compact, qr_null!, qr_null)
         @test @constinferred(default_algorithm(f, A)) == Householder()
     end
-    for f in (schur_full!, schur_full, schur_vals!, schur_vals)
+    for f in (schur_full!, schur_full)
         @test @constinferred(default_algorithm(f, A)) === QRIteration()
     end
 

test/decompositions/schur.jl

@@ -3,6 +3,7 @@ using Test
 using TestExtras
 using StableRNGs
 using LinearAlgebra: I, Diagonal
+using CUDA, AMDGPU
 
 if @isdefined(fast_tests) && fast_tests
     BLASFloats = (Float64, ComplexF64)
@@ -21,22 +22,31 @@ m = 54
 for T in (BLASFloats..., GenericFloats...)
     TestSuite.seed_rng!(123)
     if T ∈ BLASFloats
-        #=if CUDA.functional()
-            TestSuite.test_schur(CuMatrix{T}, (m, m); test_blocksize = false)
-            TestSuite.test_schur(Diagonal{T, CuVector{T}}, m; test_blocksize = false)
+        if CUDA.functional()
+            # dense GPU schur is not yet supported: there is no `gees!` for CUSOLVER
+            TestSuite.test_schur(Diagonal{T, CuVector{T}}, m)
+            TestSuite.test_schur_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),))
         end
+        #= not yet supported
         if AMDGPU.functional()
             TestSuite.test_schur(ROCMatrix{T}, (m, m); test_blocksize = false)
             TestSuite.test_schur(Diagonal{T, ROCVector{T}}, m; test_blocksize = false)
-        end=# # not yet supported
+        end=#
     end
     if !is_buildkite
         TestSuite.test_schur(T, (m, m))
         if T ∈ BLASFloats
             LAPACK_SCHUR_ALGS = (QRIteration(), QRIteration(expert = true))
             TestSuite.test_schur_algs(T, (m, m), LAPACK_SCHUR_ALGS)
         end
-        #AT = Diagonal{T, Vector{T}}
-        #TestSuite.test_schur(AT, m) # not supported yet
+        AT = Diagonal{T, Vector{T}}
+        TestSuite.test_schur(AT, m)
+        TestSuite.test_schur_algs(AT, m, (DiagonalAlgorithm(),))
     end
 end
+
+@testset "schur_vals is deprecated" begin
+    A = randn(StableRNG(123), 4, 4)
+    @test schur_vals(A) == eig_vals(A)
+    @test schur_vals!(copy(A)) == eig_vals!(copy(A))
+end

test/testsuite/decompositions/schur.jl

@@ -1,19 +1,20 @@
 using TestExtras
 using GenericSchur
 
+# `gees!` and `geev!` agree on well-scaled matrices, but not on the ordering for `Diagonal`
+sorted_vals(v) = sort!(collect(v); by = x -> (real(x), imag(x)))
+
 function test_schur(T::Type, sz; kwargs...)
     summary_str = testargs_summary(T, sz)
     return @testset "schur $summary_str" begin
         test_schur_full(T, sz; kwargs...)
-        test_schur_vals(T, sz; kwargs...)
     end
 end
 
 function test_schur_algs(T::Type, sz, algs; kwargs...)
     summary_str = testargs_summary(T, sz)
     return @testset "schur algorithms $summary_str" begin
         test_schur_full_algs(T, sz, algs; kwargs...)
-        test_schur_vals_algs(T, sz, algs; kwargs...)
     end
 end
 
@@ -26,13 +27,16 @@ function test_schur_full(
     return @testset "schur_full! $summary_str" begin
         A = instantiate_matrix(T, sz)
         Ac = deepcopy(A)
-        Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T))
+        Tc = complex(eltype(T))
 
         TA, Z, vals = @testinferred schur_full(A)
         @test eltype(TA) == eltype(Z) == eltype(T)
         @test eltype(vals) == Tc
         @test isisometric(Z)
         @test A * Z ≈ Z * TA
+        @test sorted_vals(vals) ≈ sorted_vals(eig_vals(A))
+        # a diagonal matrix is already in Schur form and is not reordered
+        A isa Diagonal && @test TA ≈ A
 
         TA2, Z2, vals2 = @testinferred schur_full!(Ac, (TA, Z, vals))
         @test TA2 === TA
@@ -51,13 +55,14 @@ function test_schur_full_algs(
     return @testset "schur_full! algorithm $alg $summary_str" for alg in algs
         A = instantiate_matrix(T, sz)
         Ac = deepcopy(A)
-        Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T))
+        Tc = complex(eltype(T))
 
         TA, Z, vals = @testinferred schur_full(A; alg)
         @test eltype(TA) == eltype(Z) == eltype(T)
         @test eltype(vals) == Tc
         @test isisometric(Z)
         @test A * Z ≈ Z * TA
+        @test sorted_vals(vals) ≈ sorted_vals(eig_vals(A))
 
         TA2, Z2, vals2 = @testinferred schur_full!(Ac, (TA, Z, vals); alg)
         @test TA2 === TA
@@ -66,47 +71,3 @@ function test_schur_full_algs(
         @test A * Z ≈ Z * TA
     end
 end
-
-function test_schur_vals(
-        T::Type, sz;
-        atol::Real = 0, rtol::Real = precision(T),
-        kwargs...
-    )
-    summary_str = testargs_summary(T, sz)
-    return @testset "schur_vals! $summary_str" begin
-        A = instantiate_matrix(T, sz)
-        Ac = deepcopy(A)
-        Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T))
-
-        valsc = @testinferred schur_vals(A)
-        @test eltype(valsc) == Tc
-        @test valsc ≈ eig_vals(A)
-
-        valsc = similar(A, Tc, size(A, 1))
-        valsc = @testinferred schur_vals!(Ac, valsc)
-        @test eltype(valsc) == Tc
-        @test valsc ≈ eig_vals(A)
-    end
-end
-
-function test_schur_vals_algs(
-        T::Type, sz, algs;
-        atol::Real = 0, rtol::Real = precision(T),
-        kwargs...
-    )
-    summary_str = testargs_summary(T, sz)
-    return @testset "schur_vals! algorithm $alg $summary_str" for alg in algs
-        A = instantiate_matrix(T, sz)
-        Ac = deepcopy(A)
-        Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T))
-
-        valsc = @testinferred schur_vals(A; alg)
-        @test eltype(valsc) == Tc
-        @test valsc ≈ eig_vals(A)
-
-        valsc = similar(A, Tc, size(A, 1))
-        valsc = @testinferred schur_vals!(Ac, valsc; alg)
-        @test eltype(valsc) == Tc
-        @test valsc ≈ eig_vals(A)
-    end
-end