use new generic zeros from MatrixPencils (#1004)

* use new generic zeros from MatrixPencils

* update tests

Author: baggepinnen

lib/ControlSystemsBase/Project.toml

@@ -38,7 +38,7 @@ ImplicitDifferentiation = "0.7.2"
 LinearAlgebra = "<0.0.1, 1"
 MacroTools = "0.5"
 MatrixEquations = "1, 2.1"
-MatrixPencils = "1.6"
+MatrixPencils = "1.8.3"
 Polynomials = "3.0, 4.0"
 PrecompileTools = "1"
 Printf = "<0.0.1, 1"

lib/ControlSystemsBase/src/analysis.jl

@@ -276,7 +276,7 @@ function tzeros(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::Abst
     tzeros(A2, B2, C2, D2)
 end
 
-function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), balance=true, kwargs...) where {T <: BlasFloat, E}
+function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), balance=true, kwargs...) where {T, E}
     if balance
         A, B, C = balance_statespace(A, B, C)
     end
@@ -289,110 +289,110 @@ function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}
     end
 end
 
-function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), kwargs...) where {T <: Union{AbstractFloat,Complex{<:AbstractFloat}}, E}
-    isempty(A) && return complex(T)[]
-
-    # Balance the system
-    A, B, C = balance_statespace(A, B, C; verbose=false)
-
-    # Compute a good tolerance
-    meps = 10*eps(real(T))*norm([A B; C D])
-
-    # Step 1:
-    A_r, B_r, C_r, D_r = reduce_sys(A, B, C, D, meps)
-
-    # Step 2: (conjugate transpose should be avoided since single complex zeros get conjugated)
-    A_rc, B_rc, C_rc, D_rc = reduce_sys(copy(transpose(A_r)), copy(transpose(C_r)), copy(transpose(B_r)), copy(transpose(D_r)), meps)
-
-    # Step 3:
-    # Compress cols of [C D] to [0 Df]
-    mat = [C_rc D_rc]
-    Wr = qr(mat').Q * I
-    W = reverse(Wr, dims=2)
-    mat = mat*W
-    if fastrank(mat', meps) > 0
-        nf = size(A_rc, 1)
-        m = size(D_rc, 2)
-        Af = ([A_rc B_rc] * W)[1:nf, 1:nf]
-        Bf = ([Matrix{T}(I, nf, nf) zeros(nf, m)] * W)[1:nf, 1:nf]
-        Z = schur(complex.(Af), complex.(Bf)) # Schur instead of eigvals to handle BigFloat
-        zs = Z.values
-        ControlSystemsBase._fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
-    else
-        zs = complex(T)[]
-    end
-    return zs
-end
-
-"""
-    reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
-Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
-A, B, C, D matrices. These are empty if there are no zeros.
-"""
-function reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
-    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
-    Cbar, Dbar = C, D
-    if isempty(A)
-        return A, B, C, D
-    end
-    while true
-        # Compress rows of D
-        U = qr(D).Q
-        D = U'D
-        C = U'C
-        sigma = fastrank(D, meps)
-        Cbar = C[1:sigma, :]
-        Dbar = D[1:sigma, :]
-        Ctilde = C[(1 + sigma):end, :]
-        if sigma == size(D, 1)
-            break
-        end
-
-        # Compress columns of Ctilde
-        V = reverse(qr(Ctilde').Q * I, dims=2)
-        Sj = Ctilde*V
-        rho = fastrank(Sj', meps)
-        nu = size(Sj, 2) - rho
-
-        if rho == 0
-            break
-        elseif nu == 0
-            # System has no zeros, return empty matrices
-            A = B = Cbar = Dbar = Matrix{T}(undef, 0,0)
-            break
-        end
-        # Update System
-        n, m = size(B)
-        Vm = [V zeros(T, n, m); zeros(T, m, n) Matrix{T}(I, m, m)] # I(m) is not used for type stability reasons (as of julia v1.7)
-        if sigma > 0
-            M = [A B; Cbar Dbar]
-            Vs = [copy(V') zeros(T, n, sigma) ; zeros(T, sigma, n) Matrix{T}(I, sigma, sigma)]
-        else
-            M = [A B]
-            Vs = copy(V')
-        end
-        sigma, rho, nu
-        M = Vs * M * Vm
-        A = M[1:nu, 1:nu]
-        B = M[1:nu, (nu + rho + 1):end]
-        C = M[(nu + 1):end, 1:nu]
-        D = M[(nu + 1):end,  (nu + rho + 1):end]
-    end
-    return A, B, Cbar, Dbar
-end
-
-# Determine the number of non-zero rows, with meps as a tolerance. For an
-# upper-triangular matrix, this is a good proxy for determining the row-rank.
-function fastrank(A::AbstractMatrix, meps::Real)
-    n, m = size(A)
-    if n*m == 0     return 0    end
-    norms = Vector{real(eltype(A))}(undef, n)
-    for i = 1:n
-        norms[i] = norm(A[i, :])
-    end
-    mrank = sum(norms .> meps)
-    return mrank
-end
+# function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), kwargs...) where {T <: Union{AbstractFloat,Complex{<:AbstractFloat}}, E}
+#     isempty(A) && return complex(T)[]
+
+#     # Balance the system
+#     A, B, C = balance_statespace(A, B, C; verbose=false)
+
+#     # Compute a good tolerance
+#     meps = 10*eps(real(T))*norm([A B; C D])
+
+#     # Step 1:
+#     A_r, B_r, C_r, D_r = reduce_sys(A, B, C, D, meps)
+
+#     # Step 2: (conjugate transpose should be avoided since single complex zeros get conjugated)
+#     A_rc, B_rc, C_rc, D_rc = reduce_sys(copy(transpose(A_r)), copy(transpose(C_r)), copy(transpose(B_r)), copy(transpose(D_r)), meps)
+
+#     # Step 3:
+#     # Compress cols of [C D] to [0 Df]
+#     mat = [C_rc D_rc]
+#     Wr = qr(mat').Q * I
+#     W = reverse(Wr, dims=2)
+#     mat = mat*W
+#     if fastrank(mat', meps) > 0
+#         nf = size(A_rc, 1)
+#         m = size(D_rc, 2)
+#         Af = ([A_rc B_rc] * W)[1:nf, 1:nf]
+#         Bf = ([Matrix{T}(I, nf, nf) zeros(nf, m)] * W)[1:nf, 1:nf]
+#         Z = schur(complex.(Af), complex.(Bf)) # Schur instead of eigvals to handle BigFloat
+#         zs = Z.values
+#         ControlSystemsBase._fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
+#     else
+#         zs = complex(T)[]
+#     end
+#     return zs
+# end
+
+# """
+#     reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
+# Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
+# A, B, C, D matrices. These are empty if there are no zeros.
+# """
+# function reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
+#     T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
+#     Cbar, Dbar = C, D
+#     if isempty(A)
+#         return A, B, C, D
+#     end
+#     while true
+#         # Compress rows of D
+#         U = qr(D).Q
+#         D = U'D
+#         C = U'C
+#         sigma = fastrank(D, meps)
+#         Cbar = C[1:sigma, :]
+#         Dbar = D[1:sigma, :]
+#         Ctilde = C[(1 + sigma):end, :]
+#         if sigma == size(D, 1)
+#             break
+#         end
+
+#         # Compress columns of Ctilde
+#         V = reverse(qr(Ctilde').Q * I, dims=2)
+#         Sj = Ctilde*V
+#         rho = fastrank(Sj', meps)
+#         nu = size(Sj, 2) - rho
+
+#         if rho == 0
+#             break
+#         elseif nu == 0
+#             # System has no zeros, return empty matrices
+#             A = B = Cbar = Dbar = Matrix{T}(undef, 0,0)
+#             break
+#         end
+#         # Update System
+#         n, m = size(B)
+#         Vm = [V zeros(T, n, m); zeros(T, m, n) Matrix{T}(I, m, m)] # I(m) is not used for type stability reasons (as of julia v1.7)
+#         if sigma > 0
+#             M = [A B; Cbar Dbar]
+#             Vs = [copy(V') zeros(T, n, sigma) ; zeros(T, sigma, n) Matrix{T}(I, sigma, sigma)]
+#         else
+#             M = [A B]
+#             Vs = copy(V')
+#         end
+#         sigma, rho, nu
+#         M = Vs * M * Vm
+#         A = M[1:nu, 1:nu]
+#         B = M[1:nu, (nu + rho + 1):end]
+#         C = M[(nu + 1):end, 1:nu]
+#         D = M[(nu + 1):end,  (nu + rho + 1):end]
+#     end
+#     return A, B, Cbar, Dbar
+# end
+
+# # Determine the number of non-zero rows, with meps as a tolerance. For an
+# # upper-triangular matrix, this is a good proxy for determining the row-rank.
+# function fastrank(A::AbstractMatrix, meps::Real)
+#     n, m = size(A)
+#     if n*m == 0     return 0    end
+#     norms = Vector{real(eltype(A))}(undef, n)
+#     for i = 1:n
+#         norms[i] = norm(A[i, :])
+#     end
+#     mrank = sum(norms .> meps)
+#     return mrank
+# end
 
 """
     relative_gain_array(G, w::AbstractVector)

lib/ControlSystemsBase/test/test_analysis.jl

@@ -1,6 +1,6 @@
 @test_throws MethodError poles(big(1.0)*ssrand(1,1,1)) # This errors before loading GenericSchur
 using GenericSchur # Required to compute eigvals (in tzeros and poles) of a matrix with exotic element types
-@testset "test_analysis" begin
+
 es(x) = sort(x, by=LinearAlgebra.eigsortby)
 ## tzeros ##
 # Examples from the Emami-Naeini & Van Dooren Paper
@@ -23,7 +23,7 @@ ex_3 = ss(A, B, C, D)
 @test es(tzeros(ex_3)) ≈ es([0.3411639019140099 + 1.161541399997252im,
                              0.3411639019140099 - 1.161541399997252im,
                              0.9999999999999999 + 0.0im,
-                             -0.6823278038280199 + 0.0im])
+                             -0.6823278038280199 + 0.0im]) atol=1e-6
 # Example 4
 A = [-0.129    0.0   0.396e-1  0.25e-1    0.191e-1;
      0.329e-2  0.0  -0.779e-4  0.122e-3  -0.621;
@@ -58,7 +58,7 @@ C = [0 -1 0]
 D = [0]
 ex_6 = ss(A, B, C, D)
 @test tzeros(ex_6) == [2] # From paper: "Our algorithm will extract the singular part of S(A) and will yield a regular pencil containing the single zero at 2."
-@test_broken tzeros(big(1.0)ex_6) == [2]
+@test tzeros(big(1.0)ex_6) == [2]
 @test ControlSystemsBase.count_integrators(ex_6) == 2
 
 @test ss(A, [0 0 1]', C, D) == ex_6
@@ -79,8 +79,8 @@ C = [0 0 0 1 0 0]
 D = [0]
 ex_8 = ss(A, B, C, D)
 # TODO : there may be a way to improve the precision of this example.
-@test tzeros(ex_8) ≈ [-1.0, -1.0] atol=1e-7
-@test tzeros(big(1)ex_8) ≈ [-1.0, -1.0] atol=1e-7
+@test tzeros(ex_8) ≈ [-1.0, -1.0] atol=1e-6
+@test tzeros(big(1)ex_8) ≈ [-1.0, -1.0] atol=1e-12
 @test ControlSystemsBase.count_integrators(ex_8) == 0
 
 # Example 9
@@ -103,7 +103,7 @@ D = [0 0;
      0 0;
      0 0]
 ex_11 = ss(A, B, C, D)
-@test tzeros(ex_11) ≈ [4.0, -3.0]
+@test tzeros(ex_11) ≈ [4.0, -3.0] atol=1e-5
 @test tzeros(big(1)ex_11) ≈ [4.0, -3.0]
 
 # Test for multiple zeros, siso tf
@@ -405,7 +405,6 @@ end
 @test length(tzeros(G)) == 3
 @test es(tzeros(G)) ≈ es(tzeros(big(1)G))
 
-end
 
 
 ## large TF poles and zeros