Generalize lyap and sylvester to abstract matrices (#46545)

Author: dkarrasch

stdlib/LinearAlgebra/src/dense.jl

@@ -1580,21 +1580,22 @@ julia> X = sylvester(A, B, C)
  -4.46667   1.93333
   3.73333  -1.8
 
-julia> A*X + X*B + C
-2×2 Matrix{Float64}:
-  2.66454e-15  1.77636e-15
- -3.77476e-15  4.44089e-16
+julia> A*X + X*B ≈ -C
+true
 ```
 """
-function sylvester(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) where T<:BlasFloat
+function sylvester(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix)
+    T = promote_type(float(eltype(A)), float(eltype(B)), float(eltype(C)))
+    return sylvester(copy_similar(A, T), copy_similar(B, T), copy_similar(C, T))
+end
+function sylvester(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}) where {T<:BlasFloat}
     RA, QA = schur(A)
     RB, QB = schur(B)
-
-    D = -(adjoint(QA) * (C*QB))
-    Y, scale = LAPACK.trsyl!('N','N', RA, RB, D)
-    rmul!(QA*(Y * adjoint(QB)), inv(scale))
+    D = QA' * C * QB
+    D .= .-D
+    Y, scale = LAPACK.trsyl!('N', 'N', RA, RB, D)
+    rmul!(QA * Y * QB', inv(scale))
 end
-sylvester(A::StridedMatrix{T}, B::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:Integer} = sylvester(float(A), float(B), float(C))
 
 Base.@propagate_inbounds function _sylvester_2x1!(A, B, C)
     b = B[1]
@@ -1652,18 +1653,19 @@ julia> X = lyap(A, B)
   0.5  -0.5
  -0.5   0.25
 
-julia> A*X + X*A' + B
-2×2 Matrix{Float64}:
- 0.0          6.66134e-16
- 6.66134e-16  8.88178e-16
+julia> A*X + X*A' ≈ -B
+true
 ```
 """
-function lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:BlasFloat}
+function lyap(A::AbstractMatrix, C::AbstractMatrix)
+    T = promote_type(float(eltype(A)), float(eltype(C)))
+    return lyap(copy_similar(A, T), copy_similar(C, T))
+end
+function lyap(A::AbstractMatrix{T}, C::AbstractMatrix{T}) where {T<:BlasFloat}
     R, Q = schur(A)
-
-    D = -(adjoint(Q) * (C*Q))
+    D = Q' * C * Q
+    D .= .-D
     Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
-    rmul!(Q*(Y * adjoint(Q)), inv(scale))
+    rmul!(Q * Y * Q', inv(scale))
 end
-lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:Integer} = lyap(float(A), float(C))
 lyap(a::Union{Real,Complex}, c::Union{Real,Complex}) = -c/(2real(a))

stdlib/LinearAlgebra/test/dense.jl

@@ -132,8 +132,20 @@ bimg  = randn(n,2)/2
         @testset "Lyapunov/Sylvester" begin
             x = lyap(a, a2)
             @test -a2 ≈ a*x + x*a'
+            y = lyap(a', a2')
+            @test y ≈ lyap(Array(a'), Array(a2'))
+            @test -a2' ≈ a'y + y*a
+            z = lyap(Tridiagonal(a)', Diagonal(a2))
+            @test z ≈ lyap(Array(Tridiagonal(a)'), Array(Diagonal(a2)))
+            @test -Diagonal(a2) ≈ Tridiagonal(a)'*z + z*Tridiagonal(a)
             x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n])
             @test -a2[1:3, 4:n] ≈ a[1:3, 1:3]*x2 + x2*a[4:n, 4:n]
+            y2 = sylvester(a[1:3, 1:3]', a[4:n, 4:n]', a2[4:n,1:3]')
+            @test y2 ≈ sylvester(Array(a[1:3, 1:3]'), Array(a[4:n, 4:n]'), Array(a2[4:n,1:3]'))
+            @test -a2[4:n, 1:3]' ≈ a[1:3, 1:3]'*y2 + y2*a[4:n, 4:n]'
+            z2 = sylvester(Tridiagonal(a[1:3, 1:3]), Diagonal(a[4:n, 4:n]), a2[1:3,4:n])
+            @test z2 ≈ sylvester(Array(Tridiagonal(a[1:3, 1:3])), Array(Diagonal(a[4:n, 4:n])), Array(a2[1:3,4:n]))
+            @test -a2[1:3, 4:n] ≈ Tridiagonal(a[1:3, 1:3])*z2 + z2*Diagonal(a[4:n, 4:n])
         end
 
         @testset "Matrix square root" begin