diff --git a/Project.toml b/Project.toml index 6e5eb80..f9f69bd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MatrixPencils" uuid = "48965c70-4690-11ea-1f13-43a2532b2fa8" authors = ["Andreas Varga "] -version = "0.3.0" +version = "0.4.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index f88ff46..2a23358 100644 --- a/README.md +++ b/README.md @@ -26,6 +26,7 @@ The current version of the package includes the following functions: **Manipulation of structured linear matrix pencils of the form [A-λE B; C D]** +* **sreduceBF** Reduction to the basic condensed form [B A-λE; D C] with E upper triangular and nonsingular. * **sklf** Computation of the Kronecker-like form exhibiting the full Kronecker structure * **sklf_left** Computation of the Kronecker-like form exhibiting the left Kronecker structure * **sklf_right** Computation of the Kronecker-like form exhibiting the right Kronecker structure @@ -44,13 +45,19 @@ The current version of the package includes the following functions: **Some applications to structured linear matrix pencils of the form [A-λE B; C D]** -* ***spkstruct** Determination of the complete Kronecker structure +* **spkstruct** Determination of the complete Kronecker structure * **sprank** Determination of the normal rank * **speigvals** Computation of the finite and infinite eigenvalues * **spzeros** Computation of the finite and infinite (Smith) zeros +**Manipulation of linearizations of the form [A-λE B; C D] of polynomial or rational matrices** + +* **lsminreal** Computation of least order strong liniarizations of rational matrices. +* **lsminreal2** Computation of least order strong liniarizations of rational matrices (potentially more efficient). +* **lsequal** Check the equivalence two linearizations. + A complete list of implemented functions is available [here](https://sites.google.com/site/andreasvargacontact/home/software/matrix-pencils-in-julia). ## Future plans -The collection of tools will be extended by adding new functionality, such as tools for the manipulation of regular pencils (e.g., reduction to a block-diagonal structure, eigenvalue assignment), manipulation of structured linear pencils arising from the linearization of polynomial and rational matrices (e.g., computation of strong linearizations), etc. +The collection of tools will be extended by adding new functionality, such as tools for the manipulation of regular pencils (e.g., reduction to a block-diagonal structure, eigenvalue assignment), building linearizations of polynomial matrices, applications of structured linear pencils manipulations to polynomial matrix problems, etc. diff --git a/docs/make.jl b/docs/make.jl index 4739bbc..9e0518e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,7 +13,8 @@ makedocs( "sklftools.md", "klfapps.md", "sklfapps.md", - "pregular.md" + "pregular.md", + "lstools.md" ], "Internal" => [ "klftools_int.md" diff --git a/docs/src/index.md b/docs/src/index.md index 30acc50..1542d43 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -39,6 +39,7 @@ The current version of the package includes the following functions: | Function | Description | | :--- | :--- | +| **sreduceBF** | Reduction to the basic condensed form [B A-λE; D C] with E upper triangular and nonsingular. | | **sklf** | Computation of the Kronecker-like form exhibiting the full Kronecker structure | | **sklf_left** | Computation of the Kronecker-like form exhibiting the left Kronecker structure | | **sklf_right** | Computation of the Kronecker-like form exhibiting the right Kronecker structure | @@ -68,14 +69,26 @@ The current version of the package includes the following functions: | **speigvals** | Computation of the finite and infinite eigenvalues | | **spzeros** | Computation of the finite and infinite (Smith) zeros | +**Manipulation of linearizations of the form [A-λE B; C D] of polynomial or rational matrices** + +| Function | Description | +| :--- | :--- | +| **lsminreal** | Computation of least order strong liniarizations of rational matrices | +| **lsminreal2** | Computation of least order strong liniarizations of rational matrices (potentially more efficient)| +| **lsequal** | Check the equivalence two linearizations | + A complete list of implemented functions is available [here](https://sites.google.com/site/andreasvargacontact/home/software/matrix-pencils-in-julia). ## Future plans -The collection of tools will be extended by adding new functionality, such as tools for the manipulation of regular pencils (e.g., reduction to a block-diagonal structure, eigenvalue assignment), manipulation of structured linear pencils arising from the linearization of polynomial and rational matrices (e.g., computation of strong linearizations), etc. +The collection of tools will be extended by adding new functionality, such as tools for the manipulation of regular pencils (e.g., reduction to a block-diagonal structure, eigenvalue assignment), building linearizations of polynomial matrices, applications of structured linear pencils manipulations to polynomial matrix problems, etc. ## Release Notes +### Version 0.4.0 + +This release includes implementations of computational procedures to determine least order linerizations, as those which may arise from the linearization of polynomial and rational matrices. The elimination of simple eigenvalues is based on a new function to compute SVD-like forms of regular matrix pencils. A new function to check the equivalence of two linearizations is also provided. + ### Version 0.3.0 This release includes implementations of several reductions of structured linear pencils to various KLFs and covers some straightforward applications such as the computation of finite and infinite eigenvalues and zeros, the determination of the normal rank, the determination of Kronecker indices and infinite eigenvalue structure. diff --git a/docs/src/klftools_int.md b/docs/src/klftools_int.md index 18ae28f..f4c21fb 100644 --- a/docs/src/klftools_int.md +++ b/docs/src/klftools_int.md @@ -8,6 +8,8 @@ MatrixPencils._preduce2! MatrixPencils._preduce4! MatrixPencils._sreduceB! MatrixPencils._sreduceBA! +MatrixPencils._sreduceBAE! MatrixPencils._sreduceC! MatrixPencils._sreduceAC! +MatrixPencils._sreduceAEC! ``` diff --git a/docs/src/lstools.md b/docs/src/lstools.md new file mode 100644 index 0000000..34012d4 --- /dev/null +++ b/docs/src/lstools.md @@ -0,0 +1,11 @@ +# Manipulation of linearizations of the form [A-λE B; C D] of polynomial or rational matrices + +```@docs +lsequal +lsminreal +lsminreal2 +sklf_right! +sklf_left! +sklf_rightfin! +sklf_leftfin! +``` diff --git a/docs/src/makeindex.md b/docs/src/makeindex.md index caeda95..b6d66b3 100644 --- a/docs/src/makeindex.md +++ b/docs/src/makeindex.md @@ -1,7 +1,7 @@ # Index ```@index -Pages = ["klftools.md", "klfapps.md", "sklftools.md", "sklfapps.md", "pregular","klftools_int.md" ] +Pages = ["klftools.md", "klfapps.md", "sklftools.md", "sklfapps.md", "pregular","lstools.md","klftools_int.md" ] Module = ["MatrixPencils"] Order = [:type, :function] ``` diff --git a/docs/src/pregular.md b/docs/src/pregular.md index 9fa854b..5a0c70b 100644 --- a/docs/src/pregular.md +++ b/docs/src/pregular.md @@ -3,6 +3,7 @@ ```@docs isregular fisplit -sklf_right! -sklf_left! +MatrixPencils.fisplit! +MatrixPencils._qrE! +MatrixPencils._svdlikeAE! ``` diff --git a/src/MatrixPencils.jl b/src/MatrixPencils.jl index e9ff1e4..c56e5a5 100644 --- a/src/MatrixPencils.jl +++ b/src/MatrixPencils.jl @@ -14,9 +14,12 @@ using .LapackUtil2: larfg!, larfgl!, larf! export klf, klf_left, klf_right, klf_rlsplit export prank, pkstruct, peigvals, pzeros, KRInfo -export isregular, fisplit -export sreduceBF, sklf, sklf_right, sklf_left, sklf_right!, sklf_left! +export isregular, fisplit, _svdlikeAE! +export sreduceBF, sklf, sklf_right, sklf_left, sklf_right!, sklf_left!, sklf_rightfin!, sklf_leftfin! export sprank, spkstruct, speigvals, spzeros +export lsminreal, lsminreal2, lsequal +export lsbuild, lpbuild, lpbuildCF1, lpbuildCF2 +export polbascoeffs, polbasparams, polbas2mon, polmon2bas, poldeg, polreverse #export dss, ss, gnrank, gzero, gpole #export LTISystem, AbstractDescriptorStateSpace #export _preduceBF!, _preduce1!, _preduce2!, _preduce3!, _preduce4! @@ -28,6 +31,8 @@ include("regtools.jl") include("klfapps.jl") include("sklftools.jl") include("sklfapps.jl") +include("lstools.jl") +#include("poltools.jl") #include("dstools.jl") include("lputil.jl") include("slputil.jl") diff --git a/src/lputil.jl b/src/lputil.jl index 3b8ea2d..629faa8 100644 --- a/src/lputil.jl +++ b/src/lputil.jl @@ -37,9 +37,9 @@ The rank decision based on the SVD-decomposition is generally more reliable, but function _preduceBF!(M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, L::Union{AbstractMatrix{T1},Missing} = missing, R::Union{AbstractMatrix{T1},Missing} = missing; - atol::Real = zero(eltype(M)), rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(atol), - fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, - withQ = true, withZ = true) where T1 <: BlasFloat + atol::Real = zero(real(T1)), rtol::Real = (min(size(M)...)*eps(real(float(one(T1)))))*iszero(atol), + fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat # In interest of performance, no dimensional checks are performed mM, nM = size(M) T = eltype(M) @@ -64,7 +64,7 @@ function _preduceBF!(M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, i12 = 1:roff+npp j23 = coff+1:nM jN23 = coff+npm+1:nM - eltype(M) <: Complex ? tran = 'C' : tran = 'T' + T <: Complex ? tran = 'C' : tran = 'T' if fast # compute in-place the QR-decomposition N22*P2 = Q2*[R2;0] with column pivoting N22 = view(N,i22,j22) @@ -109,7 +109,7 @@ function _preduceBF!(M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, end # R <- R*Z2*Pc2 ismissing(R) || (LinearAlgebra.LAPACK.ormrz!('R',tran,N22r,tau,view(R,:,j22)); R[:,j22] = R[:,jp]) - N22[:,:] = [ zeros(n,m) triu(N22[i1,i1]); zeros(p,npm) ] + N22[:,:] = [ zeros(T,n,m) triu(N22[i1,i1]); zeros(T,p,npm) ] else # compute the complete orthogonal decomposition of N22 using the SVD-decomposition U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',view(N,i22,j22)) @@ -128,9 +128,9 @@ function _preduceBF!(M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, N[i11,j22] = N[i11,jp] # N23 <- U'*N23 N[i22,jN23] = U'*N[i22,jN23] - # Q <- Q*SVD.U + # Q <- Q*U withQ && (Q[:,i22] = Q[:,i22]*U) - # Z <- Q*SVD.V + # Z <- Q*V if withZ Z[:,j22] = Z[:,j22]*Vt' Z[:,j22] = Z[:,jp] @@ -139,7 +139,7 @@ function _preduceBF!(M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, ismissing(L) || (L[i22,:] = U'*L[i22,:]) # R <- R*V ismissing(R) || (R[:,j22] = R[:,j22]*Vt'; R[:,j22] = R[:,jp] ) - N[i22,j22] = [ zeros(n,m) Diagonal(S[1:n]) ; zeros(p,npm) ] + N[i22,j22] = [ zeros(T,n,m) Diagonal(S[1:n]) ; zeros(T,p,npm) ] end return n, m, p end @@ -173,9 +173,10 @@ The performed orthogonal or unitary transformations are accumulated in `Q` (i.e. The matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the matrix `R` is overwritten by `R*Z1` unless `R = missing`. """ function _preduce1!(n::Int, m::Int, p::Int, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, - Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real, L::Union{AbstractMatrix{T1},Missing} = missing, R::Union{AbstractMatrix{T1},Missing} = missing; - fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) where T1 <: BlasFloat + fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat # Steps 1 and 2: QR- or SVD-decomposition based full column compression of [B; D] with the UT-form preservation of E # Reduce the structured pencil # @@ -281,9 +282,10 @@ The performed orthogonal or unitary transformations are accumulated in `Q` (i.e. The matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the matrix `R` is overwritten by `R*Z1` unless `R = missing`. """ function _preduce2!(n::Int, m::Int, p::Int, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, - Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real, L::Union{AbstractMatrix{T1},Missing} = missing, R::Union{AbstractMatrix{T1},Missing} = missing; - fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) where T1 <: BlasFloat + fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat # Steps 1 and 2: QR- or SVD-decomposition based full column compression of [D C] with the UT-form preservation of E # Reduce the structured pencil # @@ -403,9 +405,10 @@ The performed orthogonal or unitary transformations are accumulated in `Q` (i.e. The matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the matrix `R` is overwritten by `R*Z1` unless `R = missing`. """ function _preduce3!(n::Int, m::Int, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, - Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real, L::Union{AbstractMatrix{T1},Missing} = missing, R::Union{AbstractMatrix{T1},Missing} = missing; - fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) where T1 <: BlasFloat + fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat # Step 3: QR- or SVD-decomposition based full row compression of B and UT-form preservation of E # Reduce the structured pencil # @@ -558,9 +561,10 @@ The performed orthogonal or unitary transformations are accumulated in `Q` (i.e. The matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the matrix `R` is overwritten by `R*Z1` unless `R = missing`. """ function _preduce4!(n::Int, m::Int, p::Int, M::AbstractMatrix{T1},N::AbstractMatrix{T1}, - Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real, L::Union{AbstractMatrix{T1},Missing} = missing, R::Union{AbstractMatrix{T1},Missing} = missing; - fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) where T1 <: BlasFloat + fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat # Step 4: QR- or SVD-decomposition based full column compression of C and UT-form preservation of E # Reduce the structured pencil # diff --git a/src/lstools.jl b/src/lstools.jl new file mode 100644 index 0000000..aa07cee --- /dev/null +++ b/src/lstools.jl @@ -0,0 +1,1121 @@ +""" + lsequal(A1, E1, B1, C1, D1, A2, E2, B2, C2, D2; fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ) -> flag::Bool + +Check if two linearizations `(A1-λE1,B1,C1,D1)` and `(A2-λE2,B2,C2,D2)` satisfy the equivalence condition + + -1 -1 + C1*(A1-λE1) *B1 + D1 = C2*(A2-λE2) *B2 + D2 + +The ckeck is performed by computing the normal rank `n` of the structured linear matrix pencil `M - λN` + + | A1-λE1 0 | B1 | + | 0 A2-λE2 | B2 | + M - λN = |---------------|-------| + | C1 -C2 | D1-D2 | + +and verifying that `n = n1+n2`, where `n1` and `n2` are the orders of the square matrices `A1` and `A2`, respectively. + +If `fastrank = true`, the rank is evaluated by counting how many singular values of `M - γ N` have magnitude +greater than `max(max(atol1,atol2), rtol*σ₁)`, where `σ₁` is the largest singular value of `M - γ N` and +`γ` is a randomly generated value [1]. +If `fastrank = false`, the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums +of right and left Kronecker indices, respectively, while `ni` and `nf` are the number of finite and +infinite eigenvalues, respectively. The sums `nr+ni` and `nf+nl`, are determined from an +appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures +of the pencil `M - λN`. For efficiency purpose, the reduction to the relevant KLF is only partially performed +using rank decisions based on rank revealing SVD-decompositions. + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance +for the nonzero elements of `M` and `N`. The default relative tolerance is `k*ϵ`, where `k` is the size of +the smallest dimension of `M`, and `ϵ` is the machine epsilon of the element type of `M`. + +[1] A. Varga, On Checking Null Rank Conditions of Rational Matrices, https://arxiv.org/abs/1812.11396, 2018. +""" +function lsequal(A1::AbstractMatrix, E1::Union{AbstractMatrix,UniformScaling{Bool}}, + B1::AbstractMatrix, C1::AbstractMatrix, D1::AbstractMatrix, + A2::AbstractMatrix, E2::Union{AbstractMatrix,UniformScaling{Bool}}, + B2::AbstractMatrix, C2::AbstractMatrix, D2::AbstractMatrix; + atol1::Real = zero(real(eltype(A1))), atol2::Real = zero(real(eltype(A1))), + rtol::Real = min(size(A1)...)*eps(real(float(one(real(eltype(A1))))))*iszero(max(atol1,atol2)), + fastrank::Bool = true) + + # implicit dimensional checks are performed using the try-catch scheme + try + T = promote_type(eltype(A1), eltype(A2)) + n1 = size(A1,1) + n2 = size(A2,1) + A = [A1 zeros(T,n1,n2); + zeros(T,n2,n1) A2] + B = [B1; B2;] + C = [C1 -C2;] + D = [D1-D2;] + if E1 == I && E2 == I + E = I + else + E = [E1 zeros(T,n1,n2); + zeros(T,n2,n1) E2] + end + return (sprank(A,E,B,C,D,atol1 = atol1, atol2 = atol2, rtol = rtol, fastrank = fastrank) == n1+n2) + catch + return false + end +end +""" + lsminreal2(A, E, B, C, D; fast = true, atol1 = 0, atol2 = 0, rtol, finite = true, infinite = true, contr = true, obs = true, noseig = true) + -> Ar, Er, Br, Cr, Dr, nuc, nuo, nse + +Reduce the linearization `(A-λE,B,C,D)` of a rational matrix to a reduced form `(Ar-λEr,Br,Cr,Dr)` such that + + -1 -1 + C*(A-λE) *B + D = Cr*(Ar-λEr) *Br + Dr + +with the least possible order `nr` of `Ar-λEr` if `finite = true`, `infinite = true`, +`contr = true`, `obs = true` and `nseig = false`. +The reduced order linearization satisfies: + + (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability) + + (2) rank[Br Er] = nr (infinite controllability) + + (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability) + + (4) rank[Er; Cr] = nr (infinite observability) + + (5) Ar-λEr has no simple eigenvalues + +The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and +respectively, condition (5) are returned in `nuc`, `nuo`, `nse`. + +Some reduction steps can be skipped by appropriately selecting the keyword arguments +`contr`, `obs`, `finite`, `infinite` and `nseig`. + +If `contr = false`, then the controllability conditions (1) and (2) are not enforced. +If `contr = true` and `finite = true`, then the finite controllability condition (1) is enforced. +If `contr = true` and `infinite = true`, then the infinite controllability condition (2) is enforced. + +If `obs = false`, then observability condition (3) and (4) are not enforced. +If `obs = true` and `finite = true`, then the finite observability condition (3) is enforced. +If `obs = true` and `infinite = true`, then the infinite observability condition (4) is enforced. + +If `nseig = false`, then condition (5) on the lack of simple eigenvalues is not enforced. + +To enforce conditions (1)-(4), the `Procedure GIR` in `[1, page 328]` is employed, which performs +orthogonal similarity transformations on the matrices of the original linearization `(A-λE,B,C,D)` +to obtain an irreducible linearization using structured pencil reduction algorithms. +To enforce condition (5), residualization formulas (see, e.g., `[1, page 329]`) are employed which +involves matrix inversions. + +The underlying pencil manipulation algorithms employ rank determinations based on either the use of +rank revealing QR-decomposition with column pivoting, if `fast = true`, or in the SVD-decomposition. +The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher. + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`, +and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`. + +[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017. +""" +function lsminreal2(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}}, + B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix; + atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), + rtol::Real = min(size(A)...)*eps(real(float(one(real(eltype(A))))))*iszero(max(atol1,atol2)), + fast::Bool = true, finite::Bool = true, infinite::Bool = true, + contr::Bool = true, obs::Bool = true, noseig::Bool = true) + + eident = (typeof(E) == UniformScaling{Bool}) || isequal(E,I) + emat = (typeof(E) <: AbstractMatrix) + isa(A,Adjoint) && (A = copy(A)) + isa(E,Adjoint) && (E = copy(E)) + n = LinearAlgebra.checksquare(A) + emat && (n,n) !== size(E) && throw(DimensionMismatch("A and E must have the same dimensions")) + p, m = size(D) + (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions")) + (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions")) + T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D)) + eident || (T = promote_type(T,eltype(E))) + T <: BlasFloat || (T = promote_type(Float64,T)) + + eltype(A) !== T && (A = convert(Matrix{T},A)) + (!eident && eltype(E) !== T) && (E = convert(Matrix{T},E)) + eltype(B) !== T && (B = convert(Matrix{T},B)) + eltype(C) !== T && (C = convert(Matrix{T},C)) + eltype(D) !== T && (D = convert(Matrix{T},D)) + + + (n == 0 || m == 0 || p == 0) && (return A, E, B, C, D, 0, 0, 0) + # save system matrices + Ar = copy(A) + Br = copy(B) + Cr = copy(C) + Dr = copy(D) + Er = copy(E) + ir = 1:n + if eident + if contr + _, _, nr, nuc = sklf_right!(Ar, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol1, rtol = rtol, withQ = false) + if nuc > 0 + ir = 1:nr + # save intermediary results + A = Ar[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore original matrices + Ar = copy(A) + Br = copy(B) + Cr = copy(C) + end + else + nuc = 0 + end + if obs + _, _, no, nuo = sklf_left!(view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:); fast = fast, atol1 = atol1, atol2 = atol1, rtol = rtol, withQ = false) + if nuo > 0 + ir = ir[end-no+1:end] + else + # restore saved matrices + Ar[ir,ir] = A + Br[ir,:] = B + Cr[:,ir] = C + end + else + nuo = 0 + end + return Ar[ir,ir], emat ? Er[ir,ir] : I, Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0 + else + if finite + if contr + _, _, _, nr, nfuc = sklf_rightfin!(Ar, Er, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false) + if nfuc > 0 + ir = 1:nr + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore original matrices + Ar = copy(A) + Er = copy(E) + Br = copy(B) + Cr = copy(C) + end + else + nfuc = 0 + end + if obs + _, _, _, no, nfuo = sklf_leftfin!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Br,ir,:); + fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false) + if nfuo > 0 + ir = ir[end-no+1:end] + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore saved matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + else + nfuo = 0 + end + else + nfuc = 0 + nfuo = 0 + end + if infinite + if contr + _, _, _, nr, niuc = sklf_rightfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Br,ir,:), view(Cr,:,ir); + fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, + withQ = false, withZ = false) + if niuc > 0 + ir = ir[1:nr] + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore original matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + else + niuc = 0 + end + if obs + _, _, _, no, niuo = sklf_leftfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:); + fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, + withQ = false, withZ = false) + if niuo > 0 + ir = ir[end-no+1:end] + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore saved matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + else + niuo = 0 + end + else + niuc = 0 + niuo = 0 + end + nuc = nfuc+niuc + nuo = nfuo+niuo + if noseig + rE, rA22 = _svdlikeAE!(view(Ar,ir,ir), view(Er,ir,ir), nothing, nothing, view(Br,ir,:), view(Cr,:,ir), + fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false) + if rA22 > 0 + i1 = ir[1:rE] + i2 = ir[rE+1:rE+rA22] + # make A22 = I + fast ? (A22 = UpperTriangular(Ar[i2,i2])) : (A22 = Diagonal(Ar[i2,i2])) + ldiv!(A22,view(Ar,i2,i1)) + ldiv!(A22,view(Br,i2,:)) + # apply simplified residualization formulas + Dr -= Cr[:,i2]*Br[i2,:] + Br[i1,:] -= Ar[i1,i2]*Br[i2,:] + Cr[:,i1] -= Cr[:,i2]*Ar[i2,i1] + Ar[i1,i1] -= Ar[i1,i2]*Ar[i2,i1] + ir = [i1; ir[rE+rA22+1:end]] + else + # restore saved matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, rA22 + else + return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0 + end + end +end +""" + lsminreal(A, E, B, C, D; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true, noseig = true) -> Ar, Er, Br, Cr, Dr, nuc, nuo, nse + +Reduce the linearization `(A-λE,B,C,D)` of a rational matrix to a reduced form `(Ar-λEr,Br,Cr,Dr)` such that + + -1 -1 + C*(A-λE) *B + D = Cr*(Ar-λEr) *Br + Dr + +with the least possible order `nr` of `Ar-λEr` if `contr = true`, `obs = true` and `nseig = false`. +The reduced order linearization satisfies: + + (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability) + + (2) rank[Br Er] = nr (infinite controllability) + + (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability) + + (4) rank[Er; Cr] = nr (infinite observability) + + (5) Ar-λEr has no simple eigenvalues + +The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and +respectively, condition (5) are returned in `nuc`, `nuo`, `nse`. + +Some reduction steps can be skipped by appropriately selecting the keyword arguments +`contr`, `obs` and `nseig`. + +If `contr = false`, then the controllability conditions (1) and (2) are not enforced. + +If `obs = false`, then observability condition (3) and (4) are not enforced. + +If `nseig = false`, then condition (5) on the lack of simple eigenvalues is not enforced. + +To enforce conditions (1)-(4), orthogonal similarity transformations are performed on +the matrices of the original linearization `(A-λE,B,C,D)` to obtain an irreducible linearization using +structured pencil reduction algorithms, as the fast versions of the reduction techniques of the +full row rank pencil [B A-λE] and full column rank pencil [A-λE;C] proposed in [1]. +To enforce condition (5), residualization formulas (see, e.g., `[2, page 329]`) are employed which +involves matrix inversions. + +The underlying pencil manipulation algorithms employ rank determinations based on either the use of +rank revealing QR-decomposition with column pivoting, if `fast = true`, or in the SVD-decomposition. +The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher. + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`, +and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`. + +[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory, +IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981. + +[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017. +""" +function lsminreal(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}}, + B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix; + atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), + rtol::Real = min(size(A)...)*eps(real(float(one(real(eltype(A))))))*iszero(max(atol1,atol2)), + fast::Bool = true, contr::Bool = true, obs::Bool = true, noseig::Bool = true) + + eident = (typeof(E) == UniformScaling{Bool}) || isequal(E,I) + emat = (typeof(E) <: AbstractMatrix) + isa(A,Adjoint) && (A = copy(A)) + isa(E,Adjoint) && (E = copy(E)) + n = LinearAlgebra.checksquare(A) + emat && (n,n) !== size(E) && throw(DimensionMismatch("A and E must have the same dimensions")) + p, m = size(D) + (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions")) + (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions")) + T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D)) + eident || (T = promote_type(T,eltype(E))) + T <: BlasFloat || (T = promote_type(Float64,T)) + + eltype(A) !== T && (A = convert(Matrix{T},A)) + (!eident && eltype(E) !== T) && (E = convert(Matrix{T},E)) + eltype(B) !== T && (B = convert(Matrix{T},B)) + eltype(C) !== T && (C = convert(Matrix{T},C)) + eltype(D) !== T && (D = convert(Matrix{T},D)) + + + (n == 0 || m == 0 || p == 0) && (return A, E, B, C, D, 0, 0, 0) + # save system matrices + Ar = copy(A) + Br = copy(B) + Cr = copy(C) + Dr = copy(D) + Er = copy(E) + ir = 1:n + if eident + if contr + _, _, nr, nuc = sklf_right!(Ar, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol1, rtol = rtol, withQ = false) + if nuc > 0 + ir = 1:nr + # save intermediary results + A = Ar[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore original matrices + Ar = copy(A) + Br = copy(B) + Cr = copy(C) + end + else + nuc = 0 + end + if obs + _, _, no, nuo = sklf_left!(view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:); fast = fast, atol1 = atol1, atol2 = atol1, rtol = rtol, withQ = false) + if nuo > 0 + ir = ir[end-no+1:end] + else + # restore saved matrices + Ar[ir,ir] = A + Br[ir,:] = B + Cr[:,ir] = C + end + else + nuo = 0 + end + return Ar[ir,ir], emat ? Er[ir,ir] : I, Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0 + else + if contr + _, _, _, nr, nfuc, niuc = sklf_right!(Ar, Er, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1, rtol = rtol, withQ = false, withZ = false) + nuc = nfuc+niuc + if nuc > 0 + ir = 1:nr + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore original matrices + Ar = copy(A) + Er = copy(E) + Br = copy(B) + Cr = copy(C) + end + else + nuc = 0 + end + if obs + _, _, _, no, nfuo, niuo = sklf_left!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Br,ir,:); + fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1, + rtol = rtol, withQ = false, withZ = false) + nuo = nfuo+niuo + if nuo > 0 + ir = ir[end-no+1:end] + # save intermediary results + A = Ar[ir,ir] + E = Er[ir,ir] + B = Br[ir,:] + C = Cr[:,ir] + else + # restore saved matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + else + nuo = 0 + end + if noseig + rE, rA22 = _svdlikeAE!(view(Ar,ir,ir), view(Er,ir,ir), nothing, nothing, view(Br,ir,:), view(Cr,:,ir), + fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false) + if rA22 > 0 + i1 = ir[1:rE] + i2 = ir[rE+1:rE+rA22] + # make A22 = I + fast ? (A22 = UpperTriangular(Ar[i2,i2])) : (A22 = Diagonal(Ar[i2,i2])) + ldiv!(A22,view(Ar,i2,i1)) + ldiv!(A22,view(Br,i2,:)) + # apply simplified residualization formulas + Dr -= Cr[:,i2]*Br[i2,:] + Br[i1,:] -= Ar[i1,i2]*Br[i2,:] + Cr[:,i1] -= Cr[:,i2]*Ar[i2,i1] + Ar[i1,i1] -= Ar[i1,i2]*Ar[i2,i1] + ir = [i1; ir[rE+rA22+1:end]] + else + # restore saved matrices + Ar[ir,ir] = A + Er[ir,ir] = E + Br[ir,:] = B + Cr[:,ir] = C + end + return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, rA22 + else + return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0 + end + end +end +""" + sklf_left!(A, E, C, B; fast = true, atol1 = 0, atol2 = 0, atol3 = 0, rtol, withQ = true, withZ = true) -> (Q, Z, μl, no, nfuo, niuo) + +Reduce the partitioned full column rank linear pencil + + M - λN = | A-λE | n + | C | p + n + +to an equivalent form `F - λG = diag(Q',I)*(M - λN)*Z` using orthogonal or unitary transformation matrices `Q` and `Z` +such that `M - λN` is transformed into a Kronecker-like form exhibiting its +infinite, finite and left Kronecker structures (also known as the generalized observability staircase form): + + | Aiuo-λEiuo * | * | niuo + | 0 Afuo-λEfuo | * | nfuo + | At-λEt | | 0 0 | Ao-λEo | no + |--------| = |------------------------|--------| + | Ct | | 0 0 | Co | p + niuo nfuo no + +`Ct = C*Z`, `At = Q'*A*Z` and `Et = Q'*E*Z` are returned in `C`, `A` and `E`, respectively, +and `Q'*B` is returned in `B` (unless `B = missing'). + +The subpencil `| Ao-λEo |` has full column rank `no`, is in a staircase form, and contains the left Kronecker indices of `M - λN`. + `| Co |` +The `nl`-dimensional vector `μl` contains the row and column dimensions of the blocks +of the staircase form such that `i`-th block has dimensions `μl[nl-i] x μl[nl-i+1]` (with μl[0] = p) and +has full column rank. The difference `μl[nl-i]-μl[nl-i+1]` for `i = 1, 2, ..., nl` is the number of elementary Kronecker blocks +of size `i x (i-1)`. + +The `niuo x niuo` subpencil `Aiuo-λEiuo` contains the infinite eigenvalues of `M - λN` (also called the unobservable infinite eigenvalues of `A-λE`). + +The `nfuo x nfuo` subpencil `Afuo-λEfuo` contains the finite eigenvalues of `M - λN` (also called the unobservable finite eigenvalues of `A-λE`). + +The keyword arguments `atol1`, `atol2`, `atol3` and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, +the absolute tolerance for the nonzero elements of `C`, and +the relative tolerance for the nonzero elements of `A`, `E` and `C`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true`. +Otherwise, `Z` is set to `nothing`. +""" +function sklf_left!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, C::AbstractMatrix{T1}, B::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, + atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))), atol3::Real = zero(real(eltype(C))), + rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + n = LinearAlgebra.checksquare(A) + n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) + p, n1 = size(C) + n == n1 || throw(DimensionMismatch("A and C must have the same number of columns")) + (!ismissing(B) && n !== size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows")) + + maxpn = max(p,n) + μl = Vector{Int}(undef,maxpn) + nfu = 0 + niu = 0 + tol1 = atol1 + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing + + # fast returns for null dimensions + if p == 0 && n == 0 + return Q, Z, μl, n, nfu, niu + elseif n == 0 + μl[1] = 0 + return Q, Z, μl[1:1], n, nfu, niu + end + + tolA = max(atol1, rtol*opnorm(A,1)) + tolC = max(atol3, rtol*opnorm(C,Inf)) + + ρ1 = _sreduceC!(A, E, C, Z, tolC; fast = fast, withZ = withZ) + ρ1 == n && (return Q, Z, [ρ1], n, nfu, niu) + + # reduce to basic form + n1, m1, p1 = _preduceBF!(A, E, Q, Z, B, missing; ctrail = ρ1, atol = atol2, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ) + + mrinf = 0 + nrinf = 0 + roff = 0 + coff = 0 + rtrail = 0 + ctrail = ρ1 + niu = 0 + while m1 > 0 + # Steps 1 & 2: Standard algorithm PREDUCE + τ, ρ = _preduce1!(n1, m1, p1, A, E, Q, Z, tolA, B, missing; fast = fast, + roff = roff, coff = coff, ctrail = ctrail, withQ = withQ, withZ = withZ) + ρ+τ == m1 || error("| C' | A'-λE' | has no full row rank") + roff += m1 + coff += m1 + niu += m1 + n1 -= ρ + m1 = ρ + p1 -= τ + end + + if ρ1 > 0 + j = 1 + μl[1] = ρ1 + else + return Q, Z, μl[1:0], 0, n1, niu + end + + no = ρ1 + nfu = n1 + while p1 > 0 + # Step 3: Particular form of the dual algorithm PREDUCE + ρ = _preduce4!(nfu, 0, p1, A, E, Q, Z, tolA, B, missing; fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail, + withQ = withQ, withZ = withZ) + ρ == 0 && break + j += 1 + μl[j] = ρ + rtrail += p1 + ctrail += ρ + no += ρ + nfu -= ρ + p1 = ρ + end + + return Q, Z, reverse(μl[1:j]), no, nfu, niu +end +""" + sklf_leftfin!(A, E, C, B; fast = true, atol1 = 0, atol2 = 0, + rtol, withQ = true, withZ = true) -> (Q, Z, μl, no, nuo) + +Reduce the partitioned full column rank linear pencil + + M - λN = | A-λE | n + | C | p + n + +with `A-λE` regular, to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Z)` +using orthogonal or unitary transformation matrices `Q` and `Z` +such that `M - λN` is transformed into a staircase form exhibiting the separation of its finite +eigenvalues: + + | Auo-λEuo | * | nuo + | At-λEt | | 0 | Ao-λEo | no + |--------| = |-----------|--------| + | Ct | | 0 | Co | p + nuo no + +`Ct = C*Z`, `At = Q'*A*Z` and `Et = Q'*E*Z` are returned in `C`, `A` and `E`, respectively, +and `Q'*B` is returned in `B` (unless `B = missing'). +The resulting `Et` is upper triangular. If `E` is already upper triangular, then +the preliminary reduction of `E` to upper triangular form is not performed. + +The subpencil `| Ao-λEo |` has full column rank `no` for all finite values of `λ`, is in a staircase form, + `| Co |` +and, if `E` is nonsingular, contains the left Kronecker indices of `M - λN`. The `nl`-dimensional vector `μl` contains the row and column dimensions of the blocks +of the staircase form such that `i`-th block has dimensions `μl[nl-i] x μl[nl-i+1]` (with μl[0] = p) and +has full column rank. If `E` is nonsingular, the difference `μl[nl-i]-μl[nl-i+1]` for `i = 1, 2, ..., nl` is the number of elementary Kronecker blocks +of size `i x (i-1)`. + +The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `C`, +and the relative tolerance for the nonzero elements of `A` and `C`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true`. +Otherwise, `Z` is set to `nothing`. + +Note: This function, called with reversed input parameters `E` and `A` (i.e., instead `A` and `E`), performs the +separation all infinite and nonzero finite eigenvalues of the pencil `M - λN`. +""" +function sklf_leftfin!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, C::AbstractMatrix{T1}, B::Union{AbstractMatrix{T1},Missing}; + fast::Bool = true, atol1::Real = zero(real(T1)), atol2::Real = zero(real(T1)), + rtol::Real = (min(size(A)...)*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + + n = LinearAlgebra.checksquare(A) + n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) + p, n1 = size(C) + n == n1 || throw(DimensionMismatch("A and C must have the same number of columns")) + (!ismissing(B) && n !== size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows")) + + μl = Vector{Int}(undef,max(n,1)) + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing + + # fast returns for null dimensions + if p == 0 && n == 0 + return Q, Z, μl[1:0], 0, 0 + elseif n == 0 + μl[1] = 0 + return Q, Z, μl[1:1], 0, 0 + end + + # Reduce E to upper triangular form if necessary + istriu(E) || _qrE!(A, E, Q, B; withQ = withQ) + p == 0 && (return Q, Z, μl[1:0], 0, n) + + + tolA = max(atol1, rtol*opnorm(A,1)) + tolC = max(atol2, rtol*opnorm(C,Inf)) + + i = 0 + init = true + rtrail = 0 + ctrail = 0 + no = 0 + nuo = n + while p > 0 && no < n + init = (i == 0) + init ? tol = tolC : tol = tolA + # Step 3: Particular case of the standard algorithm PREDUCE + ρ = _sreduceAEC!(nuo, p, A, E, C, B, Q, Z, tol, fast = fast, init = init, + rtrail = rtrail, ctrail = ctrail, withQ = withQ, withZ = withZ) + ρ == 0 && break + i += 1 + μl[i] = ρ + ctrail += ρ + init || (rtrail += p) + no += ρ + nuo -= ρ + p = ρ + end + + return Q, Z, reverse(μl[1:i]), no, nuo +end +""" + sklf_right!(A, E, B, C; fast = true, atol1 = 0, atol2 = 0, atol3 = 0, rtol, withQ = true, withZ = true) -> (Q, Z, νr, nc, nfuc, niuc) + +Reduce the partitioned full row rank linear pencil + + M - λN = | B | A-λE | n + m n + +with `A-λE` regular, to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Z)` using orthogonal or unitary transformation matrices `Q` and `Z` +such that `M - λN` is transformed into a Kronecker-like form `| Bt | At-λEt |` exhibiting its +right, finite and infinite Kronecker structures (also known as the generalized controllability staircase form): + + | Bc | Ac-λEc * * | nc + | Bt | At-λEt | = |-----|------------------------------| + | 0 | 0 Afuc-λEfuc * | nfuc + | 0 | 0 0 Aiuc-λEiuc | niuc + m nc nfuc niuc + +`Bt = Q'*B`, `At = Q'*A*Z`and `Et = Q'*E*Z` are returned in `B`, `A` and `E`, respectively, and `C*Z` is returned in `C` (unless `C = missing'). + +The subpencil `| Bc | Ac-λEc |` has full row rank `nc`, is in a staircase form, and contains the right Kronecker indices of `M - λN`. +The `nr`-dimensional vector `νr` contains the row and column dimensions of the blocks +of the staircase form `| Bc | Ac-λI |` such that `i`-th block has dimensions `νr[i] x νr[i-1]` (with `νr[0] = m`) and +has full row rank. The difference `νr[i-1]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks +of size `(i-1) x i`. + +The `nfuc x nfuc` subpencil `Afuc-λEfuc` contains the finite eigenvalues of `M - λN` (also called the uncontrollable finite eigenvalues of `A - λE`). + +The `niuc x niuc` subpencil `Aiuc-λEiuc` contains the infinite eigenvalues of `M - λN` (also called the uncontrollable infinite eigenvalues of `A - λE`). + +The keyword arguments `atol1`, `atol2`, , `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `B`, +and the relative tolerance for the nonzero elements of `A`, `E` and `B`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. +The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true` or +`C` is provided. Otherwise, `Z` is set to `nothing`. +""" +function sklf_right!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, + atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))), atol3::Real = zero(real(eltype(B))), + rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + n = LinearAlgebra.checksquare(A) + n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) + n1, m = size(B) + n == n1 || throw(DimensionMismatch("A and B must have the same number of rows")) + (!ismissing(C) && n !== size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns")) + + maxmn = max(m,n) + νr = Vector{Int}(undef,maxmn) + nfu = 0 + niu = 0 + tol1 = atol1 + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing + + # fast returns for null dimensions + if m == 0 && n == 0 + return Q, Z, νr, n, nfu, niu + elseif n == 0 + νr[1] = 0 + return Q, Z, νr[1:1], n, nfu, niu + end + + tolA = max(atol1, rtol*opnorm(A,1)) + tolB = max(atol3, rtol*opnorm(B,1)) + + ρ1 = _sreduceB!(A, E, B, Q, tolB; fast = fast, withQ = withQ) + ρ1 == n && (return Q, Z, [ρ1], n, nfu, niu) + + n1, m1, p1 = _preduceBF!(A, E, Q, Z, missing, C; roff = ρ1, atol = atol2, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ) + + mrinf = ρ1 + nrinf = 0 + rtrail = 0 + ctrail = 0 + niu = 0 + while p1 > 0 + # Step 1 & 2: Dual algorithm PREDUCE + τ, ρ = _preduce2!(n1, m1, p1, A, E, Q, Z, tolA, missing, C; fast = fast, + roff = mrinf, coff = nrinf, rtrail = rtrail, ctrail = ctrail, withQ = withQ, withZ = withZ) + ρ+τ == p1 || error("| B | A-λE | has no full row rank") + ctrail += p1 + rtrail += p1 + niu += p1 + n1 -= ρ + p1 = ρ + m1 -= τ + end + + if ρ1 > 0 + i = 1 + νr[1] = ρ1 + else + return Q, Z, νr[1:0], 0, n1, niu + end + if m1 > 0 + imA11 = 1:n-rtrail + A11 = view(A,imA11,1:n) + E11 = view(E,imA11,1:n) + ismissing(C) ? C1 = missing : (C1 = view(C,:,1:n)) + end + nc = ρ1 + nfu = n1 + while m1 > 0 + # Step 3: Particular case of the standard algorithm PREDUCE + ρ = _preduce3!(nfu, m1, A11, E11, Q, Z, tolA, missing, C1, fast = fast, coff = nrinf, roff = mrinf, ctrail = ctrail, withQ = withQ, withZ = withZ) + ρ == 0 && break + i += 1 + νr[i] = ρ + mrinf += ρ + nrinf += m1 + nc += ρ + nfu -= ρ + m1 = ρ + end + + return Q, Z, νr[1:i], nc, nfu, niu +end +""" + sklf_rightfin!(A, E, B, C; fast = true, atol1 = 0, atol2 = 0, + rtol, withQ = true, withZ = true) -> (Q, Z, νr, nc, nuc) + +Reduce the partitioned full row rank linear pencil + + M - λN = | B | A-λE | n + m n + +with `A-λE` regular, to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Z)` +using orthogonal or unitary transformation matrices `Q` and `Z` +such that `M - λN` is transformed into a staircase form `| Bt | At-λEt |` exhibiting the separation of its finite +eigenvalues: + + | Bc | Ac-λEc * | nc + | Bt | At-λEt | = |-----|-----------------| + | 0 | 0 Auc-λEuc | nuc + m nc nuc + +`Bt = Q'*B`, `At = Q'*A*Z`and `Et = Q'*E*Z` are returned in `B`, `A` and `E`, respectively, and `C*Z` is returned in `C` (unless `C = missing'). +The resulting `Et` is upper triangular. If `E` is already upper triangular, then +the preliminary reduction of `E` to upper triangular form is not performed. + +The subpencil `| Bc | Ac-λEc |` has full row rank `nc` for all finite values of `λ`, is in a staircase form, and, +if E is invertible, contains the right Kronecker indices of `M - λN`. +The `nr`-dimensional vector `νr` contains the row and column dimensions of the blocks +of the staircase form `| Bc | Ac-λI |` such that `i`-th block has dimensions `νr[i] x νr[i-1]` (with `νr[0] = m`) and +has full row rank. If E is invertible, the difference `νr[i-1]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks +of size `(i-1) x i`. + +The `nuc x nuc` subpencil `Auc-λEuc` contains the finite eigenvalues of `M - λN` (also called the uncontrollable finite eigenvalues of `A - λE`). +If E is singular, `Auc-λEuc` may also contain a part of the infinite eigenvalues of `M - λN` (also called the uncontrollable infinite eigenvalues of `A - λE`). + +The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `B`, +and the relative tolerance for the nonzero elements of `A` and `B`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true`. +Otherwise, `Z` is set to `nothing`. + +Note: This function, called with reversed input parameters `E` and `A` (i.e., instead `A` and `E`), performs the +separation all infinite and nonzero finite eigenvalues of the pencil `M - λN`. +""" +function sklf_rightfin!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::Union{AbstractMatrix{T1},Missing}; + fast::Bool = true, atol1::Real = zero(real(T1)), atol2::Real = zero(real(T1)), + rtol::Real = (min(size(A)...)*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + n = LinearAlgebra.checksquare(A) + n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) + n1, m = size(B) + n == n1 || throw(DimensionMismatch("A and B must have the same number of rows")) + (!ismissing(C) && n !== size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns")) + + νr = Vector{Int}(undef,max(n,1)) + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing + + # fast returns for null dimensions + if m == 0 && n == 0 + return Q, Z, νr[1:0], 0, 0 + elseif n == 0 + νr[1] = 0 + return Q, Z, νr[1:1], 0, 0 + end + + # Reduce E to upper triangular form if necessary + istriu(E) || _qrE!(A, E, Q, B; withQ = withQ) + m == 0 && (return Q, Z, νr[1:0], 0, n) + + tolA = max(atol1, rtol*opnorm(A,1)) + tolB = max(atol2, rtol*opnorm(B,1)) + + i = 0 + init = true + roff = 0 + coff = 0 + nc = 0 + nuc = n + while m > 0 && nc < n + init = (i == 0) + init ? tol = tolB : tol = tolA + # Step 3: Particular case of the standard algorithm PREDUCE + ρ = _sreduceBAE!(nuc, m, A, E, B, C, Q, Z, tol, fast = fast, init = init, roff = roff, coff = coff, withQ = withQ, withZ = withZ) + ρ == 0 && break + i += 1 + νr[i] = ρ + roff += ρ + init || (coff += m) + nc += ρ + nuc -= ρ + m = ρ + end + + return Q, Z, νr[1:i], nc, nuc +end +""" + sklf_right!(A, B, C; fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true) -> (Q, νr, nc, nuc) + +Reduce the partitioned full row rank linear pencil + + M - λN = | B | A-λI | n + m n + +to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Q)` using orthogonal or unitary transformation matrix `Q` +such that `M - λN` is transformed into a Kronecker-like form `| Bt | At-λI |` exhibiting its +right and finite Kronecker structures (also known as the controllability staircase form): + + | Bc | Ac-λI * | nc + | Bt | At-λI | = |-----|----------------| + | 0 | 0 Auc-λI | nuc + m nc nuc + +`Bt = Q'*B` and `At = Q'*A*Q` are returned in `B` and `A`, respectively, and `C*Q` is returned in `C` (unless `C = missing'). + +The subpencil `| Bc | Ac-λI |` has full row rank `nc`, is in a staircase form, and contains the right Kronecker indices of `M - λN`. +The `nr`-dimensional vector `νr` contains the row and column dimensions of the blocks +of the staircase form `| Bc | Ac-λI |` such that `i`-th block has dimensions `νr[i] x νr[i-1]` (with `νr[0] = m`) and +has full row rank. The difference `νr[i-1]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks +of size `(i-1) x i`. + +The `nuc x nuc` matrix `Auc` contains the (finite) eigenvalues of `M - λN` (also called the uncontrollable eigenvalues of `A`). + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `B`, and the relative tolerance +for the nonzero elements of `A` and `B`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. +The performed orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +""" +function sklf_right!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, + atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), + rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)), + withQ::Bool = true) where T1 <: BlasFloat + n = LinearAlgebra.checksquare(A) + n1, m = size(B) + n == n1 || throw(DimensionMismatch("A and B must have the same number of rows")) + (!ismissing(C) && n !== size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns")) + + νr = Vector{Int}(undef,max(n,1)) + nu = 0 + tol1 = atol1 + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + + # fast returns for null dimensions + if m == 0 && n == 0 + return Q, νr[1:0], 0, 0 + elseif n == 0 + νr[1] = 0 + return Q, νr[1:1], 0, 0 + elseif m == 0 + return Q, νr[1:0], 0, n + end + + tolA = max(atol1, rtol*opnorm(A,1)) + tolB = max(atol2, rtol*opnorm(B,1)) + + i = 0 + init = true + roff = 0 + coff = 0 + nc = 0 + nu = n + while m > 0 && nc < n + init = (i == 0) + init ? tol = tolB : tol = tolA + # Step 3: Particular case of the standard algorithm PREDUCE + ρ = _sreduceBA!(nu, m, A, B, C, Q, tol, fast = fast, init = init, roff = roff, coff = coff, withQ = withQ) + ρ == 0 && break + i += 1 + νr[i] = ρ + roff += ρ + init || (coff += m) + nc += ρ + nu -= ρ + m = ρ + end + + return Q, νr[1:i], nc, nu +end +""" + sklf_left!(A, C, B; fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true) -> (Q, μl, no, nuo) + +Reduce the partitioned full column rank linear pencil + + M - λN = | A-λI | n + | C | p + n + +to an equivalent form `F - λL = diag(Q',I)*(M - λN)*Q` using orthogonal or unitary transformation matrix `Q` +such that `M - λN` is transformed into a Kronecker-like form exhibiting its +finite and left Kronecker structures (also known as the observability staircase form): + + | Auo-λI | * | nuo + | At-λI | | 0 | Ao-λI | no + |-------| = |---------|-------| + | Ct | | 0 | Co | p + nuo no + +`Ct = C*Q` and `At = Q'*A*Q` are returned in `C` and `A`, respectively, and `Q'*B` is returned in `B` (unless `B = missing'). + +The subpencil `| Ao-λI; Co |` has full column rank `no`, is in a staircase form, +and contains the left Kronecker indices of `M - λN`. + +The `nl`-dimensional vector `μl` contains the row and column dimensions of the blocks +of the staircase form such that `i`-th block has dimensions `μl[nl-i] x μl[nl-i+1]` (with μl[0] = p) and +has full column rank. The difference `μl[nl-i]-μl[nl-i+1]` for `i = 1, 2, ..., nl` is the number of elementary Kronecker blocks +of size `i x (i-1)`. + +The `nuo x nuo` matrix `Auo` contains the (finite) eigenvalues of `M - λN` (also called the unobservable eigenvalues of `A`). + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `C`, and the relative tolerance +for the nonzero elements of `A` and `C`. +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`. +Otherwise, `Q` is set to `nothing`. +""" +function sklf_left!(A::AbstractMatrix{T1}, C::AbstractMatrix{T1}, B::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), + rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)), + withQ::Bool = true) where T1 <: BlasFloat + n = LinearAlgebra.checksquare(A) + p, n1 = size(C) + n == n1 || throw(DimensionMismatch("A and C must have the same number of columns")) + (!ismissing(B) && n !== size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows")) + + μl = Vector{Int}(undef,max(n,p)) + nu = 0 + tol1 = atol1 + withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing + + # fast returns for null dimensions + if p == 0 && n == 0 + return Q, μl[1:0], 0, 0 + elseif n == 0 + μl[1] = 0 + return Q, μl[1:1], 0, 0 + elseif p == 0 + return Q, μl[1:0], 0, n + end + + tolA = max(atol1, rtol*opnorm(A,1)) + tolC = max(atol2, rtol*opnorm(C,Inf)) + + i = 0 + init = true + rtrail = 0 + ctrail = 0 + no = 0 + nu = n + while p > 0 && no < n + init = (i == 0) + init ? tol = tolC : tol = tolA + # Step 3: Particular case of the standard algorithm PREDUCE + ρ = _sreduceAC!(nu, p, A, C, B, Q, tol, fast = fast, init = init, rtrail = rtrail, ctrail = ctrail, withQ = withQ) + ρ == 0 && break + i += 1 + μl[i] = ρ + ctrail += ρ + init || (rtrail += p) + no += ρ + nu -= ρ + p = ρ + end + + return Q, reverse(μl[1:i]), no, nu +end diff --git a/src/regtools.jl b/src/regtools.jl index 68676ea..dc39f75 100644 --- a/src/regtools.jl +++ b/src/regtools.jl @@ -1,38 +1,236 @@ """ - isregular(M, N; atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> Bool + _qrE!(A, E, Q, B; withQ = true) -Test whether the linear pencil `M-λN` is regular (i.e., det(M-λN) !== 0). The underlying computational procedure -reduces the pencil `M-λN` to an appropriate Kronecker-like form (KLF), which provides information on the rank of `M-λN`. +Reduce the linear regular pencil `A - λE` to an equivalent form `A1 - λE1 = Q1'*(A - λE)` using an +orthogonal or unitary transformation matrix `Q1` such that the transformed matrix `E1` is upper triangular. +The reduction is performed using the QR-decomposition of E. + +The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1` +if `withQ = true`. Otherwise, `Q` is unchanged. + +`Q1'*B` is returned in `B` unless `B = missing`. +""" +function _qrE!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, Q::Union{AbstractMatrix{T1},Nothing}, + B::Union{AbstractMatrix{T1},Missing} = missing; withQ::Bool = true) where T1 <: BlasFloat + + # fast return for dimensions 0 or 1 + size(A,1) <= 1 && return + T = eltype(A) + T <: Complex ? tran = 'C' : tran = 'T' + # compute in-place the QR-decomposition E = Q1*E1 + _, τ = LinearAlgebra.LAPACK.geqrf!(E) + # A <- Q1'*A + LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,A) + # Q <- Q*Q1 + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E,τ,Q) + # B <- Q1'*B + ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,B) + triu!(E) +end +""" + _svdlikeAE!(A, E, Q, Z, B, C; fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (rE, rA22) + +Reduce the linear regular pencil `A - λE` to an equivalent form `A1 - λE1 = Q1'*(A - λE)*Z1` using +orthogonal or unitary transformation matrices `Q1` and `Z1` such that the transformed matrices `A1` and `E1` +are in the following SVD-like coordinate form + + | A11-λE11 | A12 | A13 | + |----------|-------|-------|, + A1 - λE1 = | A21 | A22 | 0 | + |----------|-------|-------| + | A31 | 0 | 0 | + +where the `rE x rE` matrix `E11` and `rA22 x rA22` matrix `A22` are nosingular, and `E11` and `A22` are upper triangular, +if `fast = true`, and diagonal, if `fast = false`. +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, and the relative tolerance +for the nonzero elements of `A` and `E`. + +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1` +if `withQ = true`. Otherwise, `Q` is unchanged. +The performed right orthogonal or unitary transformations Z1 are accumulated in the matrix `Z <- Z*Z1` if `withZ = true`. +Otherwise, `Z` is unchanged. + +`Q1'*B` is returned in `B` unless `B = missing` and `C*Z1` is returned in `C` unless `C = missing` . + +""" +function _svdlikeAE!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, + B::Union{AbstractMatrix{T1},Missing} = missing, C::Union{AbstractMatrix{T1},Missing} = missing; + fast::Bool = true, atol1::Real = zero(real(T1)), atol2::Real = zero(real(T1)), + rtol::Real = (size(A,1)*eps(real(float(one(T1)))))*iszero(min(atol1,atol2)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + + # fast returns for null dimensions + n = size(A,1) + if n == 0 + return 0, 0 + end + T = eltype(A) + T <: Complex ? tran = 'C' : tran = 'T' + + if fast + # compute in-place the QR-decomposition E*P1 = Q1*[E1;0] with column pivoting + _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(E) + tol = max(atol2, rtol*abs(E[1,1])) + rE = count(x -> x > tol, abs.(diag(E))) + n2 = n-rE + i1 = 1:rE + # A <- Q1'*A + LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,A) + # A <- A*P1 + A[:,:] = A[:,jpvt] + # Q <- Q*Q1 + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E,τ,Q) + # Z <- Z*P1 + withZ && (Z[:,:] = Z[:,jpvt]) + # B <- Q1'*B + ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,B) + # E = [E1;0] + E[:,:] = [ triu(E[i1,:]); zeros(T,n2,n) ] + # C <- C*P1 + ismissing(C) || (C[:,:] = C[:,jpvt]) + + # compute in-place the complete orthogonal decomposition E*Z1 = [E11 0; 0 0] with E11 nonsingular and UT + E1 = view(E,i1,:) + _, tau = LinearAlgebra.LAPACK.tzrzf!(E1) + # A <- A*Z1 + LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,A) + withZ && LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,Z) + # C <- C*Z1 + ismissing(C) || LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,C); + E1[:,:] = [ triu(E[i1,i1]) zeros(T,rE,n2) ] + if n2 > 0 + # assume + # A = [A11 A12] + # [A21 A22] + # compute in-place the QR-decomposition A22*P2 = Q2*[R2;0] with column pivoting + i22 = rE+1:n + A22 = view(A,i22,i22) + _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(A22) + tolA = max(atol2, rtol*opnorm(A,1)) + rA22 = count(x -> x > tolA, abs.(diag(A22))) + n3 = n2-rA22 + i2 = rE+1:rE+rA22 + i3 = rE+rA22+1:n + # A21 <- Q2'*A21 + LinearAlgebra.LAPACK.ormqr!('L',tran,A22,τ,view(A,i22,i1)) + # A12 <- A12*P1 + A[i1,i22] = A[i1,i22[jpvt]] + # Q <- Q*Q2 + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',A22,τ,view(Q,:,i22)) + # Z <- Z*P2 + withZ && (Z[:,i22] = Z[:,i22[jpvt]]) + # B <- Q2'*B + ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,A22,τ,view(B,i22,:)) + # A22 = [R2;0] + A22[:,:] = [ triu(A[i2,i22]); zeros(T,n3,n2) ] + # R <- R*P1 + ismissing(C) || (C[:,i22] = C[:,i22[jpvt]]) + + # compute in-place the complete orthogonal decomposition R2*Z2 = [A2 0; 0 0] with A22 nonsingular and UT + A2 = view(A,i2,i22) + _, tau = LinearAlgebra.LAPACK.tzrzf!(A2) + # A12 <- A12*Z2 + LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(A,i1,i22)) + withZ && LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(Z,:,i22)) + # C <- C*Z2 + ismissing(C) || LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(C,:,i22)); + A2[:,:] = [ triu(A[i2,i2]) zeros(T,rA22,n3) ] + else + rA22 = 0 + end + else + # compute the complete orthogonal decomposition of E using the SVD-decomposition + U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',E) + tolE = max(atol2, rtol*S[1]) + rE = count(x -> x > tolE, S) + n2 = n-rE + # A <- A*V + A[:,:] = A[:,:]*Vt' + # A <- U'*A + A[:,:] = U'*A[:,:] + # Q <- Q*U + withQ && (Q[:,:] = Q[:,:]*U) + # Z <- Q*V + withZ && (Z[:,:] = Z[:,:]*Vt') + # B <- U'*B + ismissing(B) || (B[:,:] = U'*B[:,:]) + # C <- C*V + ismissing(C) || (C[:,:] = C[:,:]*Vt') + E[:,:] = [ Diagonal(S[1:rE]) zeros(T,rE,n2) ; zeros(T,n2,n) ] + if n2 > 0 + # assume + # A = [A11 A12] + # [A21 A22] + # compute the complete orthogonal decomposition of A22 using the SVD-decomposition + i1 = + i22 = rE+1:n + A22 = view(A,i22,i22) + U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',A22) + tolA = max(atol2, rtol*opnorm(A,1)) + rA22 = count(x -> x > tolA, S) + n3 = n2-rA22 + i1 = 1:rE + i2 = rE+1:rE+rA22 + i3 = rE+rA22+1:n + # A12 <- A12*V + A[i1,i22] = A[i1,i22]*Vt' + # A21 <- U'*A21 + A[i22,i1] = U'*A[i22,i1] + # Q <- Q*U + withQ && (Q[:,i22] = Q[:,i22]*U) + # Z <- Q*V + withZ && (Z[:,i22] = Z[:,i22]*Vt') + # B <- U'*B + ismissing(B) || (B[i22,:] = U'*B[i22,:]) + # C <- C*V + ismissing(C) || (C[:,i22] = C[:,i22]*Vt') + A22[:,:] = [ Diagonal(S[1:rA22]) zeros(T,rA22,n3); zeros(T,n3,n2) ] + else + rA22 = 0 + end + end + return rE, rA22 +end +""" + isregular(A, E; atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> Bool + +Test whether the linear pencil `A-λE` is regular (i.e., det(A-λE) !== 0). The underlying computational procedure +reduces the pencil `A-λE` to an appropriate Kronecker-like form (KLF), which provides information on the rank of `M-λN`. The keyword arguements `atol1`, `atol2` and `rtol` specify the absolute tolerance for the nonzero -elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively. -The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the -machine epsilon of the element type of `M`. +elements of `A`, the absolute tolerance for the nonzero elements of `E`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively. +The default relative tolerance is `n*ϵ`, where `n` is the size of `A`, and `ϵ` is the +machine epsilon of the element type of `A`. """ -function isregular(M ::AbstractMatrix, N::AbstractMatrix; atol1::Real = zero(eltype(M)), atol2::Real = zero(eltype(M)), - rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2))) +function isregular(A::AbstractMatrix, E::AbstractMatrix; atol1::Real = zero(eltype(A)), atol2::Real = zero(eltype(A)), + rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2))) - mM, nM = size(M) - (mM,nM) == size(N) || throw(DimensionMismatch("M and N must have the same dimensions")) - mM == nM || (return false) - isa(M,Adjoint) && (M = copy(M)) - isa(N,Adjoint) && (N = copy(N)) - T = promote_type(eltype(M), eltype(N)) + mA, nA = size(A) + (mA,nA) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions")) + mA == nA || (return false) + isa(A,Adjoint) && (A = copy(A)) + isa(E,Adjoint) && (E = copy(E)) + T = promote_type(eltype(A), eltype(E)) T <: BlasFloat || (T = promote_type(Float64,T)) - eltype(M) == T || (M = convert(Matrix{T},M)) - eltype(N) == T || (N = convert(Matrix{T},N)) + eltype(A) == T || (A = convert(Matrix{T},A)) + eltype(E) == T || (E = convert(Matrix{T},E)) Q = nothing Z = nothing # Step 0: Reduce to the standard form - n, m, p = _preduceBF!(M, N, Q, Z; atol = atol2, rtol = rtol, fast = false, withQ = false, withZ = false) + n, m, p = _preduceBF!(A, E, Q, Z; atol = atol2, rtol = rtol, fast = false, withQ = false, withZ = false) mrinf = 0 nrinf = 0 - tol1 = max(atol1, rtol*opnorm(M,1)) + tol1 = max(atol1, rtol*opnorm(A,1)) while m > 0 # Steps 1 & 2: Standard algorithm PREDUCE - τ, ρ = _preduce1!(n, m, p, M, N, Q, Z, tol1; fast = false, + τ, ρ = _preduce1!(n, m, p, A, E, Q, Z, tol1; fast = false, roff = mrinf, coff = nrinf, withQ = false, withZ = false) ρ+τ == m || (return false) mrinf += ρ+τ @@ -44,16 +242,16 @@ function isregular(M ::AbstractMatrix, N::AbstractMatrix; atol1::Real = zero(elt return true end """ - fisplit(A, E, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (F, G, H, L, Q, Z, ν, nf, ni) + fisplit(A, E, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Bt, Ct, Q, Z, ν, nf, ni) -Reduce the linear regular pencil `A - λE` to an equivalent form `F - λG = Q'*(A - λE)*Z` using -orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in one of the +Reduce the linear regular pencil `A - λE` to an equivalent form `At - λEt = Q'*(A - λE)*Z` using +orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the following block upper-triangular forms: (1) if `finite_infinite = false`, then | Ai-λEi | * | - F - λG = |--------|---------|, + At - λEt = |--------|---------|, | O | Af-λEf | where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary @@ -69,7 +267,7 @@ The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite (2) if `finite_infinite = true`, then | Af-λEf | * | - F - λG = |--------|--------|, + At - λEt = |--------|--------|, | O | Ai-λEi | where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, @@ -94,7 +292,8 @@ Otherwise, `Q` is set to `nothing`. The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`. Otherwise, `Z` is set to `nothing`. -`Q'*B` is returned in `B` unless `B = missing` and `C*Z` is returned in `C` unless `C = missing` . +`Bt = Q'*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`, +unless `C = missing`, in which case `Ct = missing` is returned . """ function fisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,Missing}, C::Union{AbstractMatrix,Missing}; @@ -123,16 +322,89 @@ function fisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,M withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing) withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing) - # Step 0: Reduce to the standard form - n1, m1, p1 = _preduceBF!(A, E, Q, Z, B, C; atol = atol2, rtol = rtol, fast = fast) + ν, nf, ni = fisplit!(A, E, Q, Z, B, C; + fast = fast, finite_infinite = finite_infinite, + atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ) + + + + return A, E, B, C, Q, Z, ν, nf, ni +end +""" + fisplit!(A, E, B, C, Q, Z; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (ν, nf, ni) + +Reduce the linear regular pencil `A - λE` to an equivalent form `At - λEt = Q1'*(A - λE)*Z1` using +orthogonal or unitary transformation matrices `Q1` and `Z1` such that the transformed matrices `At` and `Et` are in one of the +following block upper-triangular forms: + +(1) if `finite_infinite = false`, then + + | Ai-λEi | * | + At - λEt = |--------|---------|, + | O | Af-λEf | + +where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary +divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`. + +The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and `Ei` nilpotent. +The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks +of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`. +The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i` +(with `ν[nb+1] = 0`). + + +(2) if `finite_infinite = true`, then + + | Af-λEf | * | + At - λEt = |--------|--------|, + | O | Ai-λEi | + +where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, +contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi` +contains the infinite elementary divisors. + +The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and `Ei` nilpotent. +The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks +of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`. +The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary +divisors of degree `j` (with `ν[0] = 0`). + +The reduced matrices `At` and `Et` are returned in `A` and `E`, respectively, +while `Q1'*B` is returned in `B`, unless `B = missing`, and `C*Z1`, is returned in `C`, +unless `C = missing`. + +The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the +nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, and the relative tolerance +for the nonzero elements of `A` and `E`. + +The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting +if `fast = true` or the more reliable SVD-decompositions if `fast = false`. + +The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1` +if `withQ = true`. Otherwise, `Q` is unchanged. +The performed right orthogonal or unitary transformations Z1 are accumulated in the matrix `Z <- Z*Z1` +if `withZ = true`. Otherwise, `Z` is unchanged. + +""" +function fisplit!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, + B::Union{AbstractMatrix{T1},Missing} = missing, C::Union{AbstractMatrix{T1},Missing} = missing; + fast::Bool = true, finite_infinite::Bool = false, + atol1::Real = zero(real(T1)), atol2::Real = zero(real(T1)), + rtol::Real = (size(A,1)*eps(real(float(one(T1)))))*iszero(min(atol1,atol2)), + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat - ν = Vector{Int}(undef,n) - # fast returns for null dimensions + n = size(A,1) + ν = Vector{Int}(undef,n) if n == 0 - return A, E, B, C, Q, Z, ν, 0, 0 + return ν, 0, 0 end + T = eltype(A) + # Step 0: Reduce to the standard form + n1, m1, p1 = _preduceBF!(A, E, Q, Z, B, C; atol = atol2, rtol = rtol, fast = fast) + tolA = max(atol1, rtol*opnorm(A,1)) if finite_infinite @@ -160,7 +432,7 @@ function fisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,M p1 = ρ m1 -= τ end - return A, E, B, C, Q, Z, reverse(ν[1:i]), nf, ni + return reverse(ν[1:i]), nf, ni else # Reduce A-λE to the to the Kronecker-like form by splitting the infinite-finite structures @@ -187,6 +459,6 @@ function fisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,M p1 -= τ end - return A, E, B, C, Q, Z, ν[1:i], nf, ni + return ν[1:i], nf, ni end end diff --git a/src/sklfapps.jl b/src/sklfapps.jl index ef0186a..1096ad2 100644 --- a/src/sklfapps.jl +++ b/src/sklfapps.jl @@ -268,7 +268,7 @@ function spkstruct(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,Uni end """ - sprank(A, E, B, C, D; fastrank = true, atol1::Real=0, atol2::Real=0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) + sprank(A, E, B, C, D; fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ) Compute the normal rank of the structured linear matrix pencil `M - λN` diff --git a/src/sklftools.jl b/src/sklftools.jl index 6ae62c3..6d78607 100644 --- a/src/sklftools.jl +++ b/src/sklftools.jl @@ -1,5 +1,5 @@ """ - sreduceBF(A, E, B, C, D; fast = true, withQ = true, withZ = true) -> F, G, N, Q, Z, n, m, p + sreduceBF(A, E, B, C, D; fast = true, atol = 0, rtol, withQ = true, withZ = true) -> F, G, N, Q, Z, n, m, p Reduce the partitioned linear pencil `M - λN` @@ -94,7 +94,7 @@ function sreduceBF(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,Uni isa(E,Adjoint) && (E = copy(E)) ndx, nx = size(A) T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D)) - eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and M must have the same dimensions")) + eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions")) ny, nu = size(D) (ndx,nu) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions")) (ny,nx) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions")) @@ -413,419 +413,3 @@ function sklf_left(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,Uni return M, N, Q, Z, ν, μ, nf, νl, μl end -""" - sklf_left!(A, E, C, L; fast = true, atol1 = 0, atol2 = 0, atol3 = 0, rtol, withQ = true, withZ = true) -> (Q, Z, νl, no, nfuo, niuo) - -Reduce the partitioned full column rank linear pencil - - M - λN = | A-λE | n - | C | p - n - -to an equivalent form `F - λG = diag(Q',I)*(M - λN)*Z` using orthogonal or unitary transformation matrices `Q` and `Z` -such that `M - λN` is transformed into a Kronecker-like form exhibiting its -infinite, finite and left Kronecker structures (also known as the generalized observability staircase form): - - | Aiuo-λEiuo * | * | - | 0 Afuo-λEfuo | * | - | At-λEt | | 0 0 | Ao-λEo | - |--------| = |------------------------|--------| - | Ct | | 0 0 | Co | - -`Ct = C*Z`, `At = Q'*A*Z` and `Et = Q'*E*Z` are returned in `C`, `A` and `E`, respectively, -and `Q'*B` is returned in `B` (unless `B = missing'). - -The subpencil `| Ao-λEo |` has full column rank `no`, is in a staircase form, and contains the left Kronecker indices of `M - λN`. - `| Co |` -The `nl`-dimensional vector `μl` contains the row and column dimensions of the blocks -of the staircase form such that `i`-th block has dimensions `μl[nl-i] x μl[nl-i+1]` (with μl[0] = p) and -has full column rank. The difference `μl[nl-i]-μl[nl-i+1]` for `i = 1, 2, ..., nl` is the number of elementary Kronecker blocks -of size `i x (i-1)`. - -The `niuo x niuo` subpencil `Aiuo-λEiuo` contains the infinite eigenvalues of `M - λN` (also called the unobservable infinite eigenvalues of `A-λE`). - -The `nfuo x nfuo` subpencil `Afuo-λEfuo` contains the finite eigenvalues of `M - λN` (also called the unobservable finite eigenvalues of `A-λE`). - -The keyword arguments `atol1`, `atol2`, `atol3` and `rtol`, specify, respectively, the absolute tolerance for the -nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, -the absolute tolerance for the nonzero elements of `C`, and -the relative tolerance for the nonzero elements of `A`, `E` and `C`. -The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting -if `fast = true` or the more reliable SVD-decompositions if `fast = false`. - -The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. -Otherwise, `Q` is set to `nothing`. -The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true` or -`C` is provided. Otherwise, `Z` is set to `nothing`. -""" -function sklf_left!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, C::AbstractMatrix{T1}, B::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, - atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))), atol3::Real = zero(real(eltype(C))), - rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)), - withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat - n = LinearAlgebra.checksquare(A) - n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) - p, n1 = size(C) - n == n1 || throw(DimensionMismatch("A and C must have the same number of columns")) - (!ismissing(B) && n !== size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows")) - - maxpn = max(p,n) - μl = Vector{Int}(undef,maxpn) - nfu = 0 - niu = 0 - tol1 = atol1 - withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing - withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing - - # fast returns for null dimensions - if p == 0 && n == 0 - return Q, Z, μl, n, nfu, niu - elseif n == 0 - μl[1] = 0 - return Q, Z, μl[1:1], n, nfu, niu - end - - tolA = max(atol1, rtol*opnorm(A,1)) - tolC = max(atol3, rtol*opnorm(C,Inf)) - - ρ1 = _sreduceC!(A, E, C, Z, tolC; fast = fast, withZ = withZ) - ρ1 == n && (return Q, Z, [ρ1], n, nfu, niu) - - # reduce to basic form - n1, m1, p1 = _preduceBF!(A, E, Q, Z, B, missing; ctrail = ρ1, atol = atol2, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ) - - mrinf = 0 - nrinf = 0 - roff = 0 - coff = 0 - rtrail = 0 - ctrail = ρ1 - niu = 0 - while m1 > 0 - # Steps 1 & 2: Standard algorithm PREDUCE - τ, ρ = _preduce1!(n1, m1, p1, A, E, Q, Z, tolA, B, missing; fast = fast, - roff = roff, coff = coff, ctrail = ctrail, withQ = withQ, withZ = withZ) - ρ+τ == m1 || error("| C' | A'-λE' | has no full row rank") - roff += m1 - coff += m1 - niu += m1 - n1 -= ρ - m1 = ρ - p1 -= τ - end - - if ρ1 > 0 - j = 1 - μl[1] = ρ1 - else - return Q, Z, μl[1:0], 0, n1, niu - end - - no = ρ1 - nfu = n1 - while p1 > 0 - # Step 3: Particular form of the dual algorithm PREDUCE - ρ = _preduce4!(nfu, 0, p1, A, E, Q, Z, tolA, B, missing; fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail, - withQ = withQ, withZ = withZ) - ρ == 0 && break - j += 1 - μl[j] = ρ - rtrail += p1 - ctrail += ρ - no += ρ - nfu -= ρ - p1 = ρ - end - - return Q, Z, reverse(μl[1:j]), no, nfu, niu -end - -""" - sklf_right!(A, E, B, C; fast = true, atol1 = 0, atol2 = 0, atol3 = 0, rtol, withQ = true, withZ = true) -> (Q, Z, νr, nc, nfuc, niuc) - -Reduce the partitioned full row rank linear pencil - - M - λN = | B | A-λE | n - m n - -with `A-λE` regular, to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Z)` using orthogonal or unitary transformation matrices `Q` and `Z` -such that `M - λN` is transformed into a Kronecker-like form `| Bt | At-λEt |` exhibiting its -right, finite and infinite Kronecker structures (also known as the generalized controllability staircase form): - - | Bc | Ac-λEc * * | - | Bt | At-λEt | = |-----|------------------------------| - | 0 | 0 Afuc-λEfuc * | - | 0 | 0 0 Aiuc-λEiuc | - -`Bt = Q'*B`, `At = Q'*A*Z`and `Et = Q'*E*Z` are returned in `B`, `A` and `E`, respectively, and `C*Z` is returned in `C` (unless `C = missing'). - -The subpencil `| Bc | Ac-λEc |` has full row rank `nc`, is in a staircase form, and contains the right Kronecker indices of `M - λN`. -The `nr`-dimensional vector `νr` contains the row and column dimensions of the blocks -of the staircase form `| Bc | Ac-λI |` such that `i`-th block has dimensions `νr[i] x νr[i-1]` (with `νr[0] = m`) and -has full row rank. The difference `νr[i-1]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks -of size `(i-1) x i`. - -The `nfuc x nfuc` subpencil `Afuc-λEfuc` contains the finite eigenvalues of `M - λN` (also called the uncontrollable finite eigenvalues of `A - λE`). - -The `niuc x niuc` subpencil `Aiuc-λEiuc` contains the infinite eigenvalues of `M - λN` (also called the uncontrollable infinite eigenvalues of `A - λE`). - -The keyword arguments `atol1`, `atol2`, , `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the -nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `B`, -and the relative tolerance for the nonzero elements of `A`, `E` and `B`. -The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting -if `fast = true` or the more reliable SVD-decompositions if `fast = false`. -The performed orthogonal or unitary left transformations are accumulated in the matrix `Q` if `withQ = true`. -Otherwise, `Q` is set to `nothing`. -The performed orthogonal or unitary right transformations are accumulated in the matrix `Z` if `withZ = true` or -`C` is provided. Otherwise, `Z` is set to `nothing`. -""" -function sklf_right!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, - atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))), atol3::Real = zero(real(eltype(B))), - rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)), - withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat - n = LinearAlgebra.checksquare(A) - n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions")) - n1, m = size(B) - n == n1 || throw(DimensionMismatch("A and B must have the same number of rows")) - (!ismissing(C) && n !== size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns")) - - maxmn = max(m,n) - νr = Vector{Int}(undef,maxmn) - nfu = 0 - niu = 0 - tol1 = atol1 - withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing - withZ ? Z = Matrix{eltype(A)}(I,n,n) : Z = nothing - - # fast returns for null dimensions - if m == 0 && n == 0 - return Q, Z, νr, n, nfu, niu - elseif n == 0 - νr[1] = 0 - return Q, Z, νr[1:1], n, nfu, niu - end - - tolA = max(atol1, rtol*opnorm(A,1)) - tolB = max(atol3, rtol*opnorm(B,1)) - - ρ1 = _sreduceB!(A, E, B, Q, tolB; fast = fast, withQ = withQ) - ρ1 == n && (return Q, Z, [ρ1], n, nfu, niu) - - n1, m1, p1 = _preduceBF!(A, E, Q, Z, missing, C; roff = ρ1, atol = atol2, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ) - - mrinf = ρ1 - nrinf = 0 - rtrail = 0 - ctrail = 0 - niu = 0 - while p1 > 0 - # Step 1 & 2: Dual algorithm PREDUCE - τ, ρ = _preduce2!(n1, m1, p1, A, E, Q, Z, tolA, missing, C; fast = fast, - roff = mrinf, coff = nrinf, rtrail = rtrail, ctrail = ctrail, withQ = withQ, withZ = withZ) - ρ+τ == p1 || error("| B | A-λE | has no full row rank") - ctrail += p1 - rtrail += p1 - niu += p1 - n1 -= ρ - p1 = ρ - m1 -= τ - end - - if ρ1 > 0 - i = 1 - νr[1] = ρ1 - else - return Q, Z, νr[1:0], 0, n1, niu - end - if m1 > 0 - imA11 = 1:n-rtrail - A11 = view(A,imA11,1:n) - E11 = view(E,imA11,1:n) - ismissing(C) ? C1 = missing : (C1 = view(C,:,1:n)) - end - nc = ρ1 - nfu = n1 - while m1 > 0 - # Step 3: Particular case of the standard algorithm PREDUCE - ρ = _preduce3!(nfu, m1, A11, E11, Q, Z, tolA, missing, C1, fast = fast, coff = nrinf, roff = mrinf, ctrail = ctrail, withQ = withQ, withZ = withZ) - ρ == 0 && break - i += 1 - νr[i] = ρ - mrinf += ρ - nrinf += m1 - nc += ρ - nfu -= ρ - m1 = ρ - end - - return Q, Z, νr[1:i], nc, nfu, niu -end -""" - sklf_right!(A, B, C; fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true) -> (Q, νr, nc, nuc) - -Reduce the partitioned full row rank linear pencil - - M - λN = | B | A-λI | n - m n - -to an equivalent form `F - λG = Q'*(M - λN)*diag(I,Q)` using orthogonal or unitary transformation matrix `Q` -such that `M - λN` is transformed into a Kronecker-like form `| Bt | At-λI |` exhibiting its -right and finite Kronecker structures (also known as the controllability staircase form): - - | Bc | Ac-λI * | - | Bt | At-λI | = |-----|----------------| - | 0 | 0 Auc-λI | - -`Bt = Q'*B` and `At = Q'*A*Q` are returned in `B` and `A`, respectively, and `C*Q` is returned in `C` (unless `C = missing'). - -The subpencil `| Bc | Ac-λI |` has full row rank `nc`, is in a staircase form, and contains the right Kronecker indices of `M - λN`. -The `nr`-dimensional vector `νr` contains the row and column dimensions of the blocks -of the staircase form `| Bc | Ac-λI |` such that `i`-th block has dimensions `νr[i] x νr[i-1]` (with `νr[0] = m`) and -has full row rank. The difference `νr[i-1]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks -of size `(i-1) x i`. - -The `nuc x nuc` matrix `Auc` contains the (finite) eigenvalues of `M - λN` (also called the uncontrollable eigenvalues of `A`). - -The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the -nonzero elements of `A`, the absolute tolerance for the nonzero elements of `B`, and the relative tolerance -for the nonzero elements of `A` and `B`. -The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting -if `fast = true` or the more reliable SVD-decompositions if `fast = false`. -The performed orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`. -Otherwise, `Q` is set to `nothing`. -""" -function sklf_right!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), - rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)), - withQ::Bool = true) where T1 <: BlasFloat - n = LinearAlgebra.checksquare(A) - n1, m = size(B) - n == n1 || throw(DimensionMismatch("A and B must have the same number of rows")) - (!ismissing(C) && n !== size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns")) - - νr = Vector{Int}(undef,max(n,1)) - nu = 0 - tol1 = atol1 - withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing - - # fast returns for null dimensions - if m == 0 && n == 0 - return Q, νr[1:0], 0, 0 - elseif n == 0 - νr[1] = 0 - return Q, νr[1:1], 0, 0 - elseif m == 0 - return Q, νr[1:0], 0, n - end - - tolA = max(atol1, rtol*opnorm(A,1)) - tolB = max(atol2, rtol*opnorm(B,1)) - - i = 0 - init = true - roff = 0 - coff = 0 - nc = 0 - nu = n - while m > 0 && nc < n - init = (i == 0) - init ? tol = tolB : tol = tolA - # Step 3: Particular case of the standard algorithm PREDUCE - ρ = _sreduceBA!(nu, m, A, B, C, Q, tol, fast = fast, init = init, roff = roff, coff = coff, withQ = withQ) - ρ == 0 && break - i += 1 - νr[i] = ρ - roff += ρ - init || (coff += m) - nc += ρ - nu -= ρ - m = ρ - end - - return Q, νr[1:i], nc, nu -end -""" - sklf_left!(A, C, B; fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true) -> (Q, μl, no, nuo) - -Reduce the partitioned full column rank linear pencil - - M - λN = | A-λI | n - | C | p - n - -to an equivalent form `F - λL = diag(Q',I)*(M - λN)*Q` using orthogonal or unitary transformation matrix `Q` -such that `M - λN` is transformed into a Kronecker-like form exhibiting its -finite and left Kronecker structures (also known as the observability staircase form): - - | Auo-λI | * | - | At-λI | | 0 | Ao-λI | - |-------| = |---------|-------| - | Ct | | 0 | Co | - -`Ct = C*Q` and `At = Q'*A*Q` are returned in `C` and `A`, respectively, and `Q'*B` is returned in `B` (unless `B = missing'). - -The subpencil `| Ao-λI; Co |` has full column rank `no`, is in a staircase form, -and contains the left Kronecker indices of `M - λN`. - -The `nl`-dimensional vector `μl` contains the row and column dimensions of the blocks -of the staircase form such that `i`-th block has dimensions `μl[nl-i] x μl[nl-i+1]` (with μl[0] = p) and -has full column rank. The difference `μl[nl-i]-μl[nl-i+1]` for `i = 1, 2, ..., nl` is the number of elementary Kronecker blocks -of size `i x (i-1)`. - -The `nuo x nuo` matrix `Auo` contains the (finite) eigenvalues of `M - λN` (also called the unobservable eigenvalues of `A`). - -The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the -nonzero elements of `A`, the absolute tolerance for the nonzero elements of `C`, and the relative tolerance -for the nonzero elements of `A` and `C`. -The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting -if `fast = true` or the more reliable SVD-decompositions if `fast = false`. -The performed orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`. -Otherwise, `Q` is set to `nothing`. -""" -function sklf_left!(A::AbstractMatrix{T1}, C::AbstractMatrix{T1}, B::Union{AbstractMatrix{T1},Missing}; fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), - rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)), - withQ::Bool = true) where T1 <: BlasFloat - n = LinearAlgebra.checksquare(A) - p, n1 = size(C) - n == n1 || throw(DimensionMismatch("A and C must have the same number of columns")) - (!ismissing(B) && n !== size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows")) - - μl = Vector{Int}(undef,max(n,p)) - nu = 0 - tol1 = atol1 - withQ ? Q = Matrix{eltype(A)}(I,n,n) : Q = nothing - - # fast returns for null dimensions - if p == 0 && n == 0 - return Q, μl[1:0], 0, 0 - elseif n == 0 - μl[1] = 0 - return Q, μl[1:1], 0, 0 - elseif p == 0 - return Q, μl[1:0], 0, n - end - - tolA = max(atol1, rtol*opnorm(A,1)) - tolC = max(atol2, rtol*opnorm(C,Inf)) - - i = 0 - init = true - rtrail = 0 - ctrail = 0 - no = 0 - nu = n - while p > 0 && no < n - init = (i == 0) - init ? tol = tolC : tol = tolA - # Step 3: Particular case of the standard algorithm PREDUCE - ρ = _sreduceAC!(nu, p, A, C, B, Q, tol, fast = fast, init = init, rtrail = rtrail, ctrail = ctrail, withQ = withQ) - ρ == 0 && break - i += 1 - μl[i] = ρ - ctrail += ρ - init || (rtrail += p) - no += ρ - nu -= ρ - p = ρ - end - - return Q, reverse(μl[1:i]), no, nu -end diff --git a/src/slputil.jl b/src/slputil.jl index 5a31f69..066c0f6 100644 --- a/src/slputil.jl +++ b/src/slputil.jl @@ -1,5 +1,83 @@ + +""" + _sreduceB!(A::AbstractMatrix{T1},E::AbstractMatrix{T1},B::AbstractMatrix{T1},Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast = true, withQ = true) + +Reduce the `n x m` matrix `B` using an orthogonal or unitary similarity transformation `Q1` to the row +compressed form + + BT = Q1'*B = [ B11 ] ρ + [ 0 ] n-ρ + m + +where `B11` has full row rank `ρ`. `Q1'*A`, `Q1'*E` and `BT` are returned in `A`, `E` and `B`, respectively. +The performed orthogonal or unitary transformations are accumulated in `Q` if `withQ = true`. +The rank decisions use the absolute tolerance `tol` for the nonzero elements of `B`. """ - _sreduceC!(A::AbstractMatrix,E::AbstractMatrix,C::AbstractMatrix,Z::Union{AbstractMatrix,Nothing}, tol; +function _sreduceB!(A::AbstractMatrix{T1},E::AbstractMatrix{T1},B::AbstractMatrix{T1}, + Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, withQ::Bool = true) where T1 <: BlasFloat + n, m = size(B) + (m == 0 || n == 0) && (return 0) + mA, nA = size(A) + T = eltype(A) + ZERO = zero(T) + if m == 1 + b = view(B,:,1) + if n == 1 + abs(b[1]) > tol ? (return 1) : (b[1] = ZERO; return 0) + else + τ, β = larfg!(b) + if abs(β) <= tol + b[:] = zeros(T,n) + return 0 + else + larf!('L', b, conj(τ), A) + larf!('L', b, conj(τ), E) + withQ && larf!('R', b, τ, Q) + b[:] = [ β; zeros(T,n-1)] + return 1 + end + end + end + if fast + _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(B) + ρ = count(x -> x > tol, abs.(diag(B))) + T <: Complex ? tran = 'C' : tran = 'T' + LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,A) + LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,E) + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',B,τ,Q) + B[:,:] = [ triu(B[1:ρ,:])[:,invperm(jpvt)]; zeros(T,n-ρ,m) ] + else + if n > m + _, τ = LinearAlgebra.LAPACK.geqrf!(B) + T <: Complex ? tran = 'C' : tran = 'T' + LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,A) + LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,E) + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',B,τ,Q) + B[:,:] = [ triu(B[1:m,:]); zeros(T,n-m,m) ] + end + mn = min(n,m) + ics = 1:mn + jcs = 1:m + SVD = svd(B[ics,jcs], full = true) + ρ = count(x -> x > tol, SVD.S) + if ρ == mn + return ρ + else + B[ics,jcs] = [ Diagonal(SVD.S[1:ρ])*SVD.Vt[1:ρ,:]; zeros(T,mn-ρ,m) ] + if ρ == 0 + return ρ + end + end + withQ && (Q[:,ics] = Q[:,ics]*SVD.U) + A[ics,:] = SVD.U'*A[ics,:] + E[ics,:] = SVD.U'*E[ics,:] + end + return ρ +end +""" + _sreduceC!(A::AbstractMatrix{T1},E::AbstractMatrix{T1},C::AbstractMatrix{T1},Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; fast = true, withZ = true) Reduce the `p x n` matrix `C` using an orthogonal or unitary similarity transformation `Z1` to the column @@ -12,8 +90,9 @@ where `C11` has full column rank `ρ`. `A*Z1`, `E*Z1` and `CT` are returned in ` The performed orthogonal or unitary transformations are accumulated in `Z` if `withZ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `C`. """ -function _sreduceC!(A::AbstractMatrix,E::AbstractMatrix,C::AbstractMatrix,Z::Union{AbstractMatrix,Nothing}, tol; - fast = true, withZ = true) +function _sreduceC!(A::AbstractMatrix{T1},E::AbstractMatrix{T1},C::AbstractMatrix{T1}, + Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, withZ::Bool = true) where T1 <: BlasFloat p, n = size(C) (p == 0 || n == 0) && (return 0) mA, nA = size(A) @@ -81,89 +160,13 @@ function _sreduceC!(A::AbstractMatrix,E::AbstractMatrix,C::AbstractMatrix,Z::Uni Z1 = reverse(SVD.V,dims=2) withZ && (Z[:,jcs] = Z[:,jcs]*Z1) A[:,jcs] = A[:,jcs]*Z1 - end - return ρ -end - -""" - _sreduceB!(A::AbstractMatrix,E::AbstractMatrix,B::AbstractMatrix,Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, withQ = true) - -Reduce the `n x m` matrix `B` using an orthogonal or unitary similarity transformation `Q1` to the row -compressed form - - BT = Q1'*B = [ B11 ] ρ - [ 0 ] n-ρ - m - -where `B11` has full row rank `ρ`. `Q1'*A`, `Q1'*E` and `BT` are returned in `A`, `E` and `B`, respectively. -The performed orthogonal or unitary transformations are accumulated in `Q` if `withQ = true`. -The rank decisions use the absolute tolerance `tol` for the nonzero elements of `B`. -""" -function _sreduceB!(A::AbstractMatrix,E::AbstractMatrix,B::AbstractMatrix,Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, withQ = true) - n, m = size(B) - (m == 0 || n == 0) && (return 0) - mA, nA = size(A) - T = eltype(A) - ZERO = zero(T) - if m == 1 - b = view(B,:,1) - if n == 1 - abs(b[1]) > tol ? (return 1) : (b[1] = ZERO; return 0) - else - τ, β = larfg!(b) - if abs(β) <= tol - b[:] = zeros(T,n) - return 0 - else - larf!('L', b, conj(τ), A) - larf!('L', b, conj(τ), E) - withQ && larf!('R', b, τ, Q) - b[:] = [ β; zeros(T,n-1)] - return 1 - end - end - end - if fast - _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(B) - ρ = count(x -> x > tol, abs.(diag(B))) - T <: Complex ? tran = 'C' : tran = 'T' - LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,A) - LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,E) - withQ && LinearAlgebra.LAPACK.ormqr!('R','N',B,τ,Q) - B[:,:] = [ triu(B[1:ρ,:])[:,invperm(jpvt)]; zeros(T,n-ρ,m) ] - else - if n > m - _, τ = LinearAlgebra.LAPACK.geqrf!(B) - T <: Complex ? tran = 'C' : tran = 'T' - LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,A) - LinearAlgebra.LAPACK.ormqr!('L',tran,B,τ,E) - withQ && LinearAlgebra.LAPACK.ormqr!('R','N',B,τ,Q) - B[:,:] = [ triu(B[1:m,:]); zeros(T,n-m,m) ] - end - mn = min(n,m) - ics = 1:mn - jcs = 1:m - SVD = svd(B[ics,jcs], full = true) - ρ = count(x -> x > tol, SVD.S) - if ρ == mn - return ρ - else - B[ics,jcs] = [ Diagonal(SVD.S[1:ρ])*SVD.Vt[1:ρ,:]; zeros(T,mn-ρ,m) ] - if ρ == 0 - return ρ - end - end - withQ && (Q[:,ics] = Q[:,ics]*SVD.U) - A[ics,:] = SVD.U'*A[ics,:] - E[ics,:] = SVD.U'*E[ics,:] + E[:,jcs] = E[:,jcs]*Z1 end return ρ end """ - _sreduceBA!(n::Int,m::Int,A::AbstractMatrix,B::AbstractMatrix,C::Union{AbstractMatrix,Missing},Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, init = true, roff = 0, coff = 0, withQ = true) + _sreduceBA!(n::Int,m::Int,A::AbstractMatrix{T1},B::AbstractMatrix{T1},C::Union{AbstractMatrix{T1},Missing},Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast = true, init = true, roff = 0, coff = 0, withQ = true) -> ρ Reduce for `init = true`, the pair `(A,B)` using an orthogonal or unitary similarity transformation on the matrices `A` and `B` of the form `H = Q1'*A*Q1`, `G = Q1'*B`, to the form @@ -194,8 +197,9 @@ where `B11` has full row rank `ρ`. `H` and `C*diag(I,Q1)` are returned in `A` a and `B` is unchanged. The performed orthogonal or unitary transformations are accumulated in `Q` if `withQ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `A`. """ -function _sreduceBA!(n::Int,m::Int,A::AbstractMatrix,B::AbstractMatrix,C::Union{AbstractMatrix,Missing},Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, init = true, roff = 0, coff = 0, withQ = true) +function _sreduceBA!(n::Int,m::Int,A::AbstractMatrix{T1},B::AbstractMatrix{T1},C::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, init::Bool = true, roff::Int = 0, coff::Int = 0, withQ::Bool = true) where T1 <: BlasFloat (m == 0 || n == 0) && (return 0) mA, nA = size(A) T = eltype(A) @@ -276,8 +280,9 @@ function _sreduceBA!(n::Int,m::Int,A::AbstractMatrix,B::AbstractMatrix,C::Union{ return ρ end """ - _sreduceAC!(n::Int,p::Int,A::AbstractMatrix,C::AbstractMatrix,B::Union{AbstractMatrix,Missing},Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, init = true, rtrail = 0, ctrail = 0, withQ = true) + _sreduceAC!(n::Int,p::Int,A::AbstractMatrix{T1},C::AbstractMatrix{T1},B::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast = true, init = true, rtrail = 0, ctrail = 0, withQ = true) -> ρ Reduce for `init = true`, the pair `(A,C)` using an orthogonal or unitary similarity transformation on the matrices `A` and `C` of the form `H = Q1'*A*Q1`, `L = C*Q1`, to the form @@ -310,8 +315,9 @@ where `C11` has full column rank `ρ`. `H` and `diag(Q1',I)*B` are returned in ` and `C` is unchanged. The performed orthogonal or unitary transformations are accumulated in `Q` if `withQ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `A`. """ -function _sreduceAC!(n::Int,p::Int,A::AbstractMatrix,C::AbstractMatrix,B::Union{AbstractMatrix,Missing},Q::Union{AbstractMatrix,Nothing}, tol; - fast = true, init = true, rtrail = 0, ctrail = 0, withQ = true) +function _sreduceAC!(n::Int,p::Int,A::AbstractMatrix{T1},C::AbstractMatrix{T1},B::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, init::Bool = true, rtrail::Int = 0, ctrail::Int = 0, withQ::Bool = true) where T1 <: BlasFloat (p == 0 || n == 0) && (return 0) mA, nA = size(A) T = eltype(A) @@ -397,3 +403,339 @@ function _sreduceAC!(n::Int,p::Int,A::AbstractMatrix,C::AbstractMatrix,B::Union{ end return ρ end +""" + _sreduceBAE!(n::Int,m::Int,A::AbstractMatrix{T1},E::AbstractMatrix{T1},B::AbstractMatrix{T1},C::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast = true, init = true, roff = 0, coff = 0, withQ = true) + +Reduce for `init = true`, the pair `(A-λE,B)`, with E upper-triangular, using an orthogonal or unitary +similarity transformations on the matrices `A`, `E` and `B` of the form `At = Q1'*A*Z1`, `Et = Q1'*E*Z1`, +`Bt = Q1'*B`, to the form + + Bt = [ B11 ] At = [ A11 A12 ] ρ Et = [ E11 E12 ] ρ + [ 0 ] B2 A22 ] n-ρ [ 0 E22 ] n-ρ + m ρ n-ρ ρ n-ρ + +where `B11` has full row rank `ρ` and `Et` is upper-triangular. `Bt`, `At`, `Et` and `C*Z1` are returned in +`B`, `A`, `E` and `C`, respectively. +The performed orthogonal or unitary transformations are accumulated in `Q` as `Q <- Q*Q1` if `withQ = true` +and in `Z` as `Z <- Z*Z1` if `withZ = true`. +The rank decisions use the absolute tolerance `tol` for the nonzero elements of `B`. + +Reduce for `init = false`, the matrices `A` and `E` of the form + + A = [ * * * ] roff E = [ * * * ] roff + [ 0 B1 A1 ] n [ 0 0 E1 ] n + coff m n coff m n + +with 'E1' upper triangular, using an orthogonal or unitary similarity transformations on the submatrices +'A1', 'E1' and 'B1' of the form 'At1 = Q1'*A1*Z1', 'Et1 = Q1'*E1*Z1', 'Bt1 = Q1'*B1', to the form + + [ * * * * ] roff + At = diag(I,Q1')*A*diag(I,Z1) = [ 0 B11 A11 A12 ] ρ + [ 0 0 B2 A22 ] n-ρ + coff m ρ n-ρ + + [ * * * * ] roff + Et = diag(I,Q1')*A*diag(I,Z1) = [ 0 0 E11 E12 ] ρ + [ 0 0 0 E22 ] n-ρ + coff m ρ n-ρ + + +where `B11` has full row rank `ρ`, and `E11` and `E22` are upper triangular. `At`, `Et` and `C*diag(I,Z1)` +are returned in `A`, `E` and `C`, respectively, and `B` is unchanged. +The performed orthogonal or unitary transformations are accumulated in `Q` as `Q <- Q*diag(I,Q1)` if `withQ = true` +and in `Z` as `Z <- Z*diag(I,Z1)` if `withZ = true`. +The rank decisions use the absolute tolerance `tol` for the nonzero elements of `A`. +""" +function _sreduceBAE!(n::Int,m::Int,A::AbstractMatrix{T1},E::AbstractMatrix{T1},B::AbstractMatrix{T1},C::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, init::Bool = true, roff::Int = 0, coff::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + (m == 0 || n == 0) && (return 0) + mA, nA = size(A) + T = eltype(A) + ZERO = zero(T) + ib = roff+1:roff+n + ia = 1:roff+n + if init + # coff and roff must be zero + (coff == 0 && roff == 0) || error("coff and roff must be zero at first call") + ja = 1:n + B1 = view(B,ib,1:m) + A1 = view(A,ib,ja) + E1 = view(E,ib,ja) + else + (coff+m+n == nA && roff+n == mA) || error("coff and roff must have compatible values with the dimensions of A") + ja = coff+m+1:coff+m+n + B1 = view(A,ib,coff+1:coff+m) + A1 = view(A,ib,ja) + E1 = view(E,ib,ja) + end + if fast + ρ = 0 + nrm = similar(real(A),m) + jp = Vector(1:m) + nm = min(n,m) + for j = 1:nm + for l = j:m + nrm[l] = norm(B1[j:n,l]) + end + nrmax, ind = findmax(nrm[j:m]) + ind += j-1 + if nrmax < tol + break + else + ρ += 1 + end + if ind !== j + (jp[j], jp[ind]) = (jp[ind], jp[j]) + (B1[:,j],B1[:,ind]) = (B1[:,ind],B1[:,j]) + end + for ii = n:-1:j+1 + iim1 = ii-1 + if B1[ii,j] !== ZERO + G, B1[iim1,j] = givens(B1[iim1,j],B1[ii,j],iim1,ii) + B1[ii,j] = ZERO + lmul!(G,view(B1,:,j+1:m)) + lmul!(G,A1) + lmul!(G,view(E1,:,iim1:n)) + withQ && rmul!(view(Q,:,ib),G') + G, r = givens(conj(E1[ii,ii]),conj(E1[ii,iim1]),ii,iim1) + E1[ii,ii] = conj(r) + E1[ii,iim1] = ZERO + rmul!(view(E,1:roff+iim1,ja),G') + withZ && rmul!(view(Z,:,ja),G') + rmul!(view(A,:,ja),G') + ismissing(C) || rmul!(view(C,:,ja),G') + end + end + end + B1[:,:] = [ B1[1:ρ,invperm(jp)]; zeros(T,n-ρ,m) ] + return ρ + else + if n > m + for j = 1:m + for ii = n:-1:j+1 + iim1 = ii-1 + if B1[ii,j] !== ZERO + G, B1[iim1,j] = givens(B1[iim1,j],B1[ii,j],iim1,ii) + B1[ii,j] = ZERO + lmul!(G,view(B1,:,j+1:m)) + lmul!(G,A1) + lmul!(G,view(E1,:,iim1:n)) + withQ && rmul!(view(Q,:,ib),G') + G, r = givens(conj(E1[ii,ii]),conj(E1[ii,iim1]),ii,iim1) + E1[ii,ii] = conj(r) + E1[ii,iim1] = ZERO + rmul!(view(E,1:roff+iim1,ja),G') + withZ && rmul!(view(Z,:,ja),G') + rmul!(view(A,:,ja),G') + ismissing(C) || rmul!(view(C,:,ja),G') + end + end + end + end + mn = min(n,m) + if mn > 0 + ics = 1:mn + jcs = 1:m + SVD = svd(B1[ics,jcs], full = true) + ρ = count(x -> x > tol, SVD.S) + if ρ == mn + return ρ + else + B1[ics,jcs] = [ Diagonal(SVD.S[1:ρ])*SVD.Vt[1:ρ,:]; zeros(T,mn-ρ,m) ] + if ρ == 0 + return ρ + end + end + ibt = roff+1:roff+mn + jt = coff+m+1:nA + withQ && (Q[:,ibt] = Q[:,ibt]*SVD.U) + E[ibt,jt] = SVD.U'*E[ibt,jt] + A[ibt,jt] = SVD.U'*A[ibt,jt] + tau = similar(E,mn) + jt1 = coff+m+1:coff+m+mn + E11 = view(E,ibt,jt1) + LinearAlgebra.LAPACK.gerqf!(E11,tau) + eltype(A) <: Complex ? tran = 'C' : tran = 'T' + LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(A,:,jt1)) + withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(Z,:,jt1)) + LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(E,1:roff,jt1)) + ismissing(C) || LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(C,:,jt1)) + triu!(E11) + else + ρ = 0 + end + end + return ρ +end +""" + _sreduceAEC!(n::Int,p::Int,A::AbstractMatrix{T1},E::AbstractMatrix{T1},C::AbstractMatrix{T1},B::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast = true, init = true, rtrail = 0, ctrail = 0, withQ = true, withZ = true) -> ρ + +Reduce for `init = true`, the pair `(A-λE,C)`, with E upper-triangular, using an orthogonal or unitary +similarity transformations on the matrices `A`, `E` and `C` of the form `At = Q1'*A*Z1`, `Et = Q1'*E*Z1`, +`Ct = C*Z1`, to the form + + Ct = [ 0 C11 ] p At = [ A11 A12 ] n-ρ Et = [ E11 E12 ] n-ρ + n-ρ ρ [ C2 A22 ] ρ [ 0 E22 ] ρ + n-ρ ρ n-ρ ρ + +where `C11` has full column rank `ρ` and `Et` is upper-triangular. `Ct`, `At`, `Et` and `Q1'*B` are returned in +`C`, `A`, `E` and `B`, respectively. +The performed orthogonal or unitary transformations are accumulated in `Q` as `Q <- Q*Q1` if `withQ = true` +and in `Z` as `Z <- Z*Z1` if `withZ = true`. +The rank decisions use the absolute tolerance `tol` for the nonzero elements of `C`. + +Reduce for `init = false`, the matrices `A` and `E` of the form + + [ A1 * ] n [ E1 * ] n + A = [ C1 * ] p E = [ 0 * ] p + [ 0 * ] rtrail [ 0 * ] rtrail + n ctrail + +with 'E1' upper triangular, using an orthogonal or unitary similarity transformations on the submatrices +'A1', 'E1' and 'C1' of the form 'At1 = Q1'*A1*Z1', 'Et1 = Q1'*E1*Z1', 'Ct1 = C1*Z1', to the form + + + [ A11 A12 * ] n-ρ + [ C2 A22 * ] ρ + At = diag(Q1',I)*A*diag(Q1,I) = [ 0 C11 * ] p + [ 0 0 * ] rtrail + n-ρ ρ ctrail + + [ E11 E12 * ] n-ρ + [ 0 E22 * ] ρ + Et = diag(Q1',I)*A*diag(Q1,I) = [ 0 0 * ] p + [ 0 0 * ] rtrail + n-ρ ρ ctrail + +where `C11` has full column rank `ρ`, and `E11` and `E22` are upper triangular. `At`, `Et` and `diag(I,Q1')*B` +are returned in `A`, `E` and `B`, respectively, and `C` is unchanged. +The performed orthogonal or unitary transformations are accumulated in `Q` as `Q <- Q*diag(I,Q1)` if `withQ = true` +and in `Z` as `Z <- Z*diag(I,Z1)` if `withZ = true`. +The rank decisions use the absolute tolerance `tol` for the nonzero elements of `A`. +""" +function _sreduceAEC!(n::Int,p::Int,A::AbstractMatrix{T1},E::AbstractMatrix{T1},C::AbstractMatrix{T1},B::Union{AbstractMatrix{T1},Missing}, + Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; + fast::Bool = true, init::Bool = true, rtrail::Int = 0, ctrail::Int = 0, + withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat + (p == 0 || n == 0) && (return 0) + mA, nA = size(A) + T = eltype(A) + ZERO = zero(T) + ia = 1:n + ja = 1:nA + A1 = view(A,ia,ia) + E1 = view(E,ia,ia) + if init + # coff and roff must be zero + (ctrail == 0 && rtrail == 0) || error("rtrail and ctrail must be zero at first call") + C1 = view(C,1:p,ia) + else + (ctrail+n == nA && rtrail+n+p == mA) || error("rtrail and ctrail must have compatible values with the dimensions of A") + C1 = view(A,n+1:n+p,ia) + end + if fast + ρ = 0 + nrm = similar(real(A),p) + jp = Vector(1:p) + np = min(n,p) + for i = 1:np + ii = p-i+1 + for l = 1:ii + nrm[l] = norm(C1[l,1:n-i+1]) + end + nrmax, ind = findmax(nrm[1:ii]) + if nrmax < tol + break + else + ρ += 1 + end + if ind !== ii + (jp[ii], jp[ind]) = (jp[ind], jp[ii]) + (C1[ii,:],C1[ind,:]) = (C1[ind,:],C1[ii,:]) + end + for jj = 1:n-i + jjp1 = jj+1 + if C1[ii,jj] !== ZERO + G, r = givens(conj(C1[ii,jjp1]),conj(C1[ii,jj]),jjp1,jj) + C1[ii,jjp1] = conj(r) + C1[ii,jj] = ZERO + rmul!(view(C1,1:ii-1,:),G') + rmul!(A1,G') + rmul!(view(E1,1:jjp1,ia),G') + withZ && rmul!(view(Z,:,ia),G') + G, E[jj,jj] = givens(E[jj,jj],E[jjp1,jj],jj,jjp1) + E[jjp1,jj] = ZERO + lmul!(G,view(E,ia,jjp1:nA)) + withQ && rmul!(view(Q,:,ia),G') + lmul!(G,view(A,ia,ja)) + ismissing(B) || lmul!(G,view(B,ia,:)) + end + end + end + C1[:,1:n] = [ zeros(T,p,n-ρ) C1[invperm(jp),n-ρ+1:n]]; + else + if n > p + for i = 1:p + ii = p-i+1 + for jj = 1:n-i + jjp1 = jj+1 + if C1[ii,jj] !== ZERO + G, r = givens(conj(C1[ii,jjp1]),conj(C1[ii,jj]),jjp1,jj) + C1[ii,jjp1] = conj(r) + C1[ii,jj] = ZERO + rmul!(view(C1,1:ii-1,:),G') + rmul!(A1,G') + rmul!(view(E,1:jjp1,ia),G') + withZ && rmul!(view(Z,:,ia),G') + G, E[jj,jj] = givens(E[jj,jj],E[jjp1,jj],jj,jjp1) + E[jjp1,jj] = ZERO + lmul!(G,view(E,ia,jjp1:nA)) + withQ && rmul!(view(Q,:,ia),G') + lmul!(G,view(A,ia,ja)) + ismissing(B) || lmul!(G,view(B,ia,:)) + end + end + end + end + pn = min(n,p) + if pn > 0 + ics = 1:p + jcs = n-pn+1:n + SVD = svd(C1[ics,jcs], full = true) + ρ = count(x -> x > tol, SVD.S) + if ρ == pn + return ρ + else + Q1 = reverse(SVD.U,dims=2) + C1[ics,jcs] = [ zeros(T,p,pn-ρ) Q1[:,p-ρ+1:end]*Diagonal(reverse(SVD.S[1:ρ])) ] + if ρ == 0 + return ρ + end + end + Z1 = reverse(SVD.V,dims=2) + jt = n-pn+1:n + withZ && (Z[:,jt] = Z[:,jt]*Z1) + A[ia,jt] = A[ia,jt]*Z1 + E[ia,jt] = E[ia,jt]*Z1 # more efficient computation possible + jt1 = n+1:nA + tau = similar(E,pn) + E22 = view(E,jt,jt) + LinearAlgebra.LAPACK.geqrf!(E22,tau) + eltype(A) <: Complex ? tran = 'C' : tran = 'T' + LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(A,jt,ja)) + withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E22,tau,view(Q,:,jt)) + LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(E,jt,jt1)) + ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(B,jt,:)) + triu!(E22) + else + ρ = 0 + end + end + return ρ +end diff --git a/test/runtests.jl b/test/runtests.jl index 5a94ba1..510118f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using Test, MatrixPencils include("test_regular.jl") include("test_sklf.jl") include("test_sklfapps.jl") + include("test_lstools.jl") end end diff --git a/test/test_lstools.jl b/test/test_lstools.jl new file mode 100644 index 0000000..2c48619 --- /dev/null +++ b/test/test_lstools.jl @@ -0,0 +1,535 @@ +module Test_linsystools + +using LinearAlgebra +using MatrixPencils +using Test + + +@testset "Linear System Tools" begin + + +@testset "lsminreal and lsminreal2" begin + +A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = zeros(0,0); D2 = zeros(0,0); +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); + +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D) + +@test lsequal(A,E,B,C,D,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == 0 + +A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = zeros(0,0); D2 = zeros(0,0); +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); + +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D) + +@test lsequal(A,E,B,C,D,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == 0 + + + +# Example 1: DeTeran, Dopico, Mackey, ELA 2009 +# P = zeros(2,2,3); +# P[:,:,1] = [0 0; 0 1.]; +# P[:,:,2] = [0 1.; 1. 0]; +# P[:,:,3] = [1. 0; 0 0]; + +# # build a strong (least order) structured linearization which preserves the finite eigenvalues, +# # infinite zeros (with multiplicities), the left and right Kronecker structure +# A2,E2,B2,C2,D2 = lsbuild(P) + +A2 = [1.0 0.0 0.0 0.0 0.0 0.0 +0.0 1.0 0.0 0.0 0.0 0.0 +0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0]; + +E2 = [0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0 +0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0]; + +B2 = [0.0 0.0 +0.0 0.0 +0.0 1.0 +1.0 0.0 +1.0 0.0 +0.0 0.0]; + +C2 = [-1.0 0.0 0.0 0.0 0.0 0.0 +0.0 -1.0 0.0 0.0 0.0 0.0]; + +D2 = [0.0 0.0 +0.0 1.0]; + + +# A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +# Q, Z, νr, nc, nuc = sklf_rightfin!(A,E,B,C,fast=false) +# @test Q'*A2*Z≈A && Q'*E2*Z≈E && Q'*B2≈B && C2*Z≈C + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = + lsminreal2(A,E,B,C,D) +@test lsequal(A,E,B,C,D,A1,E1,B1,C1,D1) && + nuc == 2 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ Float64[] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([1], [1], [1, 1], 0) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([1], [1], [1, 1], 0) + +#A2,E2,B2,C2,D2 = lsbuild(P) +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D) +@test lsequal(A,E,B,C,D,A1,E1,B1,C1,D1) && + nuc == 2 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ Float64[] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([1], [1], [1, 1], 0) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([1], [1], [1, 1], 0) + +# Example Van Dooren & Dewilde, LAA 1983. +P = zeros(3,3,3) +P[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0] +P[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2] +P[:,:,3] = [1 4 2; 0 0 0; 1 4 2] + +# observable realization with A2 = I +A2 = [ 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]; + +E2 = [ 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]; + +B2 = [ 0.0 0.0 0.0 +0.0 0.0 0.0 +0.0 0.0 0.0 +1.0 3.0 0.0 +1.0 4.0 2.0 +0.0 -1.0 -2.0 +1.0 4.0 2.0 +0.0 0.0 0.0 +1.0 4.0 2.0]; + +C2 = [ -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0]; + +D2 = [ 1.0 2.0 -2.0 +0.0 -1.0 -2.0 +0.0 0.0 0.0]; + +# build a strong (least order) structured linearization which preserves the finite eigenvalues, +# infinite zeros (with multiplicities), the left and right Kronecker structure +#A2,E2,B2,C2,D2 = lsbuild(P) +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ [1.0] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [1., Inf, Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,obs=false); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ [1.0] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [1., Inf, Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) + + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ [1.0] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [1., Inf, Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,obs=false,finite=false); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 1 +@time val, iz, info = spzeros(A1,E1,B1,C1,D1) +@test val ≈ [1.0] && iz == [] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) +@time val, info = speigvals(A1,E1,B1,C1,D1) +@test val ≈ [1., Inf, Inf, Inf] && (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1) + +# perform lsminreal with A and E interchanged +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time E1, A1, B1, C1, D1, nuc, nuo, nse = lsminreal(E,A,B,C,D,obs = false) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 0 + +# perform lsminreal2 with A and E interchanged +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time E1, A1, B1, C1, D1, nuc, nuo, nse = lsminreal2(E,A,B,C,D,obs = false) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 5 && nuo == 0 && nse == 0 + + +# controllable realization with A2 = I +A2 = [ 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]; + +E2 = [ 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]; + +B2 = [ 0.0 0.0 0.0 +0.0 0.0 0.0 +0.0 0.0 0.0 +0.0 0.0 0.0 +0.0 0.0 0.0 +0.0 0.0 0.0 +-1.0 0.0 0.0 +0.0 -1.0 0.0 +0.0 0.0 -1.0]; + +C2 = [ 1.0 4.0 2.0 1.0 3.0 0.0 0.0 0.0 0.0 +0.0 0.0 0.0 1.0 4.0 2.0 0.0 0.0 0.0 +1.0 4.0 2.0 0.0 -1.0 -2.0 0.0 0.0 0.0 +]; + +D2 = [ 1.0 2.0 -2.0 +0.0 -1.0 -2.0 +0.0 0.0 0.0]; + + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,contr = false); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,contr = false,finite=false); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + + +A2 = [ + 1. 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 1]; +E2 = [ + 0. 0 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0]; +B2 = [ + -1. 0 0 + 0 0 0 + 0 0 0 + 0 -1 0 + 0 0 0 + 0 0 0 + 0 0 -1 + 0 0 0 + 0 0 0]; +C2 = [ + 0. 1 1 0 3 4 0 0 2 + 0 1 0 0 4 0 0 2 0 + 0 0 1 0 -1 4 0 -2 2]; + +D2 = zeros(Float64,3,3); + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = + lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = + lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,contr=false); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,I,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,I,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 6 && nuo == 3 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = + lsminreal2(A,E,B,C,D,contr=false,finite=false,atol1 = 1.e-7, atol2 = 1.e-7); +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 5 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,I,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,I,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 6 && nuo == 3 && nse == 0 + + +A2 = [ + -2 -3 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 + 0 0 -2 -3 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 1]; + E2 = [ + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0]; + B2 = [ + 1 0 + 0 0 + 0 1 + 0 0 + -1 0 + 0 0 + 0 -1 + 0 0 + 0 0]; + C2 = [ + 1 0 1 -3 0 1 0 2 0 + 0 1 1 3 0 1 0 0 1]; +D2 = zeros(Int,2,2); + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,noseig=true) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 2 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,obs=false,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,contr=false,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 2 && nse == 0 + + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,noseig=true) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 2 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,obs=false,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,contr=false,finite=false,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 2 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 2 && nuo == 0 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 2 && nuo == 0 && nse == 0 + +# Example 1 - (Varga, Kybernetika, 1990) +A2 = [ +1 0 0 0 -1 0 0 0 +0 1 0 0 0 -1 0 0 +0 0 1 0 0 0 0 0 +0 0 0 1 0 0 0 0 +0 0 0 0 -1 0 0 0 +0 0 0 0 0 -1 0 0 +0 0 0 0 3 0 1 0 +0 0 0 0 0 2 0 1 +] +E2 = [ +0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 +0 1 0 0 0 0 0 0 +0 0 1 0 0 0 0 0 +0 0 0 1 0 0 0 0 +] +B2 = [ + -1 1 + 0 0 + 0 0 + 0 0 + 1 -2 + -2 3 + 0 0 + 3 -3 +] +C2 = [ + 0 0 0 0 0 0 -1 0 + 0 0 0 0 0 0 0 -1 +] +D2 = zeros(Int,2,2); + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 6 && nuo == 0 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 6 && nuo == 0 && nse == 1 + + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 2 && nuo == 4 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 2 && nuo == 4 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A',E',C',B',D',contr=false, atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 0 && nuo == 6 && nse == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A',E',C',B',D',contr=false, atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') && + nuc == 0 && nuo == 6 && nse == 1 + +A2 = [ + 0 0 0 -24 0 0 0 0 0 0 0 + 1 0 0 -50 0 0 0 0 0 0 0 + 0 1 0 -35 0 0 0 0 0 0 0 + 0 0 1 -10 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 -30 0 0 0 + 0 0 0 0 1 0 0 -61 0 0 0 + 0 0 0 0 0 1 0 -41 0 0 0 + 0 0 0 0 0 0 1 -11 0 0 0 + 0 0 0 0 0 0 0 0 0 0 -15 + 0 0 0 0 0 0 0 0 1 0 -23 + 0 0 0 0 0 0 0 0 0 1 -9 +] +E2 = I; +B2 = reshape([18; 42; 30; 6; 10;17;8;1;0;-10;-2;],11,1) +C2 = reshape([0 0 0 0 0 0 0 1 0 0 0],1,11) +D2 = zeros(Int,1,1) + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 7 && nuo == 3 && nse == 0 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 7 && nuo == 3 && nse == 0 + +fast = true; Ty = Float64 +for fast in (true, false) + +for Ty in (Float64, Complex{Float64}) + +#fast = true; Ty = Complex{Float64} +n = 10; m = 5; p = 6; +rE = 5; +rA22 = min(2,n-rE) +n2 = rE+rA22 +n3 = n-n2 +A2 = [rand(Ty,rE,n); +rand(Ty,rA22,n2) zeros(Ty,rA22,n3); +rand(Ty,n3,rE) zeros(Ty,n3,n-rE) ]; +E2 = [rand(Ty,rE,rE) zeros(Ty,rE,n-rE); +zeros(Ty,n-rE,n)] +Q = qr(rand(Ty,n,n)).Q +Z = qr(rand(Ty,n,n)).Q +A2 = Q*A2*Z; E2 = Q*E2*Z; +B2 = rand(Ty,n,m); C2 = rand(Ty,p,n); +D2 = zeros(Ty,p,m) + + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,fast=fast) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == rA22 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2); +@time A1, E1, B1, C1, D1, nuc, nuo, nse = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,fast=fast) +@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) && + nuc == 0 && nuo == 0 && nse == rA22 + +end +end + + +end # regular testset +end # module +end # module \ No newline at end of file diff --git a/test/test_regular.jl b/test/test_regular.jl index 7e16f2b..d17fc73 100644 --- a/test/test_regular.jl +++ b/test/test_regular.jl @@ -7,6 +7,257 @@ using Test @testset "Regular Matrix Pencils Utilities" begin +@testset "_svdlikeAE" begin +A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing; +A = copy(A2); E = copy(E2); C = copy(C2); B = missing; +Q = nothing +Z = nothing + + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test rE == 0 && rA22 == 0 + +A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,2); C2 = missing; +A = copy(A2); E = copy(E2); B = copy(B2); C = missing; + +Q = nothing +Z = nothing + + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test rE == 0 && rA22 == 0 + + +n = 9; +A2 = [ + 1. 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 1]; +E2 = [ + 0. 0 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0]; +B2 = [ + -1. 0 0 + 0 0 0 + 0 0 0 + 0 -1 0 + 0 0 0 + 0 0 0 + 0 0 -1 + 0 0 0 + 0 0 0]; +C2 = [ + 0. 1 1 0 3 4 0 0 2 + 0 1 0 0 4 0 0 2 0 + 0 0 1 0 -1 4 0 -2 2]; + +A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); +T = Float64; +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 7 && rA22 == 0 + +n = 9; + +A2 = [ + 1. 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 1]; +E2 = [ + 0. 0 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0]; +B2 = [ + -1. 0 0 + 0 0 0 + 0 0 0 + 0 -1 0 + 0 0 0 + 0 0 0 + 0 0 -1 + 0 0 0 + 0 0 0]; +C2 = [ + 0. 1 1 0 3 4 0 0 2 + 0 1 0 0 4 0 0 2 0 + 0 0 1 0 -1 4 0 -2 2]; + +T = Float64; +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); + +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 7 && rA22 == 0 + +# TG01ED EXAMPLE PROGRAM DATA +n = 4; +A2 = [ + -1. 0 0 3 + 0 0 1 2 + 1 1 0 4 + 0 0 0 0]; + E2 = [ + 1. 2 0 0 + 0 1 0 1 + 3 9 6 3 + 0 0 2 0]; + B2 = [ + 1. 0 + 0 0 + 0 1 + 1 1]; + C2 = [ + -1. 0 1 0 + 0 1 -1 1]; + +T = Float64; +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); + +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 3 && rA22 == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); + +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C, fast = false) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 3 && rA22 == 1 + +n = 4; +A2 = [ + -3.330669073875470e-16 0 0 0 + -5.336424998587595e-17 4.999999999999997e-01 0 0 + 7.882634225314346e-01 4.926646390821473e-02 1.017643702629568e+00 0 + -7.804970910378840e-02 -5.073231091746246e-01 -3.602232247129060e-01 1.906644727694980e+00 + ]; + E2 = [ + 9.708552809960075e-34 -4.715554847257422e-19 7.905694150420953e-01 0 + -1.002546820843656e-17 -9.625586024907403e-20 9.682458365518573e-02 4.743416490252569e-01 + -9.878387353077869e-19 2.181405497860925e-19 -3.720759787739367e-01 4.673827146373172e-02 + 1.017230340487078e-17 7.512741891975197e-20 -6.045704470706158e-02 -4.812889603890242e-01 + ]; + B2 = [ + -1.171404622872256e-16 -6.123724356957946e-01 -2.056086568456411e-17 2.775557561562891e-16 -2.844390262617513e-16 + -5.000000000000002e-01 1.250000000000000e-01 -5.000000000000003e-01 4.139921271247625e-17 4.999999999999999e-01 + -4.926646390821468e-02 -4.803480231050927e-01 -4.926646390821469e-02 7.882634225314342e-01 4.926646390821479e-02 + -4.975668955366508e-01 -7.804970910378836e-02 -4.975668955366509e-01 -7.804970910378847e-02 -5.073231091746250e-01 + ]; + C2 = [ + 9.999999999999998e-01 0 0 0 + -9.799563198300907e-17 9.999999999999998e-01 0 0 + ]; + +T = Float64; +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); + +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 2 && rA22 == 1 + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); + +Q = Matrix{T}(I,n,n) +Z = Matrix{T}(I,n,n) + +@time rE, rA22 = _svdlikeAE!(A, E, Q, Z, B, C, fast = false) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE == 2 && rA22 == 1 + +fast = true; Ty = Float64 +for fast in (true, false) + +for Ty in (Float64, Complex{Float64}) + + +n = 10; m = 3; p = 2; +rE = 5; +rA22 = min(2,n-rE) +n2 = rE+rA22 +n3 = n-n2 +A2 = [rand(Ty,rE,n); +rand(Ty,rA22,n2) zeros(Ty,rA22,n3); +rand(Ty,n3,rE) zeros(Ty,n3,n-rE) ]; +E2 = [rand(Ty,rE,rE) zeros(Ty,rE,n-rE); +zeros(Ty,n-rE,n)] +Q = qr(rand(Ty,n,n)).Q +Z = qr(rand(Ty,n,n)).Q +A2 = Q*A2*Z; E2 = Q*E2*Z; +B2 = rand(Ty,n,m); C2 = rand(Ty,p,n); + +A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); +Q = Matrix{Ty}(I,n,n) +Z = Matrix{Ty}(I,n,n) + +@time rE1, rA221 = _svdlikeAE!(A, E, Q, Z, B, C, fast = fast, atol1=1.e-7, atol2 = 1.e-7) +@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) && + rE1 == rE && rA221 == rA22 + +end +end + +end + @testset "isregular" begin @@ -89,6 +340,7 @@ end end # isregular + @testset "fisplit" begin A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing; A = copy(A2); E = copy(E2); C = copy(C2); B = missing; @@ -110,6 +362,59 @@ A = copy(A2); E = copy(E2); B = copy(B2); C = missing; (ismissing(C) || norm(C2*Z-C1) < sqrt(eps(1.))) && ν == [] && nf == 0 && ni == 0 + +A2 = [ + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 1]; +E2 = [ + 0 0 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 0]; +B2 = [ + -1 0 0 + 0 0 0 + 0 0 0 + 0 -1 0 + 0 0 0 + 0 0 0 + 0 0 -1 + 0 0 0 + 0 0 0]; +C2 = [ + 0 1 1 0 3 4 0 0 2 + 0 1 0 0 4 0 0 2 0 + 0 0 1 0 -1 4 0 -2 2]; + +A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); +@time A1, E1, B1, C1, Q, Z, ν, nf, ni = fisplit(A,E,B,C, finite_infinite=false, atol1 = 1.e-7, atol2 = 1.e-7) +@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C1) < sqrt(eps(1.))) && + ν == [2, 2, 2, 1, 1, 1] && nf == 0 && ni == 9 + +A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); +@time A1, E1, B1, C1, Q, Z, ν, nf, ni = fisplit(A,E,B,C, finite_infinite=true, atol1 = 1.e-7, atol2 = 1.e-7) +@test norm(Q'*A2*Z-A1) < sqrt(eps(1.)) && + norm(Q'*E2*Z-E1) < sqrt(eps(1.)) && + (ismissing(B) || norm(Q'*B2-B1) < sqrt(eps(1.))) && + (ismissing(C) || norm(C2*Z-C1) < sqrt(eps(1.))) && + ν == [1, 1, 1, 2, 2, 2] && nf == 0 && ni == 9 + + fast = true; finite_infinite = true; Ty = Float64 for finite_infinite in (true, false)