diff --git a/lyap_bench.jl b/lyap_bench.jl
new file mode 100644
index 000000000..b3f2c52ef
--- /dev/null
+++ b/lyap_bench.jl
@@ -0,0 +1,153 @@
+using MatrixEquations
+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]
+        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]))
+
+for k=1:4
+    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))
+    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, ["N", "LyapTest", "MatrixEquations", "LinearAlgebra"])
+end
+
+for k=1:4
+    data = reshape(sorted_results[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, ["N", "LyapTest", "MatrixEquations"])
+end
+
+
+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
new file mode 100644
index 000000000..44659c0b9
--- /dev/null
+++ b/src/LyapTest.jl
@@ -0,0 +1,495 @@
+module LyapTest
+
+#=
+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 after each block solve or at the end of the algorithm?
+=#
+
+using LinearAlgebra
+using StaticArrays
+
+
+
+"""
+    `_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.
+
+
+`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
+
+# 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
+
+# Should preferably be fixed in LinearAlgebra
+LinearAlgebra.schur(A::AbstractMatrix{T}) where T = schur!(LinearAlgebra.copy_oftype(A, LinearAlgebra.eigtype(T)))
+
+
+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))
+
+"""
+    _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), 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))
+    elseif M == 2 && N == 1
+        _sylvc!(A, B, C, Val(2), Val(1))
+    elseif M == 1 && N == 2
+        _sylvc!(A, B, C, Val(1), Val(2))
+    elseif M == 1 && N == 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}
+    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), 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
+
+    C .= reshape(Xv, M, N)
+end
+
+
+"""
+    _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), 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))
+    elseif M == 2 && N == 1
+        _sylvd!(A, B, C, Val(2), Val(1))
+    elseif M == 1 && N == 2
+        _sylvd!(A, B, C, Val(1), Val(2))
+    elseif M == 1 && N == 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}
+    As = SMatrix{M,M}(A)
+    Bs = SMatrix{N,N}(B)
+    Cvs = SMatrix{M,N}(C)[:] # vectorization of C
+
+    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
+
+    C .= reshape(Xv, M, N)
+end
+
+"""
+    sylvc(A, B, C) -> X
+
+Find the solution `X` to the continuous-time Sylvester equation
+
+`AX + XB = C`
+
+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.
+"""
+function sylvc(A, B, C)
+    _check_sylv_inputs(A, B, C)
+
+    At2, UA = schur(A')
+    B2, UB = schur(B)
+
+    C2 = UA'*C*UB # This should give C2 the right type
+
+    Y = _sylvc_schur!(Matrix(At2'), B2, C2, Val(:sylv))
+
+    X = mul!(Y, 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
+
+"""
+    sylvd(A, B, C) -> X
+
+Find the solution `X` to the discrete-time Sylvester equation
+
+`AXB - X = C`
+
+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.
+"""
+function sylvd(A, B, C)
+    _check_sylv_inputs(A, B, C)
+
+    At2, UA = schur(A')
+    B2, UB = schur(B)
+
+    C2 = UA'*C*UB
+
+    Y = _sylvd_schur!(Matrix(At2'), B2, C2, Val(:sylv))
+
+    X = mul!(Y, 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
+
+"""
+    lyapc(A, Q) -> X
+
+Computes the solution `X` of the continuous-time Lyapunov equation
+
+`AX + XA' + Q = 0`
+
+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)
+
+     _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
+
+    Y = _sylvc_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap))
+
+    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
+
+Find the solution `X` to the discrete-time Lyapunov equation
+
+`AXA' - X + Q = 0`
+
+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"
+    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)
+
+    _check_lyap_inputs(A, Q)
+
+    At2, U = schur(A')
+
+    Q2 = U'*Q*U
+
+    Y = _sylvd_schur!(Matrix(At2'), At2, lmul!(-1, Q2), Val(:lyap))
+
+    X = mul!(Y, U, Y*U')
+end
+
+
+"""
+    sylvc_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) -> X
+
+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)
+`B` is assumed to have upper Schur form
+
+See also `sylvc`
+"""
+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))
+    else
+        nblocksa = size(A, 1)
+        nblocksb = size(B, 1)
+    end
+
+    @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
+                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
+                    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; @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
+
+                _sylvc!(Aii, Bjj, Cij) # Cij now contains the solution Yij
+
+                if alg === Val(:lyap) && i > j
+                    for l=bb[j], k=ba[i]
+                        C[l,k] = conj(C[k,l])
+                    end
+                end
+            end
+        end
+    end
+    return C
+end
+
+
+"""
+    sylvd_schur!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) -> X
+
+Solve the discrete-time Sylvester equation
+
+`AXB - X = C`
+
+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 the matrix `C` has the right type, it is overwritten with the solution `X`.
+
+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))
+    else
+        nblocksa = size(A, 1)
+        nblocksb = size(B, 1)
+    end
+
+    @inbounds for j=1:nblocksb
+        i0 = (alg === Val(:lyap) ? j : 1)
+        for i=i0:nblocksa
+            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)
+
+                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
+                    C[j,i] = conj(C[i,j])
+                end
+
+                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])
+
+                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
+
+                @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
+                    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)
+            end
+        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..fedb5c57e
--- /dev/null
+++ b/test/test_matrix_equations.jl
@@ -0,0 +1,174 @@
+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 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))
+
+        @test Matrix(A*X*B - X - C) ≈ zeros(m,n) atol=1e-13
+
+        X = LyapTest._sylvc!(A, B, copy(C))
+
+        @test Matrix(A*X + X*B - C) ≈ zeros(m,n) atol=1e-13
+    end
+end
+
+
+# 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)
+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)) ≈ 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.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.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)
+
+# 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)
+
+
+
+#
+# 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
+
+
+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)
+
+        rtol = m > 50 ? 1e-10 : 1e-12
+
+        @testset "lyap/sylv ($m, $n, $TA, $TB, $TC)" begin
+            X =  LyapTest.sylvc(A, B, C)
+            @test A*X + X*B ≈ C rtol=rtol
+
+            X =  LyapTest.sylvd(A, B, C)
+            @test A*X*B - X ≈ C rtol=rtol
+
+            if m == n
+                Q = Matrix(Hermitian(randn(m,m)))
+
+                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