RFC: Lyapunov / Sylvester solver via LAPack xTRSYL

base/exports.jl

@@ -654,6 +654,7 @@ export
     lu,
     lufact!,
     lufact,
+    lyap,
     norm,
     null,
     peakflops,
@@ -677,6 +678,7 @@ export
     svdfact,
     svdvals!,
     svdvals,
+    sylvester,
     trace,
     transpose,
     tril!,

base/linalg.jl

@@ -85,6 +85,7 @@ export
     lu,
     lufact,
     lufact!,
+    lyap,
     norm,
     null,
     peakflops,
@@ -108,6 +109,7 @@ export
     svdfact,
     svdvals!,
     svdvals,
+    sylvester,
     trace,
     transpose,
     tril,

base/linalg/dense.jl

@@ -457,3 +457,28 @@ function cond(A::StridedMatrix, p::Real=2)
     end
     throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2 or Inf"))
 end
+
+## Lyapunov and Sylvester equation
+
+# AX + XB + C = 0
+function sylvester{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T})
+    RA, QA = schur(A)
+    RB, QB = schur(B)
+
+    D = -Ac_mul_B(QA,C*QB)
+    Y, scale = LAPACK.trsyl!('N','N', RA, RB, D)
+    scale!(QA*A_mul_Bc(Y,QB), inv(scale))
+end
+sylvester{T<:Integer}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) = sylvester(float(A), float(B), float(C))
+
+# AX + XA' + C = 0
+function lyap{T<:BlasFloat}(A::StridedMatrix{T},C::StridedMatrix{T})
+    R, Q = schur(A)
+
+    D = -Ac_mul_B(Q,C*Q)
+    Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
+    scale!(Q*A_mul_Bc(Y,Q), inv(scale))
+end
+lyap{T<:Integer}(A::StridedMatrix{T},C::StridedMatrix{T}) = lyap(float(A), float(C))
+lyap{T<:Number}(a::T, c::T) = -c/(2a)
+

base/linalg/lapack.jl

@@ -3665,4 +3665,36 @@ for (fn, elty) in ((:dtrttf_, :Float64),
     end
 end
 
+# Solves the real Sylvester matrix equation: op(A)*X +- X*op(B) = scale*C and A and B are both upper quasi triangular.
+for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64),
+                   (:strsyl_, :Float32, :Float32),
+                   (:ztrsyl_, :Complex128, :Float64),
+                   (:ctrsyl_, :Complex64, :Float32))
+    @eval begin
+        function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::Int=1)
+            chkstride1(A, B, C)
+            m, n = chksquare(A, B)
+            lda = max(1, stride(A, 2))
+            ldb = max(1, stride(B, 2))
+            m1, n1 = size(C)
+            if m != m1 || n != n1 throw(DimensionMismatch("")) end
+            ldc = max(1, stride(C, 2))
+
+            scale = Array($relty, 1)
+            info = Array(BlasInt, 1)
+
+            ccall(($(string(fn)), liblapack), Void, 
+                (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
+                 Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
+                 Ptr{$relty}, Ptr{BlasInt}),
+                &transa, &transb, &isgn, &m, &n,
+                A, &lda, B, &ldb, C, &ldc,
+                scale, info)
+            @lapackerror 
+            C, scale[1]
+        end
+    end
+end
+
+
 end # module

doc/stdlib/linalg.rst

@@ -463,6 +463,14 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Matrix exponential.
 
+.. function:: lyap(A, C)
+
+   Computes the solution ``X`` to the continuous Lyapunov equation ``AX + XA' + C = 0``, where no eigenvalue of ``A`` has a zero real part and no two eigenvalues are negative complex conjugates of each other. 
+
+.. function:: sylvester(A, B, C)
+
+   Computes the solution ``X`` to the Sylvester equation ``AX + XB + C = 0``, where ``A``, ``B`` and ``C`` have compatible dimensions and ``A`` and ``-B`` have no eigenvalues with equal real part.
+
 .. function:: issym(A) -> Bool
 
    Test whether a matrix is symmetric.

test/linalg1.jl

@@ -2,7 +2,7 @@ debug = false
 
 import Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted
 
-n     = 10
+n = 10
 srand(1234321)
 
 a = rand(n,n)
@@ -17,16 +17,20 @@ end
 
 areal = randn(n,n)/2
 aimg  = randn(n,n)/2
+a2real = randn(n,n)/2
+a2img  = randn(n,n)/2
 breal = randn(n,2)/2
 bimg  = randn(n,2)/2
+
 for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
+    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
+    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real)
+    asym = a'+a                  # symmetric indefinite
+    apd  = a'*a                 # symmetric positive-definite
+    ε = εa = eps(abs(float(one(eltya))))
+
     for eltyb in (Float32, Float64, Complex64, Complex128, Int)
-        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
         b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
-        asym = a'+a                  # symmetric indefinite
-        apd  = a'*a                 # symmetric positive-definite
-
-        εa = eps(abs(float(one(eltya))))
         εb = eps(abs(float(one(eltyb))))
         ε = max(εa,εb)
 
@@ -249,23 +253,36 @@ debug && println("Test null")
         @test size(null(b), 2) == 0
     end
 
+    end # for eltyb
+
+debug && println("\ntype of a: ", eltya, "\n")
+
 debug && println("Test pinv")
-    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+    if eltya != BigFloat # Revisit when implemented in julia
         pinva15 = pinv(a[:,1:5])
         @test_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
         @test_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
     end
 
     # if isreal(a)
 debug && println("Matrix square root")
-    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+    if eltya != BigFloat # Revisit when implemented in julia
         asq = sqrtm(a)
         @test_approx_eq asq*asq a
         asymsq = sqrtm(asym)
         @test_approx_eq asymsq*asymsq asym
     end
-end
-end
+
+debug && println("Lyapunov/Sylvester")
+    if eltya != BigFloat
+        let 
+            x = lyap(a, a2)
+            @test_approx_eq -a2 a*x + x*a'
+            x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n])
+            @test_approx_eq -a2[1:3, 4:n] a[1:3, 1:3]*x2 + x2*a[4:n, 4:n]
+        end
+    end
+end # for eltya
 
 #6941
 #@test (ones(10^7,4)*ones(4))[3] == 4.0