Lyapunov equation and xtrsyl wrapper

lyap.jl

using Base.LinAlg.LAPACK
import Base.LinAlg: BlasFloat, BlasChar, BlasInt, LAPACK.liblapack, LAPACK.@assertargsok, LAPACK.chkstride1

for (fn, elty) in ((:dtrsyl_, :Float64),
                   (:strsyl_, :Float32),
                   (:ztrsyl_, :Complex128),
                   (:ctrsyl_, :Complex64))
    @eval begin
        function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::BlasInt=1)

            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            m = size(A, 1)
            n = size(B, 1)
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            if lda < m throw(DimensionMismatch("")) end
            if ldb < n throw(DimensionMismatch("")) end
            m1, n1 = size(C)
            if m != m1 || n != n1 throw(DimensionMismatch("")) end
            ldc = max(1, stride(C, 2))

            scale = Array($elty, 1)
            info = Array(BlasInt, 1)
            
            #  SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO )
            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{$elty}, Ptr{BlasInt}),
                 &transa, &transb, &isgn, &m, &n,
                 A, &lda, B, &ldb, C, &ldc,
                 scale, info)
            
            info[1]>0 && throw(SingularException(info[1]))
            return C, scale[1]
        end
    end
end

# AX + XB' + C = 0
function syl{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T})
    RA, QA =schur(A)
    RB, QB =schur(B)

    D = -QA'*C*QB
    Y, scale = trsyl!('N', 'T', RA, RB, D)
    (QA*Y*QB')/scale 
end

# AX + XA' + C = 0
function lyap{T<:BlasFloat}(A::StridedMatrix{T},C::StridedMatrix{T})
    R, Q =schur(A)

    D = -Q'*C*Q
    Y, scale = trsyl!('N', 'T', R, R, D)
    (Q*Y*Q')/scale 
end
lyap(a::Float64, c::Float64) = -0.5c/a


function kronsyl(a, b, c)
    Base.LinAlg.chksquare(a, b, c)
    k = kron(eye(a), a) + kron(b', eye(b))
    xvec=k\vec(c)
    reshape(xvec,size(c))
end

kronlyap(b, a) = kronsyl(b, b', -a)


function stable(Y, d, ep)

    # convert first d*(d+1)/2 values of Y into upper triangular matrix
    # positive definite matrix
    x = zeros(d,d)
    k = 1
    for i in 1:d
        for j in i:d
            x[i,j] = Y[k]
            k = k + 1
        end
    end
    # convert next d*(d+1)/2 -d values of Y into anti symmetric matrix
    y = zeros(d,d)
    for i in 1:d
        for j in i+1:d
            y[i,j] = Y[k]
            y[j,i] = -y[i, j]
            k = k + 1
        end
    end
    assert(k -1 == d*d == length(Y))

    # return stable matrix as a sum of a antisymmetric and a positive definite matrix
    y - x'*x - ep*eye(d)
end

#% .. function:: randposdef(d)
#%
#% Random positive definite matrix of dimension ``d``.
#%

function randposdef(d)
    x = randn(d, d)
    x'* x/sqrt(d)
end

#% .. function:: randstable(d)
#%
#% Random stable matrix (matrix with eigenvalues with negative real part) with
#% dimension ``d``.

function randstable(d)
    # positive definite matrix
    x = randn(d, d)
    a = x'*x

    # anti symmetric matrix
    y = triu(randn(d, d))
    b = y - y'

    # return stable matrix
    b/sqrt(2) - a/sqrt(2d)
end