From c241e1c37318d51e98db396f3b1a869e64f737a1 Mon Sep 17 00:00:00 2001 From: olof3 Date: Tue, 5 May 2020 14:06:22 +0200 Subject: [PATCH 01/13] Sylvester and Lypunov equations a la Bartels-Stewart. --- src/LyapTest.jl | 384 ++++++++++++++++++++++++++++++++++ test/test_matrix_equations.jl | 119 +++++++++++ 2 files changed, 503 insertions(+) create mode 100644 src/LyapTest.jl create mode 100644 test/test_matrix_equations.jl diff --git a/src/LyapTest.jl b/src/LyapTest.jl new file mode 100644 index 000000000..690bc1f45 --- /dev/null +++ b/src/LyapTest.jl @@ -0,0 +1,384 @@ +module LyapTest + +using LinearAlgebra +using StaticArrays + +#= +* Symmetry is not exploited when solving diagonal entries +* Perhaps it would be better with separate methods for the real and + complex versions of the block solve algorithms +* What convention for signs and tranposes? +=# + +""" + `b, d, nblocks = _schurstructure(R::AbstractMatrix)` + + Return the block strucutre of a quasi-traingular Schur matrix `R`. + + `d` contains the block sizees of each diagonal block (`1` or `2`) + + `b` contains the indices of the blocks + + `nblocks` is the number of blocks + +""" +function _schurstructure(R::AbstractMatrix, ul=Val(:u)::Union{Val{:u}, Val{:l}}) + n = size(R,1) + + d = Vector{Int}(undef, n) # block sizes + b = Vector{UnitRange{Int64}}(undef, n) # block indices + + j = 1 # column if ul=:u, row if ul=:l + k = 0 # block number + while j <= n + k += 1 + if j == n + d[k] = 1 + else + if ul == Val(:u) + d[k] = iszero(R[j+1, j]) ? 1 : 2 + else + d[k] = iszero(R[j, j+1]) ? 1 : 2 + end + end + b[k] = j:j+d[k]-1 + j += d[k] + end + resize!(d, k) + resize!(b, k) + return d, b, k +end + +# A step range representing -Inf:1:Inf +# for efficient reuse of blocked real algorithms for scalar complex version +struct FullStepRange end +Base.getindex(::FullStepRange, k::Int) = k + +# If the C entry is a 1x1 (0-dimenional view), then +# vector vector multiplication is interpreted as tranpose(A)*B +muladdrc!(C, A, B, α, β) = mul!(C, A, B, α, β) +function muladdrc!(c::Base.SubArray{<:Any,0}, A::Base.SubArray{<:Any,1}, B::Base.SubArray{<:Any,1}, α, β) + c[1] *= β + if eltype(A) <: ComplexF64 && eltype(B) <: ComplexF64 + c[1] += α*LinearAlgebra.BLAS.dotu(A, B) + else + c[1] += α*sum(A[k]*B[k] for k=1:length(A)) + end +end +function muladdrc!(c::Base.SubArray{<:Any,0}, A::Base.SubArray{<:Any,0}, B, α, β) + c[1] *= β + c[1] += α*A[1]*B[1] +end + + + +""" + + _sylvc!(A, B, C) + + Solve the continuous-time sylvester equation + + `AX + XB = C` + + for small matrices (1x1, 1x2, 2x1, 2x2) +""" +function _sylvc!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T <: Number} # Should really mostly be used for Reals + Cv = C[:] # vectorization of C + m, n = size(C) + Xv = (kron(I(n), A) + kron(transpose(B), I(m))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) + + if any(isinf, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end + + X .= SMatrix{size(C,1),size(C,2)}(reshape(Xv, size(C,1), size(C,2))) +end +function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + if size(C) == (1,1) + _sylvc!(C, SMatrix{1,1}(A), SMatrix{1,1}(B), SMatrix{1,1}(C)) + elseif size(C) == (1,2) + _sylvc!(C, SMatrix{1,1}(A), SMatrix{2,2}(B), SMatrix{1,2}(C)) + elseif size(C) == (2,1) + _sylvc!(C, SMatrix{2,2}(A), SMatrix{1,1}(B), SMatrix{2,1}(C)) + elseif size(C) == (2,2) + _sylvc!(C, SMatrix{2,2}(A), SMatrix{2,2}(B), SMatrix{2,2}(C)) + else + error("Matrix dimensionsins should not be greater than 2") + end + return C +end +function _sylvc!(a::Number, b::Number, c::SubArray{<:Any,0}) + c[1] = c[1] / (a + b) + + if isinf(c[1]); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end +end + + + +""" + + _sylvd!(X, A, B, C) + + Solve the discrete-time sylvester equation + + `AXB - X = C` + + for small matrices (1x1, 1x2, 2x1, 2x2) +""" +function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + if size(C) == (1,1) + _sylvd!(C, SMatrix{1,1}(A), SMatrix{1,1}(B), SMatrix{1,1}(C)) + elseif size(C) == (1,2) + _sylvd!(C, SMatrix{1,1}(A), SMatrix{2,2}(B), SMatrix{1,2}(C)) + elseif size(C) == (2,1) + _sylvd!(C, SMatrix{2,2}(A), SMatrix{1,1}(B), SMatrix{2,1}(C)) + elseif size(C) == (2,2) + _sylvd!(C, SMatrix{2,2}(A), SMatrix{2,2}(B), SMatrix{2,2}(C)) + else + error("Matrix dimensionsins should not be greater than 2") + end + return C +end +function _sylvd!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T <: Number} # Should really mostly be used for Reals + Cv = C[:] # vectorization of C + Xv = (kron(transpose(B), A) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) + + if any(isinf, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end + + X .= reshape(Xv, size(C,1), size(C,2)) +end +function _sylvd!(a::Number, b::Number, c::SubArray{<:Any,0}) + c[1] /= (a * b - 1) + + if isinf(c[1]); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end +end + + + +""" + Solve the continuous-time Sylvester equation + + `AX + XB = C` + + A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λ + μ = 0. + + [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +""" +function sylvc(A, B, C) + A2, UA = schur(Matrix(A')) + B2, UB = schur(B) + + C2 = UA'*C*UB # C2 should have the right type + + Y = _sylvc_schur!(Matrix(A2'), B2, C2, Val(:sylv), isreal(A2) ? Val(:real) : Val(:complex)) + + X = UA*Y*UB' +end + +""" + Solve the discrete-time Sylvester equation + + `AXB - X = C` + + A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λμ = 1. + + [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +""" +function sylvd(A, B, C) + A2, UA = schur(Matrix(A')) + B2, UB = schur(B) + + C2 = UA'*C*UB + + Y = _sylvd_schur!(Matrix(A2'), B2, C2, Val(:sylv), isreal(A2) ? Val(:real) : Val(:complex)) + + X = UA*Y*UB' +end + +""" + Solve the continuous-time Lyapunov equation + + `AX + XA' + Q = 0` + + A solution exists if ... FIXME + + [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +""" +function lyapc(A, Q) + + if !ishermitian(Q); error("The Q matrix must be Hermitian"); end + + At2, U = schur(Matrix(A')) + + Q2 = U'*Q*U # Could use in-place, Small savings could be possible here, see [1] + + Y = _sylvc_schur!(Matrix(At2'), At2, -Q2, Val(:lyap), isreal(At2) ? Val(:real) : Val(:complex)) + + X = U*Y*U' # Small savings could be possible here, see [1] +end +""" + `X = lyapd(A, Q)` + + Find the solution `X` to the discrete-time Lyapunov equation + + `AXA' - X + Q = 0` + + A solution `X` exists if `A` has no eigenvalue λ = ±1 and no eigenvalue pair λ₁λ₂ = 1 . + + + [1] **Barraud, A.** (1977) "A numerical algorithm to solve A'XA - X = Q" + IEEE Transactions on Automatic Control + + [2] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. + +""" +function lyapd(A, Q) + + if !ishermitian(Q); error("The Q matrix must be Hermitian"); end + + At2, U = schur(Matrix(A')) + + Q2 = U'*Q*U # Some savings could be possible here, see [1] + + Y = _sylvd_schur!(Matrix(At2'), At2, -Q2, Val(:lyap), isreal(At2) ? Val(:real) : Val(:complex)) + + X = U*Y*U' # Some savings could be possible here, see [1] +end + + +""" + Solve the continuous-time Sylvetser equation + + `AX + XB = C` + + where + `A` is assumed to have upper Schur form (lower quasi-triangular, 1x1 & 2x2 blocks on the diagonal) + `B` is assumed to have lower Schur form + + If `alg == Val(:lyap)` then `C` should be Hermitian + + See also `sylvc` +""" +# It should also hold that `eltype(C) = eltype(C / (A + B))` +function _sylvc_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} + # get block indices and nbr of blocks + realschur = (schurtype == Val(:real)) + if realschur + _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar + _, bb, nblocksb = _schurstructure(B, Val(:u)) + else + n = LinearAlgebra.checksquare(A) + ba, nblocksa = (FullStepRange(), n) + bb, nblocksb = (FullStepRange(), n) + end + + for j=1:nblocksb + i0 = (alg == Val(:lyap) ? j : 1) + for i=i0:nblocksa + + Cij = view(C, ba[i], bb[j]) + + if i > 1 + @views muladdrc!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1) # views ? + end + if j > 1 + @views muladdrc!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1) + end + + @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) + # Cij now contains the solution Yij + + if alg == Val(:lyap) && i > j + # adjoint is not defined for 0-dimensional views... + view(C, ba[j], bb[i]) .= (realschur ? Cij' : conj(Cij[1])) + end + + end + end + return C +end + + +# Should have eltype(C) = eltype(C / (A + B)) +function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} + n = LinearAlgebra.checksquare(A) + + # The matrix C is gradually replaced with the solution X + G = zeros(eltype(C), n, n) # This matrix contains A*X + + # get block dimensions, block indices, nbr of blocks + realschur = (schurtype == Val(:real)) + if realschur + _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar + _, bb, nblocksb = _schurstructure(B, Val(:u)) + else + ba, nblocksa = (FullStepRange(), size(A, 1)) + bb, nblocksb = (FullStepRange(), size(B, 1)) + end + + for j=1:nblocksb + i0 = (alg == Val(:lyap) ? j : 1) + for i=i0:nblocksa + Aii = A[ba[i], ba[i]] + Bjj = B[bb[j], bb[j]] + + Cij = view(C, ba[i], bb[j]) + Gij = view(G, ba[i], bb[j]) + + if i > 1 # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration + @views mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) + end + @views mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) + + + _sylvd!(Aii, Bjj, Cij) + # Cij now contains the solution Yij + + if alg == Val(:lyap) && i > j + # adjoint is not defined for 0-dimensional views... + view(C, ba[j], bb[i]) .= (realschur ? Cij' : conj(Cij[1])) + end + + mul!(Gij, Aii, Cij, 1, 1) + end + end + return C +end + +# Alternative version that avoids views, about twice as fast.. +function _sylvd_schurc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, ::Val{:complex}) where {T <: Number} + n = LinearAlgebra.checksquare(A) + + # The matrix C is gradually replaced with the solution X + G = zeros(eltype(C), n, n) # This matrix contains A*X + + # get block dimensions, block indices, nbr of blocks + for j=1:n + i0 = (alg == Val(:lyap) ? j : 1) + for i=i0:n + if i > 1 # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration + #@views G[i,j] += transpose(A[i, 1:i-1])*C[1:i-1, j] + G[i,j] += sum(A[i, k] * C[k, j] for k=1:i-1) + end + + C[i,j] -= sum(G[i, k] * B[k, j] for k=1:j) + + C[i,j] /= (A[i, i] * B[j, j] - 1) + # Cij now contains the solution Yij + + if isinf(C[i,j]); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end + + if alg == Val(:lyap) && i > j + # adjoint is not defined for 0-dimensional views... + C[j,i] = conj(C[i,j]) + end + + G[i,j] += A[i, i] * C[i, j] + end + end + return C +end + +end diff --git a/test/test_matrix_equations.jl b/test/test_matrix_equations.jl new file mode 100644 index 000000000..fb40f4889 --- /dev/null +++ b/test/test_matrix_equations.jl @@ -0,0 +1,119 @@ +# Test _schurstructure +R = diagm(0 => [1, 1, 1, 1, 1, 1, 1, 1], -1 => [1, 0, 0, 0, 1, 0, 0]) +d0 = [2, 1, 1, 2, 1, 1] +b0 = [1:2, 3:3, 4:4, 5:6, 7:7, 8:8] +@test (d0, b0, 6) == LyapTest._schurstructure(R, Val(:u)) +@test (d0, b0, 6) == LyapTest._schurstructure(R', Val(:l)) + + + +# Solve Sylvester equations with A = diagm(a) and B = diagm(b) +sylvc_diag = (a, b, C) -> [C[i,j]/(a[i] + b[j]) for i=1:length(a), j=1:length(b)] +sylvd_diag = (a, b, C) -> [C[i,j]/(a[i]*b[j] - 1) for i=1:length(a), j=1:length(b)] + + +# Basic tests using 1x1 Matrices + +a = [2.0]; A = diagm(a) +ac = [2.0 + im]; Ac = diagm(ac) +b = [3.0]; B = diagm(b) +C = fill(5.0, 1, 1) +Cc = fill(complex(5.0), 1, 1) + +@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvc_diag(a, b, C) +@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvc_diag(a, b, C) +@test LyapTest._sylvc_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvc_diag(ac, b, Cc) +@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:real)) ≈ sylvc_diag(a, a, C) +@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:complex)) ≈ sylvc_diag(a, a, C) +@test LyapTest._sylvc_schur!(Ac', Ac, copy(Cc), Val(:lyap), Val(:complex)) ≈ sylvc_diag(conj(ac), ac, Cc) + +@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvd_diag(a, b, C) +@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvd_diag(a, b, C) +@test LyapTest._sylvd_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvd_diag(ac, b, Cc) +@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:real)) ≈ sylvd_diag(a, a, C) +@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:complex)) ≈ sylvd_diag(a, a, C) +@test LyapTest._sylvd_schur!(Ac, copy(Ac'), copy(complex(C)), Val(:lyap), Val(:complex)) ≈ sylvd_diag(ac, conj(ac), C) + +@test LyapTest.sylvc(A, B, C) ≈ sylvc_diag(a, b, C) +@test LyapTest.sylvc(Ac, B, Cc) ≈ sylvc_diag(ac, b, Cc) +@test LyapTest.lyapc(A, C) ≈ sylvc_diag(a, a, -C) +@test LyapTest.lyapc(Ac, C) ≈ sylvc_diag(ac, conj(ac), -C) + +@test LyapTest.sylvd(A, B, C) ≈ sylvd_diag(a, b, C) +@test LyapTest.sylvd(Ac, B, Cc) ≈ sylvd_diag(ac, b, Cc) +@test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, a, -C) +@test LyapTest.lyapd(Ac, C) ≈ sylvd_diag(conj(ac), ac, -C) + + +# Basic tests with diagonal 2x2 matrices + +a = [2.0, 3]; A = diagm(a) +ac = a .+ im; Ac = diagm(ac) +b = [3.0, 5]; B = diagm(b) +C = [1.0 2; 2 1] +Cc = [1.0 2+im; 1+2im 1] +Cherm = [1.0 2+im; 2-im 1] + + +@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvc_diag(a, b, C) +@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvc_diag(a, b, C) +@test LyapTest._sylvc_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvc_diag(ac, b, Cc) +@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:real)) ≈ sylvc_diag(a, a, C) +@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:complex)) ≈ sylvc_diag(a, a, C) +@test LyapTest._sylvc_schur!(copy(Ac'), Ac, copy(Cherm), Val(:lyap), Val(:complex)) ≈ sylvc_diag(ac, conj(ac), Cherm) + +@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvd_diag(a, b, C) +@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvd_diag(a, b, C) +@test LyapTest._sylvd_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvd_diag(ac, b, Cc) +@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:real)) ≈ sylvd_diag(a, a, C) +@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:complex)) ≈ sylvd_diag(a, a, C) +@test LyapTest._sylvd_schur!(Ac, copy(Ac'), copy(Cherm), Val(:lyap), Val(:complex)) ≈ sylvd_diag(ac, conj(ac), Cherm) + +@test LyapTest.sylvc(A, B, C) ≈ sylvc_diag(a, b, C) +@test LyapTest.sylvc(Ac, B, Cc) ≈ sylvc_diag(ac, b, Cc) +@test LyapTest.lyapc(A, C) ≈ sylvc_diag(a, a, -C) +@test LyapTest.lyapc(Ac, Cherm) ≈ sylvc_diag(conj(ac), ac, -Cherm) + +@test LyapTest.sylvd(A, B, C) ≈ sylvd_diag(a, b, C) +@test LyapTest.sylvd(Ac, B, Cc) ≈ sylvd_diag(ac, b, Cc) +@test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, conj(a), -C) +@test LyapTest.lyapd(Ac, Cherm) ≈ sylvd_diag(ac, conj(ac), -Cherm) + + + +# +# Further tests with non-diagonal matrices +# + +# Helper function for evaluating the C matrix +sylvc_rhs = (A, B, X) -> (A*X + X*B) +sylvd_rhs = (A, B, X) -> (A*X*B - X) +lyapc_rhs = (A, X) -> -(A*X + X*A') +lyapd_rhs = (A, X) -> -(A*X*A' - X) + + +for (A, B, X) in [(fill(2.0, 1, 1), fill(3.0, 1, 1), fill(1.0, 1, 1)), + ([diagm([2.0,3]), diagm([4.0,5]), [1.0 2; 3 4]]), + ([-2 0; 2 3], [3.0 4; 0 5], [1.0 2; 3 4])] + + println(size(A)) + + Ac = A .* (1 + im) # Keep the structure of A.. + Xsym = X + X' + Xherm = (X .+ im) + (X .+ im)' + + @test LyapTest._sylvc_schur!(A, B, sylvc_rhs(A, B, X), Val(:sylv), Val(:real)) ≈ X + @test LyapTest._sylvc_schur!(A, B, sylvc_rhs(A, B, X), Val(:sylv), Val(:complex)) ≈ X + @test LyapTest._sylvc_schur!(Ac, B, sylvc_rhs(Ac, B, X), Val(:sylv), Val(:complex)) ≈ X + @test LyapTest._sylvc_schur!(A, copy(A'), -lyapc_rhs(A, Xsym), Val(:lyap), Val(:real)) ≈ Xsym + @test LyapTest._sylvc_schur!(A, copy(A'), -lyapc_rhs(A, Xsym), Val(:lyap), Val(:complex)) ≈ Xsym + @test LyapTest._sylvc_schur!(Ac, copy(Ac'), -lyapc_rhs(Ac, Xherm), Val(:lyap), Val(:complex)) ≈ Xherm + + @test LyapTest._sylvd_schur!(A, B, sylvd_rhs(A, B, X), Val(:sylv), Val(:real)) ≈ X + @test LyapTest._sylvd_schur!(A, B, sylvd_rhs(A, B, X), Val(:sylv), Val(:complex)) ≈ X + @test LyapTest._sylvd_schur!(Ac, B, sylvd_rhs(Ac, B, X), Val(:sylv), Val(:complex)) ≈ X + @test LyapTest._sylvd_schur!(A, copy(A'), -lyapd_rhs(A, Xsym), Val(:lyap), Val(:real)) ≈ Xsym + @test LyapTest._sylvd_schur!(A, copy(A'), -lyapd_rhs(A, Xsym), Val(:lyap), Val(:complex)) ≈ Xsym + @test LyapTest._sylvd_schur!(Ac, copy(Ac'), -lyapd_rhs(Ac, Xherm), Val(:lyap), Val(:complex)) ≈ Xherm + +end From e566545ed6afbbda242a699e19b08ea2ce58a8b1 Mon Sep 17 00:00:00 2001 From: olof3 Date: Tue, 5 May 2020 16:34:17 +0200 Subject: [PATCH 02/13] Fix. --- src/LyapTest.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 690bc1f45..4ce9d4eb3 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -328,9 +328,9 @@ function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Uni Gij = view(G, ba[i], bb[j]) if i > 1 # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration - @views mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) + @views muladdrc!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) end - @views mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) + @views muladdrc!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) _sylvd!(Aii, Bjj, Cij) From b6258615704b263ddae61e15ab3d17507543e1ca Mon Sep 17 00:00:00 2001 From: olof3 Date: Thu, 7 May 2020 13:36:30 +0200 Subject: [PATCH 03/13] Switched to different implementations for real and complex case. --- src/LyapTest.jl | 168 ++++++++++++++++++------------------------------ 1 file changed, 64 insertions(+), 104 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 4ce9d4eb3..49a37c2c5 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -8,6 +8,7 @@ using StaticArrays * Perhaps it would be better with separate methods for the real and complex versions of the block solve algorithms * What convention for signs and tranposes? +* Should error checking be done? In that case where? =# """ @@ -49,28 +50,6 @@ function _schurstructure(R::AbstractMatrix, ul=Val(:u)::Union{Val{:u}, Val{:l}}) return d, b, k end -# A step range representing -Inf:1:Inf -# for efficient reuse of blocked real algorithms for scalar complex version -struct FullStepRange end -Base.getindex(::FullStepRange, k::Int) = k - -# If the C entry is a 1x1 (0-dimenional view), then -# vector vector multiplication is interpreted as tranpose(A)*B -muladdrc!(C, A, B, α, β) = mul!(C, A, B, α, β) -function muladdrc!(c::Base.SubArray{<:Any,0}, A::Base.SubArray{<:Any,1}, B::Base.SubArray{<:Any,1}, α, β) - c[1] *= β - if eltype(A) <: ComplexF64 && eltype(B) <: ComplexF64 - c[1] += α*LinearAlgebra.BLAS.dotu(A, B) - else - c[1] += α*sum(A[k]*B[k] for k=1:length(A)) - end -end -function muladdrc!(c::Base.SubArray{<:Any,0}, A::Base.SubArray{<:Any,0}, B, α, β) - c[1] *= β - c[1] += α*A[1]*B[1] -end - - """ @@ -87,7 +66,7 @@ function _sylvc!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T m, n = size(C) Xv = (kron(I(n), A) + kron(transpose(B), I(m))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) - if any(isinf, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end + if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end X .= SMatrix{size(C,1),size(C,2)}(reshape(Xv, size(C,1), size(C,2))) end @@ -105,12 +84,6 @@ function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) end return C end -function _sylvc!(a::Number, b::Number, c::SubArray{<:Any,0}) - c[1] = c[1] / (a + b) - - if isinf(c[1]); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end -end - """ @@ -141,15 +114,11 @@ function _sylvd!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T Cv = C[:] # vectorization of C Xv = (kron(transpose(B), A) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) - if any(isinf, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end + if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end X .= reshape(Xv, size(C,1), size(C,2)) end -function _sylvd!(a::Number, b::Number, c::SubArray{<:Any,0}) - c[1] /= (a * b - 1) - if isinf(c[1]); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end -end @@ -173,6 +142,13 @@ function sylvc(A, B, C) X = UA*Y*UB' end +@inline function sylvc(a::Number, b::Number, c::Number) + x = c / (a + b) + + if !isfinite(x); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end + + x +end """ Solve the discrete-time Sylvester equation @@ -194,6 +170,13 @@ function sylvd(A, B, C) X = UA*Y*UB' end +@inline function sylvd(a::Number, b::Number, c::Number) + x = c / (a * b - 1) + + if !isfinite(x); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end + + x +end """ Solve the continuous-time Lyapunov equation @@ -264,118 +247,95 @@ end # It should also hold that `eltype(C) = eltype(C / (A + B))` function _sylvc_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} # get block indices and nbr of blocks - realschur = (schurtype == Val(:real)) - if realschur + if schurtype == Val(:real) _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar _, bb, nblocksb = _schurstructure(B, Val(:u)) else - n = LinearAlgebra.checksquare(A) - ba, nblocksa = (FullStepRange(), n) - bb, nblocksb = (FullStepRange(), n) + nblocksa = size(A, 1) + nblocksb = size(B, 1) end for j=1:nblocksb i0 = (alg == Val(:lyap) ? j : 1) for i=i0:nblocksa + if schurtype == Val(:complex) + if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end + if j > 1; C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end - Cij = view(C, ba[i], bb[j]) + C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] - if i > 1 - @views muladdrc!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1) # views ? - end - if j > 1 - @views muladdrc!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1) - end + if alg == Val(:lyap) && i > j + C[j,i] = conj(C[i,j]) + end + else + Cij = view(C, ba[i], bb[j]) - @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) - # Cij now contains the solution Yij + if i > 1; @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end + if j > 1; @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end - if alg == Val(:lyap) && i > j - # adjoint is not defined for 0-dimensional views... - view(C, ba[j], bb[i]) .= (realschur ? Cij' : conj(Cij[1])) - end + @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij + if alg == Val(:lyap) && i > j + view(C, ba[j], bb[i]) .= Cij' + end + end end end return C end - # Should have eltype(C) = eltype(C / (A + B)) function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} - n = LinearAlgebra.checksquare(A) - # The matrix C is gradually replaced with the solution X - G = zeros(eltype(C), n, n) # This matrix contains A*X + + G = zeros(eltype(C), size(A,1), size(B, 1)) # This matrix contains A*X for increased performance # get block dimensions, block indices, nbr of blocks - realschur = (schurtype == Val(:real)) - if realschur + if schurtype == Val(:real) _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar _, bb, nblocksb = _schurstructure(B, Val(:u)) else - ba, nblocksa = (FullStepRange(), size(A, 1)) - bb, nblocksb = (FullStepRange(), size(B, 1)) + nblocksa = size(A, 1) + nblocksb = size(B, 1) end for j=1:nblocksb i0 = (alg == Val(:lyap) ? j : 1) for i=i0:nblocksa - Aii = A[ba[i], ba[i]] - Bjj = B[bb[j], bb[j]] - - Cij = view(C, ba[i], bb[j]) - Gij = view(G, ba[i], bb[j]) - - if i > 1 # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration - @views muladdrc!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) - end - @views muladdrc!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) + if schurtype == Val(:complex) + # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration + if i > 1; G[i,j] += sum(A[i,k] * C[k,j] for k=1:i-1); end + C[i,j] -= sum(G[i,k] * B[k,j] for k=1:j) - _sylvd!(Aii, Bjj, Cij) - # Cij now contains the solution Yij + C[i,j] = sylvd(A[i,i], B[j,j], C[i,j]) # C[i,j] now contains solution Y[i,j] - if alg == Val(:lyap) && i > j - # adjoint is not defined for 0-dimensional views... - view(C, ba[j], bb[i]) .= (realschur ? Cij' : conj(Cij[1])) - end + if alg == Val(:lyap) && i > j + C[j,i] = conj(C[i,j]) + end - mul!(Gij, Aii, Cij, 1, 1) - end - end - return C -end + G[i,j] += A[i, i] * C[i, j] + else + Aii = A[ba[i], ba[i]] + Bjj = B[bb[j], bb[j]] -# Alternative version that avoids views, about twice as fast.. -function _sylvd_schurc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, ::Val{:complex}) where {T <: Number} - n = LinearAlgebra.checksquare(A) + Cij = view(C, ba[i], bb[j]) + Gij = view(G, ba[i], bb[j]) - # The matrix C is gradually replaced with the solution X - G = zeros(eltype(C), n, n) # This matrix contains A*X + if i > 1 + @views mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) + end - # get block dimensions, block indices, nbr of blocks - for j=1:n - i0 = (alg == Val(:lyap) ? j : 1) - for i=i0:n - if i > 1 # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration - #@views G[i,j] += transpose(A[i, 1:i-1])*C[1:i-1, j] - G[i,j] += sum(A[i, k] * C[k, j] for k=1:i-1) - end + @views mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) - C[i,j] -= sum(G[i, k] * B[k, j] for k=1:j) + _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij - C[i,j] /= (A[i, i] * B[j, j] - 1) - # Cij now contains the solution Yij + if alg == Val(:lyap) && i > j + view(C, ba[j], bb[i]) .= Cij' + end - if isinf(C[i,j]); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end - - if alg == Val(:lyap) && i > j - # adjoint is not defined for 0-dimensional views... - C[j,i] = conj(C[i,j]) + mul!(Gij, Aii, Cij, 1, 1) end - - G[i,j] += A[i, i] * C[i, j] end end return C From 22fbbcbb38c95ac4817a35f495279bbeea1d4aa1 Mon Sep 17 00:00:00 2001 From: olof3 Date: Mon, 11 May 2020 16:31:48 +0200 Subject: [PATCH 04/13] Updates. --- lyap_bench.jl | 154 ++++++++++++++++++ src/LyapTest.jl | 289 ++++++++++++++++++++++------------ test/test_matrix_equations.jl | 1 - 3 files changed, 345 insertions(+), 99 deletions(-) create mode 100644 lyap_bench.jl diff --git a/lyap_bench.jl b/lyap_bench.jl new file mode 100644 index 000000000..4a3f9f3ed --- /dev/null +++ b/lyap_bench.jl @@ -0,0 +1,154 @@ +using MatrixEquations +include("src/LyapTest.jl") + +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2 + +bg = BenchmarkGroup() +for T in [Float64, ComplexF64] + for k=[5,50,500]#[5, 20, 200] + n = k + A = randn(T, n, n) + B = randn(T, n, n) + C = randn(T, n, n) + Csym = Matrix(Hermitian(randn(T, n, n))) + + bg["1. sylvc", 1, T, n] = @benchmarkable LyapTest.sylvc($A, $B, $C) + bg["1. sylvc", 2, T, n] = @benchmarkable MatrixEquations.sylvc($A, $B, $C) + bg["1. sylvc", 3, T, n] = @benchmarkable sylvester($A, $B, $C) + + bg["3. sylvd", 1, T, n] = @benchmarkable LyapTest.sylvd($A, $B, $C) + bg["3. sylvd", 2, T, n] = @benchmarkable MatrixEquations.sylvd($A, $B, $C) + + bg["2. lyapc", 1, T, n] = @benchmarkable LyapTest.lyapc($A, $Csym) + bg["2. lyapc", 2, T, n] = @benchmarkable MatrixEquations.lyapc($A, $Csym) + bg["2. lyapc", 3, T, n] = @benchmarkable lyap($A, $Csym) + + bg["4. lyapd", 1, T, n] = @benchmarkable LyapTest.lyapd($A, $Csym) + bg["4. lyapd", 2, T, n] = @benchmarkable MatrixEquations.lyapd($A, $Csym) + end +end + +results = run(bg) + + +using PrettyTables + + +sorted_results = sort(collect(results.data), by=x->(x[1][1], x[1][3]==Float64, x[1][2], x[1][4])) + +sr = [sorted_results; fill(NaN, 12)] + +for k=1:4 + data = reshape(sr[9*(k-1) .+ (1:9)], 3, 3) + + println("$(data[1][1][1]), $(data[1][1][3])") + table_content = Matrix{Any}(fill(NaN, 3, 4)) + table_content[1:3,1] = [x[1][4] for x in sorted_results[1:3]] + table_content[1:3,2:4] = [string(BenchmarkTools.prettytime(median(x[2]).time), " (", x[2].allocs, " allocs, ", BenchmarkTools.prettymemory(x[2].memory), ")") for x in data] + + pretty_table(table_content, ["", "LyapTest", "MatrixEquations", "julia"]) +end + +f = sr[1][2] +f + +table_content = Matrix{Any}(fill(NaN, 8, 4)) +table_content[1:8,1] = ["$(x[1][1]), $(x[1][3])" for x in sorted_results[1:8]] +table_content[1:8,2:3] = [string(median(x[2])) for x in sorted_results[1:16]] +table_content[1:4,4] = [string(median(x[2])) for x in sorted_results[17:20]] + +pretty_table(table_content) + + +@pt table_content + + + + +n = 100 +for (numtype, T) in [(:real, Float64), (:complex, ComplexF64)] + A = randn(T, n, n) + B = randn(T, n, n) + C = randn(T, n, n) + Q = Matrix(Hermitian(randn(T, n, n))) + + println("") + println("sylvc ($numtype)") + X = LyapTest.sylvc(A, B, C) + @printf("LyapTest err=%.2e\n", norm(A*X + X*B - C)) + X = MatrixEquations.sylvc(A, B, C) + @printf("MatrixEqs err=%.2e\n", norm(A*X + X*B - C)) + X = LinearAlgebra.sylvester(A, B, C) + @printf("LinearAlgebra err=%.2e\n", norm(A*X + X*B + C)) + + println("\nlyapc ($numtype)") + X = LyapTest.lyapc(A, Q) + @printf("LyapTest err=%.2e\n", norm(A*X + X*A' + Q)) + X = MatrixEquations.lyapc(A, Q) + @printf("MatrixEqs err=%.2e\n", norm(A*X + X*A' + Q)) + X = LinearAlgebra.lyap(A, Q) + @printf("LinearAlgebra err=%.2e\n", norm(A*X + X*A' + Q)) + + println("\nsylvd ($numtype)") + X = LyapTest.sylvd(A, B, C) + @printf("LyapTest err=%.2e\n", norm(A*X*B - X - C)) + X = MatrixEquations.sylvd(A, B, C) + @printf("MatrixEqs err=%.2e\n", norm(A*X*B + X - C)) + + println("\nlyapd ($numtype)") + X = LyapTest.lyapd(A, Q) + @printf("LyapTest err=%.2e\n", norm(A*X*A' - X + Q)) + X = MatrixEquations.lyapd(A, Q) + @printf("MatrixEqs err=%.2e\n", norm(A*X*A' - X + Q)) +end + + +using TimerOutputs +const to = TimerOutput() + +for (numtype, T) in [(:real, Float64), (:complex, ComplexF64)] + n = 100 + A = Matrix(schur(randn(T, n, n)).T') # Lower quasi triangular + B = schur(randn(T, n, n)).T + C = randn(T, n, n) + Csym = Matrix(Hermitian(randn(T, n, n))) + + @timeit to string(numtype) begin + println("_sylvc_schur (:sylv) ($numtype)") + + X_sylvc = LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(numtype)) + @printf(" err=%.2e\n", norm(A*X_sylvc + X_sylvc*B - C)) + + X_sylvc = copy(C) + MatrixEquations.sylvcs!(Matrix(A'), B, X_sylvc) + @printf(" err=%.2e\n", norm(A'*X_sylvc + X_sylvc*B - C)) + + println("_sylvd_schur (:sylv) ($numtype)") + + X_sylvd = LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(numtype)) + @printf(" err=%.2e\n", norm(A*X_sylvd*B - X_sylvd - C)) + + X_sylvd = copy(C) + MatrixEquations.sylvds!(Matrix(A'), B, X_sylvd) + @printf(" err=%.2e\n", norm(A'*X_sylvd*B + X_sylvd - C)) + + + println("_sylvc_schur (:lyap) ($numtype)") + + X_lyapc = LyapTest._sylvc_schur!(A, Matrix(A'), copy(Csym), Val(:lyap), Val(numtype)) + @printf(" err=%.2e\n", norm(A*X_lyapc + X_lyapc*Matrix(A') - Csym)) + + X_lyapc = copy(Csym) + MatrixEquations.lyapcs!(B, X_lyapc) + @printf(" err=%.2e\n", norm(B*X_lyapc + X_lyapc*Matrix(B') + Csym)) + + println("_sylvd_schur (:lyap) ($numtype)") + + X_lyapd = LyapTest._sylvd_schur!(Matrix(B'), B, copy(Csym), Val(:lyap), Val(numtype)) + @printf(" err=%.2e\n", norm(B'*X_lyapd*B - X_lyapd - Csym)) + + X_lyapd = copy(Csym) + MatrixEquations.lyapds!(B, X_lyapd) + @printf(" err=%.2e\n", norm(B*X_lyapd*B' - X_lyapd + Csym)) + end +end diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 49a37c2c5..79a4fc6be 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -1,8 +1,5 @@ module LyapTest -using LinearAlgebra -using StaticArrays - #= * Symmetry is not exploited when solving diagonal entries * Perhaps it would be better with separate methods for the real and @@ -11,19 +8,23 @@ using StaticArrays * Should error checking be done? In that case where? =# + +using LinearAlgebra +using StaticArrays + """ `b, d, nblocks = _schurstructure(R::AbstractMatrix)` - Return the block strucutre of a quasi-traingular Schur matrix `R`. +Return the block strucutre of a quasi-traingular Schur matrix `R`. - `d` contains the block sizees of each diagonal block (`1` or `2`) +`d` contains the block sizees of each diagonal block (`1` or `2`) - `b` contains the indices of the blocks +`b` contains the indices of the blocks - `nblocks` is the number of blocks +`nblocks` is the number of blocks """ -function _schurstructure(R::AbstractMatrix, ul=Val(:u)::Union{Val{:u}, Val{:l}}) +function _schurstructure(R::AbstractMatrix, ul=Val(:U)::Union{Val{:U}, Val{:L}}) n = size(R,1) d = Vector{Int}(undef, n) # block sizes @@ -36,7 +37,7 @@ function _schurstructure(R::AbstractMatrix, ul=Val(:u)::Union{Val{:u}, Val{:l}}) if j == n d[k] = 1 else - if ul == Val(:u) + if ul == Val(:U) d[k] = iszero(R[j+1, j]) ? 1 : 2 else d[k] = iszero(R[j, j+1]) ? 1 : 2 @@ -50,95 +51,123 @@ function _schurstructure(R::AbstractMatrix, ul=Val(:u)::Union{Val{:u}, Val{:l}}) return d, b, k end +# FIXME: better handling of uniform scaling?! +issquare(A::Number) = true +issquare(A::AbstractMatrix) = size(A,1) == size(A,2) +function _check_lyap_inputs(A, Q) + if !issquare(A); error("The A matrix must be square"); end + if !ishermitian(Q); error("The Q matrix must be Hermitian"); end + if size(Q, 1) != size(A, 1); error("The A and Q matrices must have the same dimensions"); end +end -""" +function _check_sylv_inputs(A, B, C) + if !issquare(A); error("The A matrix must be square"); end + if !issquare(B); error("The B matrix must be square"); end + if size(C, 1) != size(A, 1); error("The A and C matrices have inconsistent dimensions"); end + if size(C, 2) != size(B, 2); error("The B and C matrices have inconsistent dimensions"); end +end - _sylvc!(A, B, C) +# Should preferably be fixed in LinearAlgebra +LinearAlgebra.schur(A::AbstractMatrix{T}) where T = schur!(LinearAlgebra.copy_oftype(A, LinearAlgebra.eigtype(T))) - Solve the continuous-time sylvester equation - `AX + XB = C` +sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C)) +sylvdsoltype(A, B, C) = Base.promote_op(sylvd, eltype(A), eltype(B), eltype(C)) - for small matrices (1x1, 1x2, 2x1, 2x2) """ -function _sylvc!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T <: Number} # Should really mostly be used for Reals - Cv = C[:] # vectorization of C - m, n = size(C) - Xv = (kron(I(n), A) + kron(transpose(B), I(m))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) + _sylvc!(A, B, C) - if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end +Solve the continuous-time sylvester equation - X .= SMatrix{size(C,1),size(C,2)}(reshape(Xv, size(C,1), size(C,2))) -end +`AX + XB = C` + +for small matrices (1x1, 1x2, 2x1, 2x2) +""" function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) - if size(C) == (1,1) - _sylvc!(C, SMatrix{1,1}(A), SMatrix{1,1}(B), SMatrix{1,1}(C)) - elseif size(C) == (1,2) - _sylvc!(C, SMatrix{1,1}(A), SMatrix{2,2}(B), SMatrix{1,2}(C)) - elseif size(C) == (2,1) - _sylvc!(C, SMatrix{2,2}(A), SMatrix{1,1}(B), SMatrix{2,1}(C)) - elseif size(C) == (2,2) - _sylvc!(C, SMatrix{2,2}(A), SMatrix{2,2}(B), SMatrix{2,2}(C)) + M, N = size(C) + if M == 2 && N == 2 + _sylvc!(C, A, B, C, Val{2}(), Val{2}()) + elseif M == 2 && N == 1 + _sylvc!(C, A, B, C, Val{2}(), Val{1}()) + elseif M == 1 && N == 2 + _sylvc!(C, A, B, C, Val{1}(), Val{2}()) + elseif M == 1 && N == 1 + _sylvc!(C, A, B, C, Val{1}(), Val{1}()) else error("Matrix dimensionsins should not be greater than 2") end return C end +function _sylvc!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} + A′ = SMatrix{M,M}(A) + B′ = SMatrix{N,N}(B) + Cv = SVector{M*N}(C[:]) # vectorization of C + Xv = lu(kron(SMatrix{N,N}(I), A′) + kron(transpose(B′), SMatrix{M,M}(I))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) (with A = I or B = I) -""" + if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end + + X .= reshape(Xv, M, N) +end + +""" _sylvd!(X, A, B, C) - Solve the discrete-time sylvester equation +Solve the discrete-time sylvester equation - `AXB - X = C` +`AXB - X = C` - for small matrices (1x1, 1x2, 2x1, 2x2) +for small matrices (1x1, 1x2, 2x1, 2x2). """ function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) - if size(C) == (1,1) - _sylvd!(C, SMatrix{1,1}(A), SMatrix{1,1}(B), SMatrix{1,1}(C)) - elseif size(C) == (1,2) - _sylvd!(C, SMatrix{1,1}(A), SMatrix{2,2}(B), SMatrix{1,2}(C)) - elseif size(C) == (2,1) - _sylvd!(C, SMatrix{2,2}(A), SMatrix{1,1}(B), SMatrix{2,1}(C)) - elseif size(C) == (2,2) - _sylvd!(C, SMatrix{2,2}(A), SMatrix{2,2}(B), SMatrix{2,2}(C)) + M, N = size(C) + if M == 2 && N == 2 + _sylvd!(C, A, B, C, Val{2}(), Val{2}()) + elseif M == 2 && N == 1 + _sylvd!(C, A, B, C, Val{2}(), Val{1}()) + elseif M == 1 && N == 2 + _sylvd!(C, A, B, C, Val{1}(), Val{2}()) + elseif M == 1 && N == 1 + _sylvd!(C, A, B, C, Val{1}(), Val{1}()) else error("Matrix dimensionsins should not be greater than 2") end return C end -function _sylvd!(X::AbstractMatrix, A::SMatrix, B::SMatrix, C::SMatrix) where {T <: Number} # Should really mostly be used for Reals - Cv = C[:] # vectorization of C - Xv = (kron(transpose(B), A) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) +function _sylvd!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} + A′ = SMatrix{M,M}(A) + B′ = SMatrix{N,N}(B) + Cv = SVector{M*N}(C[:]) # vectorization of C + + Xv = (kron(transpose(B′), A′) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end - X .= reshape(Xv, size(C,1), size(C,2)) + X .= reshape(Xv, M, N) end - - - """ - Solve the continuous-time Sylvester equation + sylvc(A, B, C) - `AX + XB = C` +Solve the continuous-time Sylvester equation - A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λ + μ = 0. +`AX + XB = C` - [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix - equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λ + μ = 0. + +[1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. """ function sylvc(A, B, C) - A2, UA = schur(Matrix(A')) + _check_sylv_inputs(A, B, C) + + At2, UA = schur(A') B2, UB = schur(B) C2 = UA'*C*UB # C2 should have the right type - Y = _sylvc_schur!(Matrix(A2'), B2, C2, Val(:sylv), isreal(A2) ? Val(:real) : Val(:complex)) + Y = _sylvc_schur!(Matrix(At2'), B2, C2, Val(:sylv)) X = UA*Y*UB' end @@ -151,22 +180,26 @@ end end """ - Solve the discrete-time Sylvester equation + sylvd(A, B, C) - `AXB - X = C` +Solve the discrete-time Sylvester equation - A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λμ = 1. +`AXB - X = C` - [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix - equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λμ = 1. + +[1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. """ function sylvd(A, B, C) - A2, UA = schur(Matrix(A')) + _check_sylv_inputs(A, B, C) + + At2, UA = schur(A') B2, UB = schur(B) C2 = UA'*C*UB - Y = _sylvd_schur!(Matrix(A2'), B2, C2, Val(:sylv), isreal(A2) ? Val(:real) : Val(:complex)) + Y = _sylvd_schur!(Matrix(At2'), B2, C2, Val(:sylv)) X = UA*Y*UB' end @@ -179,77 +212,100 @@ end end """ - Solve the continuous-time Lyapunov equation + lyapc(A, Q) - `AX + XA' + Q = 0` +Computes the solution `X` of the continuous-time Lyapunov equation - A solution exists if ... FIXME +`AX + XA' + Q = 0` - [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix - equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +A solution exists unless `A` has an eigenvalue λ = ±1 or an eigenvalue pair λ₁λ₂ = 1. + +[1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. """ function lyapc(A, Q) - if !ishermitian(Q); error("The Q matrix must be Hermitian"); end + _check_lyap_inputs(A, Q) - At2, U = schur(Matrix(A')) + At2, U = schur(A') Q2 = U'*Q*U # Could use in-place, Small savings could be possible here, see [1] - Y = _sylvc_schur!(Matrix(At2'), At2, -Q2, Val(:lyap), isreal(At2) ? Val(:real) : Val(:complex)) + Y = _sylvc_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap)) X = U*Y*U' # Small savings could be possible here, see [1] end """ `X = lyapd(A, Q)` - Find the solution `X` to the discrete-time Lyapunov equation +Find the solution `X` to the discrete-time Lyapunov equation - `AXA' - X + Q = 0` +`AXA' - X + Q = 0` - A solution `X` exists if `A` has no eigenvalue λ = ±1 and no eigenvalue pair λ₁λ₂ = 1 . +A solution `X` unless `A` has an eigenvalue λ = ±1 or an eigenvalue pair λ₁λ₂ = 1 . - [1] **Barraud, A.** (1977) "A numerical algorithm to solve A'XA - X = Q" - IEEE Transactions on Automatic Control +[1] **Barraud, A.** (1977) "A numerical algorithm to solve A'XA - X = Q" + IEEE Transactions on Automatic Control - [2] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix - equation AX + XB = C" Communications of the ACM, 15(9), 820-826. +[2] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix + equation AX + XB = C" Communications of the ACM, 15(9), 820-826. """ function lyapd(A, Q) - if !ishermitian(Q); error("The Q matrix must be Hermitian"); end + _check_lyap_inputs(A, Q) - At2, U = schur(Matrix(A')) + At2, U = schur(A') + #tmp = similar(Q, sylvdsoltype(A,A,Q)) + #Q2 = U'*mul!(tmp, Q, U) # Some savings could be possible here, see [1] Q2 = U'*Q*U # Some savings could be possible here, see [1] - Y = _sylvd_schur!(Matrix(At2'), At2, -Q2, Val(:lyap), isreal(At2) ? Val(:real) : Val(:complex)) + Y = _sylvd_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap)) - X = U*Y*U' # Some savings could be possible here, see [1] + X = U*Y*U' # mul!(Y.data, U, mul!(tmp, Y, U')) # Y*U*U' # Some savings could be possible here, see [1] end """ - Solve the continuous-time Sylvetser equation +Solve the continuous-time Sylvester equation - `AX + XB = C` +`AX + XB = C` - where - `A` is assumed to have upper Schur form (lower quasi-triangular, 1x1 & 2x2 blocks on the diagonal) - `B` is assumed to have lower Schur form +where +`A` is assumed to have lower Schur form (quasi-triangular, 1x1 & 2x2 blocks on the diagonal) +`B` is assumed to have upper Schur form - If `alg == Val(:lyap)` then `C` should be Hermitian +If `alg == Val(:lyap)` then `C` should be Hermitian - See also `sylvc` +The input matrix `C` is overwritten with the solution `X`, +so this matrix should have the same type as the solution, typically `typeof(C[1] / (A[1] + B[1]))` + +See also `sylvc` """ -# It should also hold that `eltype(C) = eltype(C / (A + B))` -function _sylvc_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} +function sylvc_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + _check_sylv_inputs(A, B, C) + T = sylvcsoltype(A, B, C) + _sylvc_schur!(convert(Matrix, A), convert(Matrix, B), convert(Matrix{T}, C), Val(:sylv)) +end +function lyapc_schur!(A::AbstractMatrix, Q::AbstractMatrix) + _check_lyap_inputs(A, Q) + T = sylvcsoltype(A, A, Q) + _sylvc_schur!(convert(Matrix, A), convert(Matrix, A'), lmul!(-1, Matrix{T}(Q)), Val(:lyap)) +end +function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Val{:lyap}}, + schurtype::Union{Val{:real},Val{:complex}} = isreal(A) || isreal(B) ? Val(:real) : Val(:complex)) where {T <: Number} + # The user should preferably use sylvc_schur! and lyapc_schur! + # I.e., this method does not check whether C is hermitian + # The matrix C is successively replaced with the solution X + # if alg == Val(:lyap), only the lower triangle of C is computed + # after which an Hermitian view is applied + # get block indices and nbr of blocks if schurtype == Val(:real) - _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar - _, bb, nblocksb = _schurstructure(B, Val(:u)) + _, ba, nblocksa = _schurstructure(A, Val(:L)) # A is assumed upper triangualar + _, bb, nblocksb = _schurstructure(B, Val(:U)) else nblocksa = size(A, 1) nblocksb = size(B, 1) @@ -276,6 +332,7 @@ function _sylvc_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Uni @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j + #copyto_noalias!(view(C, ba[j], bb[i]), Cij') view(C, ba[j], bb[i]) .= Cij' end end @@ -284,16 +341,40 @@ function _sylvc_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Uni return C end -# Should have eltype(C) = eltype(C / (A + B)) -function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Union{Val{:sylv},Val{:lyap}}, schurtype::Union{Val{:real},Val{:complex}}) where {T <: Number} - # The matrix C is gradually replaced with the solution X +""" +Solve the discrete-time Sylvester equation + +`AXB - X = C` - G = zeros(eltype(C), size(A,1), size(B, 1)) # This matrix contains A*X for increased performance +where `A` is assumed to have lower Schur form (quasi-triangular, 1x1 & 2x2 blocks on the diagonal) +`B` is assumed to have upper Schur form + +If `alg == Val(:lyap)` then `C` should be Hermitian + +The input matrix `C` is overwritten with the solution `X`, +so this matrix should have the same type as the solution, typically `typeof(C[1] / (A[1] * B[1] - 1)` + +See also `sylvd` +""" +function sylvd_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + _check_sylv_inputs(A, B, C) + T = sylvdsoltype(A, B, C) + _sylvd_schur!(convert(Matrix, A), convert(Matrix, B), convert(Matrix{T}, C), Val(:sylv)) +end +function lyapd_schur!(A::AbstractMatrix, Q::AbstractMatrix) + _check_lyap_inputs(A, Q) + T = sylvdsoltype(A, A, Q) + _sylvd_schur!(convert(Matrix, A), convert(Matrix, A'), lmul!(-1, Matrix{T}(Q)), Val(:lyap)) +end +function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Val{:lyap}}, + schurtype::Union{Val{:real},Val{:complex}} = isreal(A) || isreal(B) ? Val(:real) : Val(:complex)) where {T <: Number} + + G = zeros(eltype(C), size(A,1), size(B, 1)) # Keep track of A*X for improved performance # get block dimensions, block indices, nbr of blocks if schurtype == Val(:real) - _, ba, nblocksa = _schurstructure(A, Val(:l)) # A is assumed upper triangualar - _, bb, nblocksb = _schurstructure(B, Val(:u)) + _, ba, nblocksa = _schurstructure(A, Val(:L)) # A is assumed upper triangualar + _, bb, nblocksb = _schurstructure(B, Val(:U)) else nblocksa = size(A, 1) nblocksb = size(B, 1) @@ -331,6 +412,7 @@ function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Uni _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j + #copyto_noalias!(view(C, ba[j], bb[i]), Cij') view(C, ba[j], bb[i]) .= Cij' end @@ -341,4 +423,15 @@ function _sylvd_schur!(A::AbstractMatrix, B::Matrix, C::AbstractMatrix, alg::Uni return C end + + + +function copyto_noalias!(dest::AbstractArray{T1,N}, src::AbstractArray{T2,N}) where {T1,T2,N} + for I in eachindex(IndexStyle(src,dest), src) + @inbounds dest[I] = src[I] + end + dest +end + + end diff --git a/test/test_matrix_equations.jl b/test/test_matrix_equations.jl index fb40f4889..cbbaf8b38 100644 --- a/test/test_matrix_equations.jl +++ b/test/test_matrix_equations.jl @@ -80,7 +80,6 @@ Cherm = [1.0 2+im; 2-im 1] @test LyapTest.lyapd(Ac, Cherm) ≈ sylvd_diag(ac, conj(ac), -Cherm) - # # Further tests with non-diagonal matrices # From 45360fcbe279443a60bf5d7e2c47f4a9f22995dc Mon Sep 17 00:00:00 2001 From: olof3 Date: Mon, 11 May 2020 22:05:00 +0200 Subject: [PATCH 05/13] Updates. --- lyap_bench.jl | 18 ++- src/LyapTest.jl | 39 +++---- test/test_matrix_equations.jl | 207 +++++++++++++++++++++------------- 3 files changed, 165 insertions(+), 99 deletions(-) diff --git a/lyap_bench.jl b/lyap_bench.jl index 4a3f9f3ed..e6131bcae 100644 --- a/lyap_bench.jl +++ b/lyap_bench.jl @@ -1,7 +1,7 @@ using MatrixEquations include("src/LyapTest.jl") -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2 +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5 bg = BenchmarkGroup() for T in [Float64, ComplexF64] @@ -34,9 +34,9 @@ results = run(bg) using PrettyTables -sorted_results = sort(collect(results.data), by=x->(x[1][1], x[1][3]==Float64, x[1][2], x[1][4])) +sr = sort(collect(results.data), by=x->(x[1][1], x[1][3]==Float64, x[1][2], x[1][4])) + -sr = [sorted_results; fill(NaN, 12)] for k=1:4 data = reshape(sr[9*(k-1) .+ (1:9)], 3, 3) @@ -49,6 +49,18 @@ for k=1:4 pretty_table(table_content, ["", "LyapTest", "MatrixEquations", "julia"]) end +for k=1:4 + data = reshape(sr[36 + 6*(k-1) .+ (1:6)], 3, 2) + + #println(data) + println("$(data[1][1][1]), $(data[1][1][3])") + table_content = Matrix{Any}(fill(NaN, 3, 3)) + table_content[1:3,1] = [x[1][4] for x in sorted_results[1:3]] + table_content[1:3,2:3] = [string(BenchmarkTools.prettytime(median(x[2]).time), " (", x[2].allocs, " allocs, ", BenchmarkTools.prettymemory(x[2].memory), ")") for x in data] + + pretty_table(table_content, ["", "LyapTest", "MatrixEquations", "julia"]) +end + f = sr[1][2] f diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 79a4fc6be..1902342d4 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -1,6 +1,7 @@ module LyapTest #= + * Symmetry is not exploited when solving diagonal entries * Perhaps it would be better with separate methods for the real and complex versions of the block solve algorithms @@ -8,10 +9,11 @@ module LyapTest * Should error checking be done? In that case where? =# - using LinearAlgebra using StaticArrays + + """ `b, d, nblocks = _schurstructure(R::AbstractMatrix)` @@ -83,7 +85,7 @@ Solve the continuous-time sylvester equation for small matrices (1x1, 1x2, 2x1, 2x2) """ -function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) +@inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) if M == 2 && N == 2 _sylvc!(C, A, B, C, Val{2}(), Val{2}()) @@ -98,10 +100,10 @@ function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) end return C end -function _sylvc!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} +@inline function _sylvc!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} A′ = SMatrix{M,M}(A) B′ = SMatrix{N,N}(B) - Cv = SVector{M*N}(C[:]) # vectorization of C + Cv = SMatrix{M,N}(C)[:] # vectorization of C Xv = lu(kron(SMatrix{N,N}(I), A′) + kron(transpose(B′), SMatrix{M,M}(I))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) (with A = I or B = I) @@ -138,7 +140,7 @@ end function _sylvd!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} A′ = SMatrix{M,M}(A) B′ = SMatrix{N,N}(B) - Cv = SVector{M*N}(C[:]) # vectorization of C + Cv = SMatrix{M,N}(C)[:] # vectorization of C Xv = (kron(transpose(B′), A′) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) @@ -165,11 +167,11 @@ function sylvc(A, B, C) At2, UA = schur(A') B2, UB = schur(B) - C2 = UA'*C*UB # C2 should have the right type + C2 = UA'*C*UB # This should give C2 the right type Y = _sylvc_schur!(Matrix(At2'), B2, C2, Val(:sylv)) - X = UA*Y*UB' + X = mul!(Y, UA, Y*UB') end @inline function sylvc(a::Number, b::Number, c::Number) x = c / (a + b) @@ -201,7 +203,7 @@ function sylvd(A, B, C) Y = _sylvd_schur!(Matrix(At2'), B2, C2, Val(:sylv)) - X = UA*Y*UB' + X = mul!(Y, UA, Y*UB') end @inline function sylvd(a::Number, b::Number, c::Number) x = c / (a * b - 1) @@ -233,7 +235,7 @@ function lyapc(A, Q) Y = _sylvc_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap)) - X = U*Y*U' # Small savings could be possible here, see [1] + X = mul!(Y, U, Y*U') # Small savings could be possible here, see [1] end """ `X = lyapd(A, Q)` @@ -258,13 +260,11 @@ function lyapd(A, Q) At2, U = schur(A') - #tmp = similar(Q, sylvdsoltype(A,A,Q)) - #Q2 = U'*mul!(tmp, Q, U) # Some savings could be possible here, see [1] - Q2 = U'*Q*U # Some savings could be possible here, see [1] + Q2 = U'*Q*U Y = _sylvd_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap)) - X = U*Y*U' # mul!(Y.data, U, mul!(tmp, Y, U')) # Y*U*U' # Some savings could be possible here, see [1] + X = mul!(Y, U, Y*U') end @@ -326,14 +326,15 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va else Cij = view(C, ba[i], bb[j]) - if i > 1; @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end - if j > 1; @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end + + if i > 1; @inbounds @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end + if j > 1; @inbounds @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij + if alg == Val(:lyap) && i > j - #copyto_noalias!(view(C, ba[j], bb[i]), Cij') - view(C, ba[j], bb[i]) .= Cij' + copyto_noalias!(view(C, ba[j], bb[i]), Cij') end end end @@ -341,6 +342,7 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va return C end + """ Solve the discrete-time Sylvester equation @@ -412,8 +414,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - #copyto_noalias!(view(C, ba[j], bb[i]), Cij') - view(C, ba[j], bb[i]) .= Cij' + copyto_noalias!(view(C, ba[j], bb[i]), Cij') end mul!(Gij, Aii, Cij, 1, 1) diff --git a/test/test_matrix_equations.jl b/test/test_matrix_equations.jl index cbbaf8b38..679a3a6df 100644 --- a/test/test_matrix_equations.jl +++ b/test/test_matrix_equations.jl @@ -1,52 +1,54 @@ +using StaticArrays +using Test +using Printf +using Random + +using MatrixEquations + +Random.seed!(0) + + +include("src/LyapTest.jl") + # Test _schurstructure R = diagm(0 => [1, 1, 1, 1, 1, 1, 1, 1], -1 => [1, 0, 0, 0, 1, 0, 0]) d0 = [2, 1, 1, 2, 1, 1] b0 = [1:2, 3:3, 4:4, 5:6, 7:7, 8:8] -@test (d0, b0, 6) == LyapTest._schurstructure(R, Val(:u)) -@test (d0, b0, 6) == LyapTest._schurstructure(R', Val(:l)) - +@test (d0, b0, 6) == LyapTest._schurstructure(R, Val(:U)) +@test (d0, b0, 6) == LyapTest._schurstructure(R', Val(:L)) -# Solve Sylvester equations with A = diagm(a) and B = diagm(b) -sylvc_diag = (a, b, C) -> [C[i,j]/(a[i] + b[j]) for i=1:length(a), j=1:length(b)] -sylvd_diag = (a, b, C) -> [C[i,j]/(a[i]*b[j] - 1) for i=1:length(a), j=1:length(b)] +# Test inner small dimension sylv solver +for T = [Float64, ComplexF64] + for m=1:2, n=1:2 + A = randn(T, m, m) + B = randn(T, n, n) + C = randn(T, m, n) + X = LyapTest._sylvd!(A, B, copy(C)) -# Basic tests using 1x1 Matrices + @test Matrix(A*X*B - X - C) ≈ zeros(m,n) atol=1e-13 -a = [2.0]; A = diagm(a) -ac = [2.0 + im]; Ac = diagm(ac) -b = [3.0]; B = diagm(b) -C = fill(5.0, 1, 1) -Cc = fill(complex(5.0), 1, 1) + X = LyapTest._sylvc!(A, B, copy(C)) -@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvc_diag(a, b, C) -@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvc_diag(a, b, C) -@test LyapTest._sylvc_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvc_diag(ac, b, Cc) -@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:real)) ≈ sylvc_diag(a, a, C) -@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:complex)) ≈ sylvc_diag(a, a, C) -@test LyapTest._sylvc_schur!(Ac', Ac, copy(Cc), Val(:lyap), Val(:complex)) ≈ sylvc_diag(conj(ac), ac, Cc) + @test Matrix(A*X + X*B - C) ≈ zeros(m,n) atol=1e-13 + end +end -@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvd_diag(a, b, C) -@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvd_diag(a, b, C) -@test LyapTest._sylvd_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvd_diag(ac, b, Cc) -@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:real)) ≈ sylvd_diag(a, a, C) -@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:complex)) ≈ sylvd_diag(a, a, C) -@test LyapTest._sylvd_schur!(Ac, copy(Ac'), copy(complex(C)), Val(:lyap), Val(:complex)) ≈ sylvd_diag(ac, conj(ac), C) -@test LyapTest.sylvc(A, B, C) ≈ sylvc_diag(a, b, C) -@test LyapTest.sylvc(Ac, B, Cc) ≈ sylvc_diag(ac, b, Cc) -@test LyapTest.lyapc(A, C) ≈ sylvc_diag(a, a, -C) -@test LyapTest.lyapc(Ac, C) ≈ sylvc_diag(ac, conj(ac), -C) - -@test LyapTest.sylvd(A, B, C) ≈ sylvd_diag(a, b, C) -@test LyapTest.sylvd(Ac, B, Cc) ≈ sylvd_diag(ac, b, Cc) -@test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, a, -C) -@test LyapTest.lyapd(Ac, C) ≈ sylvd_diag(conj(ac), ac, -C) +# Solve Sylvester equations with A = diagm(a) and B = diagm(b) +sylvc_diag = (a, b, C) -> [C[i,j]/(a[i] + b[j]) for i=1:length(a), j=1:length(b)] +sylvd_diag = (a, b, C) -> [C[i,j]/(a[i]*b[j] - 1) for i=1:length(a), j=1:length(b)] +# Compute C from X to allow convenient check of solution +sylvc_rhs = (A, B, X) -> (A*X + X*B) +sylvd_rhs = (A, B, X) -> (A*X*B - X) +lyapc_rhs = (A, X) -> -(A*X + X*A') +lyapd_rhs = (A, X) -> -(A*X*A' - X) +# # Basic tests with diagonal 2x2 matrices - +# a = [2.0, 3]; A = diagm(a) ac = a .+ im; Ac = diagm(ac) b = [3.0, 5]; B = diagm(b) @@ -54,65 +56,116 @@ C = [1.0 2; 2 1] Cc = [1.0 2+im; 1+2im 1] Cherm = [1.0 2+im; 2-im 1] - -@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvc_diag(a, b, C) -@test LyapTest._sylvc_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvc_diag(a, b, C) -@test LyapTest._sylvc_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvc_diag(ac, b, Cc) -@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:real)) ≈ sylvc_diag(a, a, C) -@test LyapTest._sylvc_schur!(A, A, copy(C), Val(:lyap), Val(:complex)) ≈ sylvc_diag(a, a, C) -@test LyapTest._sylvc_schur!(copy(Ac'), Ac, copy(Cherm), Val(:lyap), Val(:complex)) ≈ sylvc_diag(ac, conj(ac), Cherm) - -@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:real)) ≈ sylvd_diag(a, b, C) -@test LyapTest._sylvd_schur!(A, B, copy(C), Val(:sylv), Val(:complex)) ≈ sylvd_diag(a, b, C) -@test LyapTest._sylvd_schur!(Ac, B, copy(Cc), Val(:sylv), Val(:complex)) ≈ sylvd_diag(ac, b, Cc) -@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:real)) ≈ sylvd_diag(a, a, C) -@test LyapTest._sylvd_schur!(A, copy(A'), copy(C), Val(:lyap), Val(:complex)) ≈ sylvd_diag(a, a, C) -@test LyapTest._sylvd_schur!(Ac, copy(Ac'), copy(Cherm), Val(:lyap), Val(:complex)) ≈ sylvd_diag(ac, conj(ac), Cherm) - +@test LyapTest.sylvc_schur!(A, B, copy(C)) ≈ sylvc_diag(a, b, C) +@test LyapTest.sylvc_schur!(Ac, B, copy(Cc)) ≈ sylvc_diag(ac, b, Cc) @test LyapTest.sylvc(A, B, C) ≈ sylvc_diag(a, b, C) @test LyapTest.sylvc(Ac, B, Cc) ≈ sylvc_diag(ac, b, Cc) -@test LyapTest.lyapc(A, C) ≈ sylvc_diag(a, a, -C) -@test LyapTest.lyapc(Ac, Cherm) ≈ sylvc_diag(conj(ac), ac, -Cherm) +@test LyapTest.sylvd_schur!(A, B, copy(C)) ≈ sylvd_diag(a, b, C) +@test LyapTest.sylvd_schur!(Ac, B, copy(Cc)) ≈ sylvd_diag(ac, b, Cc) @test LyapTest.sylvd(A, B, C) ≈ sylvd_diag(a, b, C) @test LyapTest.sylvd(Ac, B, Cc) ≈ sylvd_diag(ac, b, Cc) -@test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, conj(a), -C) + +@test LyapTest.lyapc_schur!(A, copy(C)) ≈ sylvc_diag(a, a, -C) +@test LyapTest.lyapc_schur!(Ac, copy(Cherm)) ≈ sylvc_diag(conj(ac), ac, -Cherm) +@test LyapTest.lyapc(A, C) ≈ sylvc_diag(a, a, -C) +@test LyapTest.lyapc(Ac, Cherm) ≈ sylvc_diag(conj(ac), ac, -Cherm) + +@test LyapTest.lyapd_schur!(A, copy(C)) ≈ sylvd_diag(a, a, -C) +@test LyapTest.lyapd_schur!(Ac, copy(Cherm)) ≈ sylvd_diag(ac, conj(ac), -Cherm) +@test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, a, -C) @test LyapTest.lyapd(Ac, Cherm) ≈ sylvd_diag(ac, conj(ac), -Cherm) + + # -# Further tests with non-diagonal matrices +# Tests with small trinagular matrices of various types # +A3 = [2 0 0; 3 4 0; 5 6 7] +B3 = [2 3 4; 0 1 1; 0 0 2] +X3 = reshape(1:9, 3, 3) + +for (A, B, X, tsname) in [(fill(2.0, 1, 1), fill(3.0, 1, 1), fill(1.0, 1, 1), "1x1"), + (diagm([2.0,3]), diagm([4.0,5]), [1.0 2; 3 4], "2x2 diagonal"), + ([-2 0; 2 3], [3.0 4; 0 5], [1.0 2; 3 4], "2x2 tringular"), + ([-2.0 0; 2 3], fill(5.0, 1, 1), [1.0; 2][:,:], "2x1"), + ([-2 0; 2 3], fill(5, 1, 1), [1; 2][:,:], "2x1"), + (big.([-2.0 0; 2 3]), fill(5//2, 1, 1), [1; 2][:,:], "2x1"), + (fill(2.0, 1, 1), [3.0 4; 0 5], [1.0 2], "1x2"), + (big.(A3), Rational.(B3), X3, "3x3"), + ] + + # Generate complex version, while keeping the structure + Ac = A .* (1 + im) + Bc = B .* (1 + im) + Xc = X .* (1 + im) + + @testset "sylv(c/d)_schur!, $(tsname) $(eltype.((A, B, X)))" begin + @test LyapTest.sylvc_schur!(A, B, sylvc_rhs(A, B, X)) ≈ X + @test LyapTest.sylvc_schur!(Ac, B, sylvc_rhs(Ac, B, X)) ≈ X + @test LyapTest.sylvc_schur!(A, B, sylvc_rhs(A, B, Xc)) ≈ Xc + @test LyapTest.sylvc_schur!(Ac, Bc, sylvc_rhs(Ac, Bc, Xc)) ≈ Xc + + @test LyapTest.sylvd_schur!(A, B, sylvd_rhs(A, B, X)) ≈ X + @test LyapTest.sylvd_schur!(Ac, B, sylvd_rhs(Ac, B, X)) ≈ X + @test LyapTest.sylvd_schur!(A, B, sylvd_rhs(A, B, Xc)) ≈ Xc + @test LyapTest.sylvd_schur!(Ac, Bc, sylvd_rhs(Ac, Bc, Xc)) ≈ Xc + + end + + if size(X,1) == size(X,2) + Xherm = X + X' + Xcherm = Xc + Xc' + + @testset "lyap(c/d)_schur!, $(tsname) $(eltype.((A, C)))" begin + @test LyapTest.lyapc_schur!(A, lyapc_rhs(A, Xherm)) ≈ Xherm + @test LyapTest.lyapc_schur!(Ac, lyapc_rhs(Ac, Xherm)) ≈ Xherm + @test LyapTest.lyapc_schur!(A, lyapc_rhs(A, Xcherm)) ≈ Xcherm + @test LyapTest.lyapc_schur!(Ac, lyapc_rhs(Ac, Xcherm)) ≈ Xcherm + + @test LyapTest.lyapd_schur!(A, lyapd_rhs(A, Xherm)) ≈ Xherm + @test LyapTest.lyapd_schur!(Ac, lyapd_rhs(Ac, Xherm)) ≈ Xherm + @test LyapTest.lyapd_schur!(A, lyapd_rhs(A, Xcherm)) ≈ Xcherm + @test LyapTest.lyapd_schur!(Ac, lyapd_rhs(Ac, Xcherm)) ≈ Xcherm + end + end +end -# Helper function for evaluating the C matrix -sylvc_rhs = (A, B, X) -> (A*X + X*B) -sylvd_rhs = (A, B, X) -> (A*X*B - X) -lyapc_rhs = (A, X) -> -(A*X + X*A') -lyapd_rhs = (A, X) -> -(A*X*A' - X) +Random.seed!(0) +for (m,n) in [(1,1), (1,2), (2,1), (2,2), (3,3), (3,2), (2,3), (5,5), (7,8), (100,100)] + for (TA, TB, TC) in [(Float64, Float64, Float64), + (ComplexF64, Float64, Float64), + (Float64, ComplexF64, Float64), + (Float64, ComplexF64, ComplexF64), + (ComplexF64, ComplexF64, ComplexF64), + (Int64, ComplexF64, ComplexF64), + (Int64, Int64, Int64)] + + # Generate complex version, while keeping the structure + A = TA <: Int ? rand(1:50, m, m) : randn(TA, m, m) + B = TB <: Int ? rand(1:50, n, n) : randn(TB, n, n) + C = TC <: Int ? rand(1:50, m, n) : randn(TC, m, n) -for (A, B, X) in [(fill(2.0, 1, 1), fill(3.0, 1, 1), fill(1.0, 1, 1)), - ([diagm([2.0,3]), diagm([4.0,5]), [1.0 2; 3 4]]), - ([-2 0; 2 3], [3.0 4; 0 5], [1.0 2; 3 4])] + rtol = m > 50 ? 1e-10 : 1e-12 - println(size(A)) + @testset "lyap/sylv ($m, $n, $TA, $TB, $TC)" begin + X = LyapTest.sylvc(A, B, C) + @test A*X + X*B ≈ C rtol=rtol - Ac = A .* (1 + im) # Keep the structure of A.. - Xsym = X + X' - Xherm = (X .+ im) + (X .+ im)' + X = LyapTest.sylvd(A, B, C) + @test A*X*B - X ≈ C rtol=rtol - @test LyapTest._sylvc_schur!(A, B, sylvc_rhs(A, B, X), Val(:sylv), Val(:real)) ≈ X - @test LyapTest._sylvc_schur!(A, B, sylvc_rhs(A, B, X), Val(:sylv), Val(:complex)) ≈ X - @test LyapTest._sylvc_schur!(Ac, B, sylvc_rhs(Ac, B, X), Val(:sylv), Val(:complex)) ≈ X - @test LyapTest._sylvc_schur!(A, copy(A'), -lyapc_rhs(A, Xsym), Val(:lyap), Val(:real)) ≈ Xsym - @test LyapTest._sylvc_schur!(A, copy(A'), -lyapc_rhs(A, Xsym), Val(:lyap), Val(:complex)) ≈ Xsym - @test LyapTest._sylvc_schur!(Ac, copy(Ac'), -lyapc_rhs(Ac, Xherm), Val(:lyap), Val(:complex)) ≈ Xherm + if m == n + Q = Matrix(Hermitian(randn(m,m))) - @test LyapTest._sylvd_schur!(A, B, sylvd_rhs(A, B, X), Val(:sylv), Val(:real)) ≈ X - @test LyapTest._sylvd_schur!(A, B, sylvd_rhs(A, B, X), Val(:sylv), Val(:complex)) ≈ X - @test LyapTest._sylvd_schur!(Ac, B, sylvd_rhs(Ac, B, X), Val(:sylv), Val(:complex)) ≈ X - @test LyapTest._sylvd_schur!(A, copy(A'), -lyapd_rhs(A, Xsym), Val(:lyap), Val(:real)) ≈ Xsym - @test LyapTest._sylvd_schur!(A, copy(A'), -lyapd_rhs(A, Xsym), Val(:lyap), Val(:complex)) ≈ Xsym - @test LyapTest._sylvd_schur!(Ac, copy(Ac'), -lyapd_rhs(Ac, Xherm), Val(:lyap), Val(:complex)) ≈ Xherm + X = LyapTest.lyapc(A, Q) + @test A*X + X*A' ≈ -Q rtol=rtol + X = LyapTest.lyapd(A, Q) + @test A*X*A' - X ≈ -Q rtol=rtol + end + end + end end From 5965620626a26759f5d598800f1eb759013d3f5f Mon Sep 17 00:00:00 2001 From: olof3 Date: Tue, 12 May 2020 11:15:08 +0200 Subject: [PATCH 06/13] Some minor fixes. --- lyap_bench.jl | 17 +---------------- src/LyapTest.jl | 42 ++++++++++++++++++++---------------------- 2 files changed, 21 insertions(+), 38 deletions(-) diff --git a/lyap_bench.jl b/lyap_bench.jl index e6131bcae..ca073a8f7 100644 --- a/lyap_bench.jl +++ b/lyap_bench.jl @@ -58,24 +58,9 @@ for k=1:4 table_content[1:3,1] = [x[1][4] for x in sorted_results[1:3]] table_content[1:3,2:3] = [string(BenchmarkTools.prettytime(median(x[2]).time), " (", x[2].allocs, " allocs, ", BenchmarkTools.prettymemory(x[2].memory), ")") for x in data] - pretty_table(table_content, ["", "LyapTest", "MatrixEquations", "julia"]) + pretty_table(table_content, ["", "LyapTest", "MatrixEquations"]) end -f = sr[1][2] -f - -table_content = Matrix{Any}(fill(NaN, 8, 4)) -table_content[1:8,1] = ["$(x[1][1]), $(x[1][3])" for x in sorted_results[1:8]] -table_content[1:8,2:3] = [string(median(x[2])) for x in sorted_results[1:16]] -table_content[1:4,4] = [string(median(x[2])) for x in sorted_results[17:20]] - -pretty_table(table_content) - - -@pt table_content - - - n = 100 for (numtype, T) in [(:real, Float64), (:complex, ComplexF64)] diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 1902342d4..cf98997c5 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -79,7 +79,7 @@ sylvdsoltype(A, B, C) = Base.promote_op(sylvd, eltype(A), eltype(B), eltype(C)) """ _sylvc!(A, B, C) -Solve the continuous-time sylvester equation +Find the solution `X` to the continuous-time Sylvester equation `AX + XB = C` @@ -88,19 +88,19 @@ for small matrices (1x1, 1x2, 2x1, 2x2) @inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) if M == 2 && N == 2 - _sylvc!(C, A, B, C, Val{2}(), Val{2}()) + _sylvc!(A, B, C, Val{2}(), Val{2}()) elseif M == 2 && N == 1 - _sylvc!(C, A, B, C, Val{2}(), Val{1}()) + _sylvc!(A, B, C, Val{2}(), Val{1}()) elseif M == 1 && N == 2 - _sylvc!(C, A, B, C, Val{1}(), Val{2}()) + _sylvc!(A, B, C, Val{1}(), Val{2}()) elseif M == 1 && N == 1 - _sylvc!(C, A, B, C, Val{1}(), Val{1}()) + _sylvc!(A, B, C, Val{1}(), Val{1}()) else error("Matrix dimensionsins should not be greater than 2") end return C end -@inline function _sylvc!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} +@inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} A′ = SMatrix{M,M}(A) B′ = SMatrix{N,N}(B) Cv = SMatrix{M,N}(C)[:] # vectorization of C @@ -109,14 +109,14 @@ end if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end - X .= reshape(Xv, M, N) + C .= reshape(Xv, M, N) end """ _sylvd!(X, A, B, C) -Solve the discrete-time sylvester equation +Find the solution `X` to the discrete-time Sylvester equation `AXB - X = C` @@ -125,19 +125,19 @@ for small matrices (1x1, 1x2, 2x1, 2x2). function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) if M == 2 && N == 2 - _sylvd!(C, A, B, C, Val{2}(), Val{2}()) + _sylvd!(A, B, C, Val{2}(), Val{2}()) elseif M == 2 && N == 1 - _sylvd!(C, A, B, C, Val{2}(), Val{1}()) + _sylvd!(A, B, C, Val{2}(), Val{1}()) elseif M == 1 && N == 2 - _sylvd!(C, A, B, C, Val{1}(), Val{2}()) + _sylvd!(A, B, C, Val{1}(), Val{2}()) elseif M == 1 && N == 1 - _sylvd!(C, A, B, C, Val{1}(), Val{1}()) + _sylvd!(A, B, C, Val{1}(), Val{1}()) else error("Matrix dimensionsins should not be greater than 2") end return C end -function _sylvd!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} +function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} A′ = SMatrix{M,M}(A) B′ = SMatrix{N,N}(B) Cv = SMatrix{M,N}(C)[:] # vectorization of C @@ -146,17 +146,17 @@ function _sylvd!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::Abs if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end - X .= reshape(Xv, M, N) + C .= reshape(Xv, M, N) end """ sylvc(A, B, C) -Solve the continuous-time Sylvester equation +Find the solution `X` to the continuous-time Sylvester equation `AX + XB = C` -A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λ + μ = 0. +A solution exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λ + μ = 0. [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix equation AX + XB = C" Communications of the ACM, 15(9), 820-826. @@ -184,11 +184,11 @@ end """ sylvd(A, B, C) -Solve the discrete-time Sylvester equation +Find the solution `X` to the discrete-time Sylvester equation `AXB - X = C` -A solution `X` exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λμ = 1. +A solution exists unless `A` and `B` have eigenvalues `λ` and `μ` such that λμ = 1. [1] **Bartels, R. H., & Stewart, G. W.** (1972). "Solution of the matrix equation AX + XB = C" Communications of the ACM, 15(9), 820-826. @@ -244,7 +244,7 @@ Find the solution `X` to the discrete-time Lyapunov equation `AXA' - X + Q = 0` -A solution `X` unless `A` has an eigenvalue λ = ±1 or an eigenvalue pair λ₁λ₂ = 1 . +A solution exists unless `A` has an eigenvalue λ = ±1 or an eigenvalue pair λ₁λ₂ = 1 . [1] **Barraud, A.** (1977) "A numerical algorithm to solve A'XA - X = Q" @@ -326,13 +326,11 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va else Cij = view(C, ba[i], bb[j]) - if i > 1; @inbounds @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end if j > 1; @inbounds @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij - if alg == Val(:lyap) && i > j copyto_noalias!(view(C, ba[j], bb[i]), Cij') end @@ -414,7 +412,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - copyto_noalias!(view(C, ba[j], bb[i]), Cij') + copyto_noalias!(view(C, ba[j], bb[i]), Cij') end mul!(Gij, Aii, Cij, 1, 1) From d0d4f211ced40102f94def2380b28488626e039d Mon Sep 17 00:00:00 2001 From: olof3 Date: Wed, 13 May 2020 09:38:54 +0200 Subject: [PATCH 07/13] Minor updates. --- src/LyapTest.jl | 62 ++++++++++++++++++++++++------------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index cf98997c5..9fd9888c9 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -1,12 +1,17 @@ module LyapTest #= - -* Symmetry is not exploited when solving diagonal entries -* Perhaps it would be better with separate methods for the real and - complex versions of the block solve algorithms +The implementaiton is intended to be simple and fast. However in order to maintain + readability, a few small optimizations have not been persued, for example: +* Symmetry is not exploited when solving diagonal entries of Lyapuniv equations +* Could use one more in-place multiplication in the top level functions +* Small savings could be possible when compting U*Q*U' where Q is symmetric, see [1] +* Could make the traiangular solve methods accepts various combinations of upper/lower + Schur matrices but this seems to add quite a bit of code at quite a bit of code for quite a small gain. + +# Some other questions * What convention for signs and tranposes? -* Should error checking be done? In that case where? +* Should error checking be done after each block solve or at the end of the algorithm? =# using LinearAlgebra @@ -15,9 +20,11 @@ using StaticArrays """ - `b, d, nblocks = _schurstructure(R::AbstractMatrix)` + `_schurstructure(R::AbstractMatrix, ul=Union{Val{:U}, Val{:L}}) -> (b, d, nblocks)` + +Return the block strucutre of an upper quasi-traingular Schur matrix `R`. +`ul` indicates if R is upper (`ul=Val(:U)`) or lower (`ul=Val(:L)`) triangular. -Return the block strucutre of a quasi-traingular Schur matrix `R`. `d` contains the block sizees of each diagonal block (`1` or `2`) @@ -32,7 +39,7 @@ function _schurstructure(R::AbstractMatrix, ul=Val(:U)::Union{Val{:U}, Val{:L}}) d = Vector{Int}(undef, n) # block sizes b = Vector{UnitRange{Int64}}(undef, n) # block indices - j = 1 # column if ul=:u, row if ul=:l + j = 1 # column if ul=:U, row if ul=:L k = 0 # block number while j <= n k += 1 @@ -231,11 +238,11 @@ function lyapc(A, Q) At2, U = schur(A') - Q2 = U'*Q*U # Could use in-place, Small savings could be possible here, see [1] + Q2 = U'*Q*U Y = _sylvc_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap)) - X = mul!(Y, U, Y*U') # Small savings could be possible here, see [1] + X = mul!(Y, U, Y*U') end """ `X = lyapd(A, Q)` @@ -269,19 +276,13 @@ end """ -Solve the continuous-time Sylvester equation +Find the solution `X` to the continuous-time Sylvester equation `AX + XB = C` -where -`A` is assumed to have lower Schur form (quasi-triangular, 1x1 & 2x2 blocks on the diagonal) +where `A` is assumed to have lower Schur form (quasi-triangular, 1x1 & 2x2 blocks on the diagonal) `B` is assumed to have upper Schur form -If `alg == Val(:lyap)` then `C` should be Hermitian - -The input matrix `C` is overwritten with the solution `X`, -so this matrix should have the same type as the solution, typically `typeof(C[1] / (A[1] + B[1]))` - See also `sylvc` """ function sylvc_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) @@ -315,10 +316,10 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va i0 = (alg == Val(:lyap) ? j : 1) for i=i0:nblocksa if schurtype == Val(:complex) - if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end - if j > 1; C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end + if i > 1; @inbounds C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end + if j > 1; @inbounds C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end - C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] + @inbounds C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] if alg == Val(:lyap) && i > j C[j,i] = conj(C[i,j]) @@ -342,6 +343,8 @@ end """ + sylvd_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + Solve the discrete-time Sylvester equation `AXB - X = C` @@ -349,10 +352,7 @@ Solve the discrete-time Sylvester equation where `A` is assumed to have lower Schur form (quasi-triangular, 1x1 & 2x2 blocks on the diagonal) `B` is assumed to have upper Schur form -If `alg == Val(:lyap)` then `C` should be Hermitian - -The input matrix `C` is overwritten with the solution `X`, -so this matrix should have the same type as the solution, typically `typeof(C[1] / (A[1] * B[1] - 1)` +If the matrix `C` has the right type, it is overwritten with the solution `X`. See also `sylvd` """ @@ -397,19 +397,19 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va G[i,j] += A[i, i] * C[i, j] else - Aii = A[ba[i], ba[i]] - Bjj = B[bb[j], bb[j]] - Cij = view(C, ba[i], bb[j]) Gij = view(G, ba[i], bb[j]) + Aii = A[ba[i], ba[i]] + Bjj = B[bb[j], bb[j]] + if i > 1 - @views mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) + @views @inbounds mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) end - @views mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) + @views @inbounds mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) - _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij + @inbounds _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j copyto_noalias!(view(C, ba[j], bb[i]), Cij') From 7667539e97a2e3c847b5ffe6f293ee855cbd7eb5 Mon Sep 17 00:00:00 2001 From: olof3 Date: Wed, 13 May 2020 12:50:56 +0200 Subject: [PATCH 08/13] Removed copyto_noalias! --- src/LyapTest.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 9fd9888c9..1c014d48a 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -333,7 +333,9 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - copyto_noalias!(view(C, ba[j], bb[i]), Cij') + @inbounds for l=bb[j], k=ba[i] + C[l,k] = conj(C[k,l]) + end end end end @@ -412,7 +414,9 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va @inbounds _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - copyto_noalias!(view(C, ba[j], bb[i]), Cij') + @inbounds for l=bb[j], k=ba[i] # Avoids aliasing of copyto! + C[l,k] = conj(C[k,l]) + end end mul!(Gij, Aii, Cij, 1, 1) @@ -424,13 +428,4 @@ end - -function copyto_noalias!(dest::AbstractArray{T1,N}, src::AbstractArray{T2,N}) where {T1,T2,N} - for I in eachindex(IndexStyle(src,dest), src) - @inbounds dest[I] = src[I] - end - dest -end - - end From 67af15acc356317a93fb6e7758ea8601e50b79d7 Mon Sep 17 00:00:00 2001 From: olof3 Date: Wed, 13 May 2020 17:58:05 +0200 Subject: [PATCH 09/13] Updates. --- lyap_bench.jl | 9 ++++----- src/LyapTest.jl | 34 ++++++++++++++++++---------------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/lyap_bench.jl b/lyap_bench.jl index ca073a8f7..ef2480226 100644 --- a/lyap_bench.jl +++ b/lyap_bench.jl @@ -34,12 +34,10 @@ results = run(bg) using PrettyTables -sr = sort(collect(results.data), by=x->(x[1][1], x[1][3]==Float64, x[1][2], x[1][4])) - - +sorted_results = sort(collect(results.data), by=x->(x[1][1], x[1][3]==Float64, x[1][2], x[1][4])) for k=1:4 - data = reshape(sr[9*(k-1) .+ (1:9)], 3, 3) + data = reshape(sorted_results[9*(k-1) .+ (1:9)], 3, 3) println("$(data[1][1][1]), $(data[1][1][3])") table_content = Matrix{Any}(fill(NaN, 3, 4)) @@ -50,7 +48,7 @@ for k=1:4 end for k=1:4 - data = reshape(sr[36 + 6*(k-1) .+ (1:6)], 3, 2) + data = reshape(sorted_results[36 + 6*(k-1) .+ (1:6)], 3, 2) #println(data) println("$(data[1][1][1]), $(data[1][1][3])") @@ -63,6 +61,7 @@ end n = 100 + for (numtype, T) in [(:real, Float64), (:complex, ComplexF64)] A = randn(T, n, n) B = randn(T, n, n) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 1c014d48a..f8957f57e 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -312,28 +312,30 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va nblocksb = size(B, 1) end - for j=1:nblocksb + @inbounds for j=1:nblocksb i0 = (alg == Val(:lyap) ? j : 1) for i=i0:nblocksa if schurtype == Val(:complex) - if i > 1; @inbounds C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end - if j > 1; @inbounds C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end + if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end + if j > 1; C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end - @inbounds C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] + C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] if alg == Val(:lyap) && i > j C[j,i] = conj(C[i,j]) end else + Aii = view(A, ba[i], ba[i]) + Bjj = view(B, bb[j], bb[j]) Cij = view(C, ba[i], bb[j]) - if i > 1; @inbounds @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end - if j > 1; @inbounds @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end + if i > 1; @views mul!(Cij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], -1, 1); end + if j > 1; @views mul!(Cij, C[ba[i], 1:bb[j-1][end]], B[1:bb[j-1][end], bb[j]], -1, 1); end - @views _sylvc!(A[ba[i], ba[i]], B[bb[j], bb[j]], Cij) # Cij now contains the solution Yij + _sylvc!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - @inbounds for l=bb[j], k=ba[i] + for l=bb[j], k=ba[i] C[l,k] = conj(C[k,l]) end end @@ -382,7 +384,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va nblocksb = size(B, 1) end - for j=1:nblocksb + @inbounds for j=1:nblocksb i0 = (alg == Val(:lyap) ? j : 1) for i=i0:nblocksa if schurtype == Val(:complex) @@ -399,22 +401,22 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va G[i,j] += A[i, i] * C[i, j] else + Aii = view(A, ba[i], ba[i]) + Bjj = view(B, bb[j], bb[j]) Cij = view(C, ba[i], bb[j]) - Gij = view(G, ba[i], bb[j]) - Aii = A[ba[i], ba[i]] - Bjj = B[bb[j], bb[j]] + Gij = view(G, ba[i], bb[j]) if i > 1 - @views @inbounds mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) + @views mul!(Gij, A[ba[i], 1:ba[i-1][end]], C[1:ba[i-1][end], bb[j]], 1, 1) end - @views @inbounds mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) + @views mul!(Cij, G[ba[i], 1:bb[j][end]], B[1:bb[j][end], bb[j]], -1, 1) - @inbounds _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij + _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij if alg == Val(:lyap) && i > j - @inbounds for l=bb[j], k=ba[i] # Avoids aliasing of copyto! + for l=bb[j], k=ba[i] # Avoids aliasing of copyto! C[l,k] = conj(C[k,l]) end end From f5bd0f8966716ca9cfb3cd1cc3c657d9176c5c1d Mon Sep 17 00:00:00 2001 From: olof3 Date: Thu, 14 May 2020 10:58:37 +0200 Subject: [PATCH 10/13] Updates to benchmark. --- lyap_bench.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/lyap_bench.jl b/lyap_bench.jl index ef2480226..b3f2c52ef 100644 --- a/lyap_bench.jl +++ b/lyap_bench.jl @@ -3,6 +3,9 @@ include("src/LyapTest.jl") BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5 +using Random +Random.seed!(0) + bg = BenchmarkGroup() for T in [Float64, ComplexF64] for k=[5,50,500]#[5, 20, 200] @@ -44,7 +47,7 @@ for k=1:4 table_content[1:3,1] = [x[1][4] for x in sorted_results[1:3]] table_content[1:3,2:4] = [string(BenchmarkTools.prettytime(median(x[2]).time), " (", x[2].allocs, " allocs, ", BenchmarkTools.prettymemory(x[2].memory), ")") for x in data] - pretty_table(table_content, ["", "LyapTest", "MatrixEquations", "julia"]) + pretty_table(table_content, ["N", "LyapTest", "MatrixEquations", "LinearAlgebra"]) end for k=1:4 @@ -56,7 +59,7 @@ for k=1:4 table_content[1:3,1] = [x[1][4] for x in sorted_results[1:3]] table_content[1:3,2:3] = [string(BenchmarkTools.prettytime(median(x[2]).time), " (", x[2].allocs, " allocs, ", BenchmarkTools.prettymemory(x[2].memory), ")") for x in data] - pretty_table(table_content, ["", "LyapTest", "MatrixEquations"]) + pretty_table(table_content, ["N", "LyapTest", "MatrixEquations"]) end From 6615224488802cc3b5ec42df3c2ddfa9839e8f0c Mon Sep 17 00:00:00 2001 From: olof3 Date: Thu, 14 May 2020 16:03:22 +0200 Subject: [PATCH 11/13] Fixes to docs. --- src/LyapTest.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index f8957f57e..5e3abdab3 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -84,13 +84,13 @@ sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C)) sylvdsoltype(A, B, C) = Base.promote_op(sylvd, eltype(A), eltype(B), eltype(C)) """ - _sylvc!(A, B, C) + _sylvc!(A, B, C) -> X Find the solution `X` to the continuous-time Sylvester equation `AX + XB = C` -for small matrices (1x1, 1x2, 2x1, 2x2) +for small matrices (1x1, 1x2, 2x1, 2x2), overwriting the input `C`. """ @inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) @@ -121,13 +121,13 @@ end """ - _sylvd!(X, A, B, C) + _sylvd!(A, B, C) -> X Find the solution `X` to the discrete-time Sylvester equation `AXB - X = C` -for small matrices (1x1, 1x2, 2x1, 2x2). +for small matrices (1x1, 1x2, 2x1, 2x2), overwriting the input `C`. """ function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) @@ -157,7 +157,7 @@ function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{ end """ - sylvc(A, B, C) + sylvc(A, B, C) -> X Find the solution `X` to the continuous-time Sylvester equation @@ -189,7 +189,7 @@ end end """ - sylvd(A, B, C) + sylvd(A, B, C) -> X Find the solution `X` to the discrete-time Sylvester equation @@ -221,7 +221,7 @@ end end """ - lyapc(A, Q) + lyapc(A, Q) -> X Computes the solution `X` of the continuous-time Lyapunov equation @@ -245,7 +245,7 @@ function lyapc(A, Q) X = mul!(Y, U, Y*U') end """ - `X = lyapd(A, Q)` + X = lyapd(A, Q) -> X Find the solution `X` to the discrete-time Lyapunov equation @@ -276,6 +276,8 @@ end """ + sylvc_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) -> X + Find the solution `X` to the continuous-time Sylvester equation `AX + XB = C` @@ -347,7 +349,7 @@ end """ - sylvd_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + sylvd_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) -> X Solve the discrete-time Sylvester equation From 7781ff9c06abe317ba7729a9028bf17faea9c019 Mon Sep 17 00:00:00 2001 From: olof3 Date: Fri, 15 May 2020 11:30:30 +0200 Subject: [PATCH 12/13] Minor fixes. --- src/LyapTest.jl | 56 +++++++++++++++++------------------ test/test_matrix_equations.jl | 2 +- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 5e3abdab3..4f4c04fca 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -46,7 +46,7 @@ function _schurstructure(R::AbstractMatrix, ul=Val(:U)::Union{Val{:U}, Val{:L}}) if j == n d[k] = 1 else - if ul == Val(:U) + if ul === Val(:U) d[k] = iszero(R[j+1, j]) ? 1 : 2 else d[k] = iszero(R[j, j+1]) ? 1 : 2 @@ -95,24 +95,24 @@ for small matrices (1x1, 1x2, 2x1, 2x2), overwriting the input `C`. @inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) if M == 2 && N == 2 - _sylvc!(A, B, C, Val{2}(), Val{2}()) + _sylvc!(A, B, C, Val(2), Val(2)) elseif M == 2 && N == 1 - _sylvc!(A, B, C, Val{2}(), Val{1}()) + _sylvc!(A, B, C, Val(2), Val(1)) elseif M == 1 && N == 2 - _sylvc!(A, B, C, Val{1}(), Val{2}()) + _sylvc!(A, B, C, Val(1), Val(2)) elseif M == 1 && N == 1 - _sylvc!(A, B, C, Val{1}(), Val{1}()) + _sylvc!(A, B, C, Val(1), Val(1)) else error("Matrix dimensionsins should not be greater than 2") end return C end @inline function _sylvc!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} - A′ = SMatrix{M,M}(A) - B′ = SMatrix{N,N}(B) - Cv = SMatrix{M,N}(C)[:] # vectorization of C + As = SMatrix{M,M}(A) + Bs = SMatrix{N,N}(B) + Cvs = SMatrix{M,N}(C)[:] # vectorization of C - Xv = lu(kron(SMatrix{N,N}(I), A′) + kron(transpose(B′), SMatrix{M,M}(I))) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) (with A = I or B = I) + Xv = lu(kron(SMatrix{N,N}(I), As) + kron(transpose(Bs), SMatrix{M,M}(I))) \ Cvs # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) (with A = I or B = I) if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvc or ?lyapc"); end @@ -132,24 +132,24 @@ for small matrices (1x1, 1x2, 2x1, 2x2), overwriting the input `C`. function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) M, N = size(C) if M == 2 && N == 2 - _sylvd!(A, B, C, Val{2}(), Val{2}()) + _sylvd!(A, B, C, Val(2), Val(2)) elseif M == 2 && N == 1 - _sylvd!(A, B, C, Val{2}(), Val{1}()) + _sylvd!(A, B, C, Val(2), Val(1)) elseif M == 1 && N == 2 - _sylvd!(A, B, C, Val{1}(), Val{2}()) + _sylvd!(A, B, C, Val(1), Val(2)) elseif M == 1 && N == 1 - _sylvd!(A, B, C, Val{1}(), Val{1}()) + _sylvd!(A, B, C, Val(1), Val(1)) else error("Matrix dimensionsins should not be greater than 2") end return C end function _sylvd!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, ::Val{M}, ::Val{N}) where {T <: Number, M, N} - A′ = SMatrix{M,M}(A) - B′ = SMatrix{N,N}(B) - Cv = SMatrix{M,N}(C)[:] # vectorization of C + As = SMatrix{M,M}(A) + Bs = SMatrix{N,N}(B) + Cvs = SMatrix{M,N}(C)[:] # vectorization of C - Xv = (kron(transpose(B′), A′) - I) \ Cv # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) + Xv = (kron(transpose(Bs), As) - SMatrix{M*N,M*N}(I)) \ Cvs # using the vectorization identity vec(AXB) = kron(B'*A)*vec(X) if any(!isfinite, Xv); error("Matrix equation has no solution, see ?sylvd or ?lyapd"); end @@ -302,11 +302,11 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va # The user should preferably use sylvc_schur! and lyapc_schur! # I.e., this method does not check whether C is hermitian # The matrix C is successively replaced with the solution X - # if alg == Val(:lyap), only the lower triangle of C is computed + # if alg === Val(:lyap), only the lower triangle of C is computed # after which an Hermitian view is applied # get block indices and nbr of blocks - if schurtype == Val(:real) + if schurtype === Val(:real) _, ba, nblocksa = _schurstructure(A, Val(:L)) # A is assumed upper triangualar _, bb, nblocksb = _schurstructure(B, Val(:U)) else @@ -315,15 +315,15 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va end @inbounds for j=1:nblocksb - i0 = (alg == Val(:lyap) ? j : 1) + i0 = (alg === Val(:lyap) ? j : 1) for i=i0:nblocksa - if schurtype == Val(:complex) + if schurtype === Val(:complex) if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end if j > 1; C[i,j] -= sum(C[i, k] * B[k, j] for k=1:j-1); end C[i,j] = sylvc(A[i, i], B[j, j], C[i, j]) # C[i,j] now contains solution Y[i,j] - if alg == Val(:lyap) && i > j + if alg === Val(:lyap) && i > j C[j,i] = conj(C[i,j]) end else @@ -336,7 +336,7 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va _sylvc!(Aii, Bjj, Cij) # Cij now contains the solution Yij - if alg == Val(:lyap) && i > j + if alg === Val(:lyap) && i > j for l=bb[j], k=ba[i] C[l,k] = conj(C[k,l]) end @@ -378,7 +378,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va G = zeros(eltype(C), size(A,1), size(B, 1)) # Keep track of A*X for improved performance # get block dimensions, block indices, nbr of blocks - if schurtype == Val(:real) + if schurtype === Val(:real) _, ba, nblocksa = _schurstructure(A, Val(:L)) # A is assumed upper triangualar _, bb, nblocksb = _schurstructure(B, Val(:U)) else @@ -387,9 +387,9 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va end @inbounds for j=1:nblocksb - i0 = (alg == Val(:lyap) ? j : 1) + i0 = (alg === Val(:lyap) ? j : 1) for i=i0:nblocksa - if schurtype == Val(:complex) + if schurtype === Val(:complex) # Compute Gij up to the contribution from Aii*Yij which is added at the end of each iteration if i > 1; G[i,j] += sum(A[i,k] * C[k,j] for k=1:i-1); end @@ -397,7 +397,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va C[i,j] = sylvd(A[i,i], B[j,j], C[i,j]) # C[i,j] now contains solution Y[i,j] - if alg == Val(:lyap) && i > j + if alg === Val(:lyap) && i > j C[j,i] = conj(C[i,j]) end @@ -417,7 +417,7 @@ function _sylvd_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va _sylvd!(Aii, Bjj, Cij) # Cij now contains the solution Yij - if alg == Val(:lyap) && i > j + if alg === Val(:lyap) && i > j for l=bb[j], k=ba[i] # Avoids aliasing of copyto! C[l,k] = conj(C[k,l]) end diff --git a/test/test_matrix_equations.jl b/test/test_matrix_equations.jl index 679a3a6df..09974d906 100644 --- a/test/test_matrix_equations.jl +++ b/test/test_matrix_equations.jl @@ -8,7 +8,7 @@ using MatrixEquations Random.seed!(0) -include("src/LyapTest.jl") +include("../src/LyapTest.jl") # Test _schurstructure R = diagm(0 => [1, 1, 1, 1, 1, 1, 1, 1], -1 => [1, 0, 0, 0, 1, 0, 0]) From eb4e945b0f2b817d2e97f2fae27765c04428a29b Mon Sep 17 00:00:00 2001 From: olof3 Date: Tue, 2 Jun 2020 14:42:50 +0200 Subject: [PATCH 13/13] Added naive method for generic inputs. --- src/LyapTest.jl | 62 ++++++++++++++++++++++++++++++++++- test/test_matrix_equations.jl | 3 ++ 2 files changed, 64 insertions(+), 1 deletion(-) diff --git a/src/LyapTest.jl b/src/LyapTest.jl index 4f4c04fca..44659c0b9 100644 --- a/src/LyapTest.jl +++ b/src/LyapTest.jl @@ -80,7 +80,8 @@ end LinearAlgebra.schur(A::AbstractMatrix{T}) where T = schur!(LinearAlgebra.copy_oftype(A, LinearAlgebra.eigtype(T))) -sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C)) +sylvcsoltype(A, B, C) = Base.promote_op((a,b,c) -> c / (a + b), eltype(A), eltype(B), eltype(C)) +#sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C)) sylvdsoltype(A, B, C) = Base.promote_op(sylvd, eltype(A), eltype(B), eltype(C)) """ @@ -236,6 +237,10 @@ function lyapc(A, Q) _check_lyap_inputs(A, Q) + if !hasmethod(schur!, (typeof(A),)) + return lyapc(A, Q, Val(:naive)) + end + At2, U = schur(A') Q2 = U'*Q*U @@ -244,6 +249,57 @@ function lyapc(A, Q) X = mul!(Y, U, Y*U') end + + +function sylvc(A, B, C, ::Val{:naive}) + Xv = kron(I(size(B,1)), A) + kron(transpose(B), I(size(A,1))) \ C[:] + return reshape(Xv, size(C)) +end + +# Mapping from Cartesian index into vector represenentation of Symmetric matrix +@inline sub2triidx(i,j) = (j >= i ? (j*(j-1))>>>1 + i : (i*(i-1))>>>1 + j) + +function lyapc(A, Q, ::Val{:naive}) # Only works for real matrices A + # Sets up and solves a system of equations for the upper triangular part of X + # and solves that equation. This gives an n(n+1)/2 system instead of n^2 + # ONLY WORKS FOR REAL A! + # Should be able to to base the discrete time version on this as well + _check_lyap_inputs(A, Q) + + if !isreal(A); error("Only works for real A matrices"); end # Should call sylvc + + n = size(Q, 1) + nR = (n*(n+1)) >>> 1 # Ssize of the system to be solved + + R = zeros(eltype(A), nR, nR) + minusQv = Vector{eltype(Q)}(undef, nR) + + for j=1:n + for i=1:j + m = sub2triidx(i,j) # Set up equation for X[i,j] + minusQv[m] = -Q[i,j] + # (AX + XA')[i,j] = sum(A[i,k]*X[k,j]) + sum(X[i,k]*A'[k,j]) + # the r = trinum[i]+r gives the kth element in the upper traingle + # which correpsonds to X[i,j] + for k=1:n + R[m, sub2triidx(k,j)] += A[i,k] + R[m, sub2triidx(i,k)] += A[j,k] + end + end + end + + Xv = R \ minusQv + + # Fill the upper traingle of X from Xv + X = [Xv[sub2triidx(i,j)] for i=1:n, j=1:n] + + return X +end + + + + + """ X = lyapd(A, Q) -> X @@ -316,6 +372,10 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va @inbounds for j=1:nblocksb i0 = (alg === Val(:lyap) ? j : 1) + + # if j > 1; mul!(C[i0:nblocksa,j], C[i0:nblocksa, 1:j-1], B[1:j-1, j], 1, -1); end # Could move this out? + # figure out the row indexing + for i=i0:nblocksa if schurtype === Val(:complex) if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end diff --git a/test/test_matrix_equations.jl b/test/test_matrix_equations.jl index 09974d906..fedb5c57e 100644 --- a/test/test_matrix_equations.jl +++ b/test/test_matrix_equations.jl @@ -76,6 +76,9 @@ Cherm = [1.0 2+im; 2-im 1] @test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, a, -C) @test LyapTest.lyapd(Ac, Cherm) ≈ sylvd_diag(ac, conj(ac), -Cherm) +# A few tests for the naive version, should have more +@test LyapTest.sylvc(A, B, C, Val(:naive)) ≈ sylvc_diag(a, a, -C) +@test LyapTest.lyapc(A, C, Val(:naive)) ≈ sylvc_diag(a, a, -C)