text:
export LaplacianOptModel

"""
    LaplacianOptModel

The composite mutable struct, `LaplacianOptModel`, holds dictionaries for input data, abstract JuMP model for optimization,
variable references and result from solving the JuMP model.
"""
mutable struct LaplacianOptModel
    data::Dict{String,Any}
    model::JuMP.Model
    variables::Dict{Symbol,Any}
    result::Dict{String,Any}

    "Contructor for struct `LaplacianOptModel`"
    function LaplacianOptModel(data::Dict{String,Any})

        data = data
        model = JuMP.Model()
        variables = Dict{Symbol,Any}()
        result = Dict{String,Any}()
        lom = new(data, model, variables, result)

        return lom
    end

end

"""
    GraphData

The composite mutable struct, `GraphData`, holds matrices of adjacency matrix, laplacian matrix,
fiedler vector and algebraic connectivity.
"""
mutable struct GraphData
    adjacency::Array{Float64}
    laplacian::Array{Float64}
    fiedler::Vector{Float64}
    ac::Number

    function GraphData(adjacency::Array{Float64})
        adj_matrix = adjacency
        laplacian = LOpt.laplacian_matrix(adjacency)
        fiedler = LOpt.fiedler_vector(adjacency)
        ac = LOpt.algebraic_connectivity(adjacency)

        graph_data = new(adj_matrix, laplacian, fiedler, ac)

        return graph_data

    end
end


text:
module EMLheater
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("heater/heater.jl")
	include("heater/cooler.jl")
end

ext: jl


text:
import DynamicEnergyBudgets: get_environment, temp

soiltemperature = [15.4u"°C"]
soilwatercontent = [0.7]
soilwaterpotential = [100.0u"Pa"]
airtemperature = [23.7u"°C"]
windspeed = [3.7u"m*s^-1"]
relhumidity = [0.74]
radiation = [985.0u"W*m^-2"]
par = [3130.0u"mol*m^-2*s^-1"]

environment = Env(soiltemperature, soilwatercontent, soilwaterpotential,
                  airtemperature, windspeed, relhumidity, radiation, par)

u = [9.0u"mol",8.0u"mol",7.0u"mol",6.0u"mol",5.0u"mol",4.0u"mol"]

@testset "apply environment to deb vars" begin
    global o = DynamicEnergyBudgets.Plant(environment=environment)
    global o1, o2 = define_organs(o, 1)
    global v1 = o1.vars; v2 = o2.vars
    global va1 = v1.assimilation
    global va2 = v2.assimilation

    apply(apply_environment!, (o1, o2), u, :not_nothing, 1)
    @test temp(v1) == 23.7u"°C"
    @test va1.J_L_F == 3130.0u"mol*m^-2*s^-1"

    @test temp(v2) == 15.4u"°C"
    @test upreferred(va2.X_H) ≈ 0.7 * 1.0u"m^3*m^-3" * 1u"kg*L^-1" / 18.0u"g*mol^-1"
end

@testset "apply environment to deb vars" begin
    global o = DynamicEnergyBudgets.FvCBPlant(environment=environment)
    global o1, o2 = define_organs(o, 1)
    global v = o1.vars;
    global va = v.assimilation;
    apply_environment!(o1, u, :no_env, 1)
    @test va.tair == 23.7u"°C"
    @test va.tleaf > va.tair
    @test temp(v) == va.tleaf
    @test va.windspeed == 3.7u"m*s^-1"
    @test va.rh == 0.74
    @test va.rnet == 985.0u"W*m^-2"
    @test va.par == 3130.0u"mol*m^-2*s^-1"
    @test va.soilmoist == 0.7
    @test va.swp == 100.0u"Pa"
end

nothing


text:
function evaluateSS(ssfunction::Array{Function,1},varvalues::Array{Float64,1},parametervalues::Array{Float64,1})
    out = Array{Float64}(length(ssfunction))
    for (i,f) in enumerate(ssfunction)
        out[i] = f(varvalues...,parametervalues...)
    end
    return out
end

function lambdifyss(ss::Array{SymPy.Sym,1},vars::Array{SymPy.Sym,1},parameters::Array{SymPy.Sym,1})
    # Takes in sympy expressions, and symbolic object versions of the variables and parameters
    # Converts to a function that takes in arrays of variable values and parameter values
    # and returns equation errors
    functionarray = Array{Function}(length(ss))
    for (i,expression) in enumerate(ss)
        functionarray[i] = SymPy.lambdify(expression,[vars;parameters])
    end
    return functionarray
end

function createSSsystem(numparameters,expressions,leads,lags,contemps,statesindices,shocks,
                        parameters,barversions,parameter_affectsSS,switching,parameter_names)

    parametersSS = createSSparameters(numparameters,parameter_affectsSS,barversions,parameter_names,parameters,switching)
    subsargs = createpairedtuplelist([leads; lags; shocks; reshape(parameters[1,:],numparameters); reshape(parameters[2,switching],countnz(switching))],
    [contemps; contemps[statesindices]; zeros(Float64,length(shocks)); parametersSS; parametersSS[switching]])
    expressionsSS = SymPy.subs(expressions,subsargs...)
    lambdifiedSS = lambdifyss(expressionsSS,contemps,parametersSS)
    return lambdifiedSS
end

function findparametersaffectSS(numparameters,parameter_names,parameter_sym,parameter_switching,
    parameter_values,expressions,leads,lags,contemps,statesindices,shocks,ss_guess,ergodic)

    println("Checking which of the switching parameters affect the steady state:")
    affectSS = BitArray(numparameters)
    fill!(affectSS,false)

    # generate bar versions for _all_ parameters
    baselinevalues = squeeze(permutedims(parameter_values,[2 1])*reshape(ergodic,(length(ergodic),1)),2)

    # Create and solve the SS using those baseline values
    SSsystem = createSSsystem(numparameters,expressions,leads,lags,contemps,statesindices,shocks,
    parameter_sym,Array{SymPy.Sym}(0),affectSS,parameter_switching,parameter_names)

    sswrapper = x->evaluateSS(SSsystem,x,baselinevalues)

    print("Solving for steady state...")
    initialeval = sswrapper(ss_guess)
    if any(isnan.(initialeval)) || any(initialeval.== Inf) || any(initialeval.== -Inf)
        error("Steady state guess does not evaluate (e.g. due to divide-by-zero). Supply a different initial guess.")
    end
    nlout = NLsolve.nlsolve(NLsolve.not_in_place(sswrapper), ss_guess)
    if !NLsolve.converged(nlout)
        println("FAILED.")
        println("Residuals at initial guess:")
        println(sswrapper(ss_guess))
        error("Could not find steady state.")
    end
    SSbaseline = nlout.zero
    println("done.")

    for (p,switches) in enumerate(parameter_switching)
        if switches
            print(string(parameter_names[p],"..."))
            values = parameter_values[:,p]
            # Find a value in the regime switching parameter values that is not equal to the baseline
            newindex = find(values.!=baselinevalues[p])
            if length(newindex) == 0
                println("does not.")
                warn("You declared $(parameter_names[p]) to be switching by giving it multiple values, but the values are all identical. This is very inefficient.")
            else
                newvalue = values[minimum(newindex)]
                newvalues = copy(baselinevalues)
                newvalues[p] = newvalue
                errors = evaluateSS(SSsystem,SSbaseline,newvalues)
                if any(abs.(errors).>nlout.ftol)
                    affectSS[p] = true
                    println("does.")
                else
                    println("does not.")
                end
            end
        else
            affectSS[p] = false
        end
    end

    return affectSS,SSbaseline
end

function SteadyState(meta,equations,vars,parameters,shocks,ss_guess)
    if length(ss_guess) != meta.numvars
        error("Your supplied guess for the steady state has $(length(ss_guess)) entries, but you have $(meta.numvars) variables.")
    end
    ss_system = createSSsystem(meta.numparameters,equations.syms,vars.leads_sym,vars.states_sym,vars.contemps_sym,
    vars.states_indices,shocks.syms,parameters.generic_sym,parameters.bar_sym,parameters.affectsSS,parameters.isswitching,parameters.names)

    return SteadyState(ss_system,ss_guess)
end


text:
using Franklin
if !(@isdefined msg)
    msg = "franklin-update"
end
publish(prerender=false, final=lunr, message=msg)

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
"""
    X = sylvc(A,B,C)

Solve the continuous Sylvester matrix equation

                AX + XB = C

using the Bartels-Stewart Schur form based approach. `A` and `B` are
square matrices, and `A` and `-B` must not have common eigenvalues.

The following particular cases are also adressed:

    X = sylvc(α*I,B,C)  or  X = sylvc(α,B,C)

Solve the matrix equation `X(αI+B)  = C`.

    X = sylvc(A,β*I,C)  or  X = sylvc(A,β,C)

Solve the matrix equation `(A+βI)X = C`.

    X = sylvc(α*I,β*I,C)  or  sylvc(α,β,C)

Solve the matrix equation `(α+β)X = C`.

    x = sylvc(α,β,γ)

Solve the equation `(α+β)x = γ`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0

julia> X = sylvc(A, B, C)
2×2 Array{Float64,2}:
 -4.46667   1.93333
  3.73333  -1.8

julia> A*X + X*B - C
2×2 Array{Float64,2}:
  2.66454e-15  1.77636e-15
 -3.77476e-15  4.44089e-16
```
"""
function sylvc(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix)
   """
   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """

   m, n = size(C);
   [m; n] == LinearAlgebra.checksquare(A,B) || throw(DimensionMismatch("A, B and C have incompatible dimensions"))

   T2 = promote_type(eltype(A), eltype(B), eltype(C))
   T2 <: BlasFloat || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (A = convert(Matrix{T2},A))
   eltype(B) == T2 || (B = convert(Matrix{T2},B))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   adjA = isa(A,Adjoint)
   adjB = isa(B,Adjoint)
   if adjA
      RA, QA = schur(A.parent)
   else
      RA, QA = schur(A)
   end
   if adjB
      RB, QB = schur(B.parent)
   else
      RB, QB = schur(B)
   end

   Y = QA' * C * QB

   sylvcs!(RA, RB, Y, adjA = adjA, adjB = adjB)

   mul!(Y, QA, Y*QB')

   return Y
end
# solve X(B+α) = C or (α+β)X = C
sylvc(A::Union{Real,Complex,UniformScaling},B::Union{AbstractMatrix,UniformScaling},C::AbstractMatrix) = C/(A*I+B)
# solve (A+β)X = C
sylvc(A::AbstractMatrix,B::Union{Real,Complex,UniformScaling},C::AbstractMatrix) = (A+B*I)\C
# solve (α+β)X = C
sylvc(A::Union{Real,Complex},B::Union{Real,Complex},C::AbstractMatrix) = A+B == 0 ? throw(SingularException(1)) : C/(A+B)
# solve (α+β)x = γ
sylvc(A::Union{Real,Complex}, B::Union{Real,Complex}, C::Union{Real,Complex}) = A+B == 0 ? throw(SingularException(1)) : C/(A+B)
"""
    X = sylvd(A,B,C)

Solve the discrete Sylvester matrix equation

                AXB + X = C

using an extension of the Bartels-Stewart Schur form based approach.
`A` and `B` are square matrices, and `A` and `-B` must not have
common reciprocal eigenvalues.

The following particular cases are also adressed:

    X = sylvd(α*I,B,C)  or  X = sylvd(α,B,C)

Solve the matrix equation `X(αB+I)  = C`.

    X = sylvd(A,β*I,C)   or  X = sylvd(A,β,C)

Solve the matrix equation `(βA+I)X = C`.

    X = sylvd(α*I,β*I,C)  or  X = sylvd(α,β,C)

Solve the matrix equation `(αβ+1)X = C`.

    x = sylvd(α,β,γ)

Solve the equation `(αβ+1)x = γ`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0

julia> X = sylvd(A, B, C)
2×2 Array{Float64,2}:
 -2.46667  -2.73333
  2.4       1.86667

julia> A*X*B + X - C
2×2 Array{Float64,2}:
  8.88178e-16   8.88178e-16
 -3.9968e-15   -5.55112e-15
```
"""
function sylvd(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix)
   """
   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """

   m, n = size(C);
   [m; n] == LinearAlgebra.checksquare(A,B) || throw(DimensionMismatch("A, B and C have incompatible dimensions"))
   T2 = promote_type(eltype(A), eltype(B), eltype(C))
   T2 <: BlasFloat || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (A = convert(Matrix{T2},A))
   eltype(B) == T2 || (B = convert(Matrix{T2},B))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   adjA = isa(A,Adjoint)
   adjB = isa(B,Adjoint)
   if adjA
      RA, QA = schur(A.parent)
   else
      RA, QA = schur(A)
   end
   if adjB
      RB, QB = schur(B.parent)
   else
      RB, QB = schur(B)
   end

   Y = QA' * C * QB

   sylvds!(RA, RB, Y, adjA = adjA, adjB = adjB)

   mul!(Y, QA, Y*QB')

   return Y
end
# solve X(αB+I) = C or X(αβ+1) = C
sylvd(A::Union{Real,Complex,UniformScaling},B::Union{AbstractMatrix,UniformScaling},C::AbstractMatrix) = C/(A*B+I)
# solve (Aβ+I)X = C
sylvd(A::AbstractMatrix,B::Union{Real,Complex,UniformScaling},C::AbstractMatrix) = (A*B+I)\C
# solve (αβ+1)X = C
sylvd(A::Union{Real,Complex},B::Union{Real,Complex},C::AbstractMatrix) = A*B+1 == 0 ? throw(SingularException(1)) : C/(A*B+1)
# solve (αβ+1)x = γ
sylvd(A::Union{Real,Complex}, B::Union{Real,Complex}, C::Union{Real,Complex}) = A*B+1 == 0 ? throw(SingularException(1)) : C/(A*B+one(C))
"""
    X = gsylv(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

              AXB + CXD = E

using a generalized Schur form based approach. `A`, `B`, `C` and `D` are
square matrices. The pencils `A-λC` and `D+λB` must be regular and
must not have common eigenvalues.

The following particular cases are also adressed:

    X = gsylv(A,B,E)

Solve the generalized Sylvester matrix equation `AXB  = E`.

    X = gsylv(A,B,γ*I,E)  or  X = gsylv(A,B,γ,E)

Solve the generalized Sylvester matrix equation `AXB +γX = E`.

    X = gsylv(A,B,γ*I,D,E)  or  X = gsylv(A,B,γ,D,E)

Solve the generalized Sylvester matrix equation `AXB +γXD = E`.

    X = gsylv(A,B,C,δ*I,E)  or  X = gsylv(A,B,C,δ,E)

Solve the generalized Sylvester matrix equation `AXB +CXδ = E`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0

julia> D = [1. -2.; -2. -1.]
2×2 Array{Float64,2}:
  1.0  -2.0
 -2.0  -1.0

julia> E = [1. -1.; -2. 2.]
2×2 Array{Float64,2}:
  1.0  -1.0
 -2.0   2.0

julia> X = gsylv(A, B, C, D, E)
2×2 Array{Float64,2}:
 -0.52094   -0.0275792
 -0.168539   0.314607

julia> A*X*B + C*X*D - E
2×2 Array{Float64,2}:
 4.44089e-16  8.88178e-16
 6.66134e-16  0.0
```
"""
function gsylv(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix,D::AbstractMatrix,E::AbstractMatrix)

    m, n = size(E);
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,C,D) ||
       throw(DimensionMismatch("A, B, C, D and E have incompatible dimensions"))
    T2 = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E))
    T2 <: BlasFloat || (T2 = promote_type(Float64,T2))
    eltype(A) == T2 || (A = convert(Matrix{T2},A))
    eltype(B) == T2 || (B = convert(Matrix{T2},B))
    eltype(C) == T2 || (C = convert(Matrix{T2},C))
    eltype(D) == T2 || (D = convert(Matrix{T2},D))
    eltype(E) == T2 || (E = convert(Matrix{T2},E))

    adjA = isa(A,Adjoint)
    adjB = isa(B,Adjoint)
    adjC = isa(C,Adjoint)
    adjD = isa(D,Adjoint)
    adjAC = adjA && adjC
    adjBD = adjB && adjD

    if adjAC
       AS, CS, Z1, Q1 = schur(A.parent,C.parent)
    else
       adjA && (A = copy(A))
       adjC && (C = copy(C))
       AS, CS, Q1, Z1 = schur(A,C)
    end
    if adjBD
       BS, DS, Z2, Q2 = schur(B.parent,D.parent)
    else
       adjB && (B = copy(B))
       adjD && (D = copy(D))
       BS, DS, Q2, Z2 = schur(B,D)
    end

    Y = Q1' * E *Z2

    gsylvs!(AS, BS, CS, DS, Y, adjAC = adjAC, adjBD = adjBD)

    mul!(Y, Z1, Y*Q2')

   return Y
end
# solve AXB = C
gsylv(A::Union{AbstractMatrix,UniformScaling,Real,Complex},B::Union{AbstractMatrix,UniformScaling,Real,Complex},E::AbstractMatrix) = (A\E)/B
# solve AXB+γX = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},E::AbstractMatrix) =
size(A,1) == size(A,2) && size(B,1) == size(B,2) ? gsylv(A,B,Matrix{eltype(C)}(C*I,size(A)),Matrix{eltype(C)}(I,size(B)),E) :
throw(DimensionMismatch("A and B must be square matrices"))
# solve AXB+γXδ = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},D::Union{UniformScaling,Real,Complex},E::AbstractMatrix) = gsylv(A,B,C*D,E)
# solve AXB+γXD = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},D::AbstractMatrix,E::AbstractMatrix) =
size(A,1) == size(A,2) ? gsylv(A,B,Matrix{eltype(C)}(C*I,size(A)),D,E) :
throw(DimensionMismatch("A must be a square matrix"))
# solve AXB+CXδ = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix,D::Union{UniformScaling,Real,Complex},E::AbstractMatrix) =
size(B,1) == size(B,2) ? gsylv(A,B,C,Matrix{eltype(D)}(D*I,size(B)),E) :
throw(DimensionMismatch("B must be a square matrix"))

"""
    (X,Y) = sylvsys(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F,

where `(A,D)`, `(B,E)` are pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.
# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0

julia> D = [1. -2.; -2. -1.]
2×2 Array{Float64,2}:
  1.0  -2.0
 -2.0  -1.0

julia> E = [1. -1.; -2. 2.]
2×2 Array{Float64,2}:
  1.0  -1.0
 -2.0   2.0

julia> F = [1. -1.; -2. 2.]
2×2 Array{Float64,2}:
  1.0  -1.0
 -2.0   2.0

julia> X, Y = sylvsys(A, B, C, D, E, F);

julia> X
2×2 Array{Float64,2}:
  1.388  -1.388
 -0.892   0.892

julia> Y
2×2 Array{Float64,2}:
 -1.788  0.192
  0.236  0.176

julia> A*X + Y*B - C
2×2 Array{Float64,2}:
  6.66134e-16  2.22045e-15
 -3.10862e-15  2.66454e-15

julia> D*X + Y*E - F
2×2 Array{Float64,2}:
  1.33227e-15  -2.22045e-15
 -4.44089e-16   4.44089e-16
```
"""
function sylvsys(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix,D::AbstractMatrix,E::AbstractMatrix,F::AbstractMatrix)

    m, n = size(C);
    (m == size(F,1) && n == size(F,2)) ||
      throw(DimensionMismatch("C and F must have the same dimensions"))
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,D,E) ||
       throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))
    T2 = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E), eltype(F))
    T2 <: BlasFloat || (T2 = promote_type(Float64,T2))
    eltype(A) == T2 || (A = convert(Matrix{T2},A))
    eltype(B) == T2 || (B = convert(Matrix{T2},B))
    eltype(C) == T2 || (C = convert(Matrix{T2},C))
    eltype(D) == T2 || (D = convert(Matrix{T2},D))
    eltype(E) == T2 || (E = convert(Matrix{T2},E))
    eltype(F) == T2 || (F = convert(Matrix{T2},F))

    isa(A,Adjoint) && (A = copy(A))
    isa(B,Adjoint) && (B = copy(B))
    isa(D,Adjoint) && (D = copy(D))
    isa(E,Adjoint) && (E = copy(E))

    AS, DS, Q1, Z1 = schur(A,D)
    BS, ES, Q2, Z2 = schur(B,E)

    CS = adjoint(Q1) * (C*Z2)
    FS = adjoint(Q1) * (F*Z2)

    X, Y, scale =  tgsyl!('N',AS,BS,CS,DS,ES,FS)

    (rmul!(Z1*(X * adjoint(Z2)), inv(scale)), rmul!(Q1*(Y * adjoint(Q2)), inv(-scale)) )
end
"""
    (X,Y) = dsylvsys(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

       AX + DY = C
       XB + YE = F ,

where `(A,D)`, `(B,E)` are pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0

julia> D = [1. -2.; -2. -1.]
2×2 Array{Float64,2}:
  1.0  -2.0
 -2.0  -1.0

julia> E = [1. -1.; -2. 2.]
2×2 Array{Float64,2}:
  1.0  -1.0
 -2.0   2.0

julia> F = [1. -1.; -2. 2.]
2×2 Array{Float64,2}:
  1.0  -1.0
 -2.0   2.0

julia> X, Y = dsylvsys(A, B, C, D, E, F);

julia> X
2×2 Array{Float64,2}:
  2.472  -1.648
 -1.848   1.232

julia> Y
2×2 Array{Float64,2}:
 -0.496  -0.336
  0.264   0.824

julia> A*X + D*Y - C
2×2 Array{Float64,2}:
  4.44089e-16  0.0
 -3.55271e-15  1.55431e-15

julia> X*B + Y*E - F
2×2 Array{Float64,2}:
 -8.88178e-16   0.0
  8.88178e-16  -4.44089e-16
```
"""
function dsylvsys(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix,D::AbstractMatrix,E::AbstractMatrix,F::AbstractMatrix)

    m, n = size(C);
    (m == size(F,1) && n == size(F,2)) ||
      throw(DimensionMismatch("C and F must have the same dimensions"))
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,D,E) ||
       throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))

    T2 = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E), eltype(F))
    T2 <: BlasFloat || (T2 = promote_type(Float64,T2))
    eltype(A) == T2 || (A = convert(Matrix{T2},A))
    eltype(B) == T2 || (B = convert(Matrix{T2},B))
    eltype(C) == T2 || (C = convert(Matrix{T2},C))
    eltype(D) == T2 || (D = convert(Matrix{T2},D))
    eltype(E) == T2 || (E = convert(Matrix{T2},E))
    eltype(F) == T2 || (F = convert(Matrix{T2},F))
    realcase = T2 <: AbstractFloat
    transsylv = isa(A,Adjoint) && isa(B,Adjoint) && isa(D,Adjoint) && isa(E,Adjoint)
    realcase ? trans = 'T' : trans = 'C'
    if transsylv
       AS, DS, Q1, Z1 = schur(A.parent,D.parent)
       BS, ES, Q2, Z2 = schur(B.parent,E.parent)
       CS = adjoint(Z1) * (C*Z2)
       FS = adjoint(Q1) * (F*Q2)
       X, Y, scale =  tgsyl!(trans,AS,BS,CS,DS,ES,-FS)
       (rmul!(Q1*(X * adjoint(Z2)), inv(scale)), rmul!(Q1*(Y * adjoint(Z2)), inv(scale)) )
    else
       AS, DS, Q1, Z1 = schur(copy(A'),copy(D'))
       BS, ES, Q2, Z2 = schur(copy(B'),copy(E'))
       CS = adjoint(Z1) * (C*Z2)
       FS = adjoint(Q1) * (F*Q2)

       X, Y, scale =  tgsyl!(trans,AS,BS,CS,DS,ES,-FS)

       (rmul!(Q1*(X * adjoint(Z2)), inv(scale)), rmul!(Q1*(Y * adjoint(Z2)), inv(scale)) )
    end
end
"""
    sylvcs!(A,B,C; adjA = false, adjB = false)

Solve the continuous Sylvester matrix equation

                op(A)X + Xop(B) =  C,

where `op(A) = A` or `op(A) = A'` if `adjA = false` or `adjA = true`, respectively,
and `op(B) = B` or `op(B) = B'` if `adjB = false` or `adjB = true`, respectively.
`A` and `B` are square matrices in Schur forms, and `A` and `-B` must not have
common eigenvalues. `C` contains on output the solution `X`.
"""
function sylvcs!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::AbstractMatrix{T1}; adjA::Bool = false, adjB::Bool = false) where  T1<:BlasFloat
   """
   This is a wrapper to the LAPACK.trsylv! function, based on the Bartels-Stewart Schur form based approach.
   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   try
      trans = T1 <: Complex ? 'C' : 'T'
      C, scale = LAPACK.trsyl!(adjA ? trans : 'N', adjB ? trans : 'N', A, B, C)
      rmul!(C, inv(scale))
      return C[:,:]
   catch err
      findfirst("LAPACKException(1)",string(err)) === nothing ? rethrow() :
               throw("ME:SingularException: A has eigenvalue(s) α and B has eigenvalues(s) β such that α+β = 0")
   end
end
function sylvd2!(adjA::Bool, adjB::Bool, C::AbstractMatrix{T}, na::Int, nb::Int, A::AbstractMatrix{T}, B::AbstractMatrix{T}, Xw::AbstractMatrix{T}) where {T <:BlasReal}
   # speed and reduced allocation oriented implementation of a solver for 1x1 and 2x2 Sylvester equations
   # encountered in solving discrete Lyapunov equations:
   # A*X*B + X = C   if adjA = false and adjB = false -> R = kron(B',A) + I
   # A'*X*B + X = C  if adjA = true  and adjB = false -> R = kron(B',A') + I
   # A*X*B' + X = C  if adjA = false and adjB = true  -> R = kron(B,A) + I
   # A'*X*B' + X = C if adjA = true  and adjB = true  -> R = kron(B,A') + I
   ONE = one(T)
   if na == 1 && nb == 1
      temp = A[1,1]*B[1,1] + ONE
      rmul!(C,inv(temp))
      any(!isfinite, C) && throw("ME:SingularException: A has eigenvalue(s) α and B has eigenvalues(s) β such that αβ = -1")
      return C
   end
   n = na*nb
   i1 = 1:n
   R = view(Xw, i1, i1)
   Y = reshape(C, n)
   if adjA && !adjB
      if na == 1
         # R12 =
         # [ a11*b11+1      a11*b21]
         # [     a11*b12  a11*b22+1]
         # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,1]*B[2,1];
         #                 A[1,1]*B[1,2]  A[1,1]*B[2,2]+ONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
         @inbounds  R[1,2] = A[1,1]*B[2,1]
         @inbounds  R[2,1] = A[1,1]*B[1,2]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+ONE
         # @inbounds  Y[1] = C[1,1]
         # @inbounds  Y[2] = C[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11+1      a21*b11]
            # [     a12*b11  a22*b11+1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[2,1]*B[1,1];
            #                 A[1,2]*B[1,1]  A[2,2]*B[1,1]+ONE ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[2,1]*B[1,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
         else
            # R =
            # [ a11*b11+1      a21*b11      a11*b21      a21*b21]
            # [     a12*b11  a22*b11+1      a12*b21      a22*b21]
            # [     a11*b12      a21*b12  a11*b22+1      a21*b22]
            # [     a12*b12      a22*b12      a12*b22  a22*b22+1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[2,1]*B[1,1]      A[1,1]*B[2,1]      A[2,1]*B[2,1];
            # A[1,2]*B[1,1]  A[2,2]*B[1,1]+ONE      A[1,2]*B[2,1]      A[2,2]*B[2,1];
            # A[1,1]*B[1,2]      A[2,1]*B[1,2]  A[1,1]*B[2,2]+ONE      A[2,1]*B[2,2];
            # A[1,2]*B[1,2]      A[2,2]*B[1,2]      A[1,2]*B[2,2]  A[2,2]*B[2,2]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[2,1]*B[1,1]
            @inbounds  R[1,3] = A[1,1]*B[2,1]
            @inbounds  R[1,4] = A[2,1]*B[2,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            @inbounds  R[2,3] = A[1,2]*B[2,1]
            @inbounds  R[2,4] = A[2,2]*B[2,1]
            @inbounds  R[3,1] = A[1,1]*B[1,2]
            @inbounds  R[3,2] = A[2,1]*B[1,2]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+ONE
            @inbounds  R[3,4] = A[2,1]*B[2,2]
            @inbounds  R[4,1] = A[1,2]*B[1,2]
            @inbounds  R[4,2] = A[2,2]*B[1,2]
            @inbounds  R[4,3] = A[1,2]*B[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
            # @inbounds  Y[3] = C[1,2]
            # @inbounds  Y[4] = C[2,2]
         end
      end
      #R = kron(transpose(B),A) + I
   elseif !adjA && adjB
      if na == 1
         # R12 =
         # [ a11*b11+1      a11*b12]
         # [     a11*b21  a11*b22+1]
         # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,1]*B[1,2];
         #                 A[1,1]*B[2,1]  A[1,1]*B[2,2]+ONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
         @inbounds  R[1,2] = A[1,1]*B[1,2]
         @inbounds  R[2,1] = A[1,1]*B[2,1]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+ONE
         # @inbounds  Y[1] = C[1,1]
         # @inbounds  Y[2] = C[1,2]
      else
         if nb == 1
            # R21 =
            #    [ a11*b11+1      a12*b11]
            #    [     a21*b11  a22*b11+1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,2]*B[1,1];
            #                 A[2,1]*B[1,1]  A[2,2]*B[1,1]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[1,2]*B[1,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
         else
            # R =
            # [ a11*b11+1      a12*b11      a11*b12      a12*b12]
            # [     a21*b11  a22*b11+1      a21*b12      a22*b12]
            # [     a11*b21      a12*b21  a11*b22+1      a12*b22]
            # [     a21*b21      a22*b21      a21*b22  a22*b22+1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,2]*B[1,1]      A[1,1]*B[1,2]      A[1,2]*B[1,2];
            # A[2,1]*B[1,1]  A[2,2]*B[1,1]+ONE      A[2,1]*B[1,2]      A[2,2]*B[1,2];
            # A[1,1]*B[2,1]      A[1,2]*B[2,1]  A[1,1]*B[2,2]+ONE      A[1,2]*B[2,2];
            # A[2,1]*B[2,1]      A[2,2]*B[2,1]      A[2,1]*B[2,2]  A[2,2]*B[2,2]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[1,2]*B[1,1]
            @inbounds  R[1,3] = A[1,1]*B[1,2]
            @inbounds  R[1,4] = A[1,2]*B[1,2]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            @inbounds  R[2,3] = A[2,1]*B[1,2]
            @inbounds  R[2,4] = A[2,2]*B[1,2]
            @inbounds  R[3,1] = A[1,1]*B[2,1]
            @inbounds  R[3,2] = A[1,2]*B[2,1]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+ONE
            @inbounds  R[3,4] = A[1,2]*B[2,2]
            @inbounds  R[4,1] = A[2,1]*B[2,1]
            @inbounds  R[4,2] = A[2,2]*B[2,1]
            @inbounds  R[4,3] = A[2,1]*B[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
            # @inbounds  Y[3] = C[1,2]
            # @inbounds  Y[4] = C[2,2]
         end
      end
      #R = kron(transpose(B),transpose(A)) + I
   elseif !adjA && !adjB
      if na == 1
         # R12 =
         # [ a11*b11 + 1,     a11*b21]
         # [     a11*b12, a11*b22 + 1]
         # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,1]*B[2,1];
         #                 A[1,1]*B[1,2]  A[1,1]*B[2,2]+ONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
         @inbounds  R[1,2] = A[1,1]*B[2,1]
         @inbounds  R[2,1] = A[1,1]*B[1,2]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+ONE
         # @inbounds  Y[1] = C[1,1]
         # @inbounds  Y[2] = C[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + 1,     a12*b11]
            # [     a21*b11, a22*b11 + 1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,2]*B[1,1];
            #                 A[2,1]*B[1,1]  A[2,2]*B[1,1]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[1,2]*B[1,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
         else
            # R =
            # [ a11*b11 + 1,     a12*b11,     a11*b21,     a12*b21]
            # [     a21*b11, a22*b11 + 1,     a21*b21,     a22*b21]
            # [     a11*b12,     a12*b12, a11*b22 + 1,     a12*b22]
            # [     a21*b12,     a22*b12,     a21*b22, a22*b22 + 1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,2]*B[1,1]      A[1,1]*B[2,1]      A[1,2]*B[2,1];
            # A[2,1]*B[1,1]  A[2,2]*B[1,1]+ONE      A[2,1]*B[2,1]      A[2,2]*B[2,1];
            # A[1,1]*B[1,2]      A[1,2]*B[1,2]  A[1,1]*B[2,2]+ONE      A[1,2]*B[2,2];
            # A[2,1]*B[1,2]      A[2,2]*B[1,2]      A[2,1]*B[2,2]  A[2,2]*B[2,2]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[1,2]*B[1,1]
            @inbounds  R[1,3] = A[1,1]*B[2,1]
            @inbounds  R[1,4] = A[1,2]*B[2,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            @inbounds  R[2,3] = A[2,1]*B[2,1]
            @inbounds  R[2,4] = A[2,2]*B[2,1]
            @inbounds  R[3,1] = A[1,1]*B[1,2]
            @inbounds  R[3,2] = A[1,2]*B[1,2]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+ONE
            @inbounds  R[3,4] = A[1,2]*B[2,2]
            @inbounds  R[4,1] = A[2,1]*B[1,2]
            @inbounds  R[4,2] = A[2,2]*B[1,2]
            @inbounds  R[4,3] = A[2,1]*B[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
            # @inbounds  Y[3] = C[1,2]
            # @inbounds  Y[4] = C[2,2]
         end
      end
      #R = kron(B,A) + I
   else
      if na == 1
         # R12 =
         # [ a11*b11 + 1,     a11*b12]
         # [     a11*b21, a11*b22 + 1]
         # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[1,1]*B[1,2];
         #                 A[1,1]*B[2,1]  A[1,1]*B[2,2]+ONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
         @inbounds  R[1,2] = A[1,1]*B[1,2]
         @inbounds  R[2,1] = A[1,1]*B[2,1]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+ONE
         # @inbounds  Y[1] = C[1,1]
         # @inbounds  Y[2] = C[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + 1,     a21*b11]
            # [     a12*b11, a22*b11 + 1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[2,1]*B[1,1];
            #                 A[1,2]*B[1,1]  A[2,2]*B[1,1]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[2,1]*B[1,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
         else
            # R =
            # [ a11*b11 + 1,     a21*b11,     a11*b12,     a21*b12]
            # [     a12*b11, a22*b11 + 1,     a12*b12,     a22*b12]
            # [     a11*b21,     a21*b21, a11*b22 + 1,     a21*b22]
            # [     a12*b21,     a22*b21,     a12*b22, a22*b22 + 1]
            # @inbounds R = [ A[1,1]*B[1,1]+ONE      A[2,1]*B[1,1]      A[1,1]*B[1,2]      A[2,1]*B[1,2];
            # A[1,2]*B[1,1]  A[2,2]*B[1,1]+ONE      A[1,2]*B[1,2]      A[2,2]*B[1,2];
            # A[1,1]*B[2,1]      A[2,1]*B[2,1]  A[1,1]*B[2,2]+ONE      A[2,1]*B[2,2];
            # A[1,2]*B[2,1]      A[2,2]*B[2,1]      A[1,2]*B[2,2]  A[2,2]*B[2,2]+ONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+ONE
            @inbounds  R[1,2] = A[2,1]*B[1,1]
            @inbounds  R[1,3] = A[1,1]*B[1,2]
            @inbounds  R[1,4] = A[2,1]*B[1,2]
            @inbounds  R[2,1] = A[1,2]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+ONE
            @inbounds  R[2,3] = A[1,2]*B[1,2]
            @inbounds  R[2,4] = A[2,2]*B[1,2]
            @inbounds  R[3,1] = A[1,1]*B[2,1]
            @inbounds  R[3,2] = A[2,1]*B[2,1]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+ONE
            @inbounds  R[3,4] = A[2,1]*B[2,2]
            @inbounds  R[4,1] = A[1,2]*B[2,1]
            @inbounds  R[4,2] = A[2,2]*B[2,1]
            @inbounds  R[4,3] = A[1,2]*B[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+ONE
            # @inbounds  Y[1] = C[1,1]
            # @inbounds  Y[2] = C[2,1]
            # @inbounds  Y[3] = C[1,2]
            # @inbounds  Y[4] = C[2,2]
         end
         #R = kron(B,transpose(A)) + I
      end
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalue(s) α and B has eingenvalu(s) β such that αβ = -1")
   #C[:,:] = Y
   return C
end
"""
    sylvds!(A,B,C; adjA = false, adjB = false)

Solve the discrete Sylvester matrix equation

                op(A)Xop(B) + X =  C,

where `op(A) = A` or `op(A) = A'` if `adjA = false` or `adjA = true`, respectively,
and `op(B) = B` or `op(B) = B'` if `adjB = false` or `adjB = true`, respectively.
`A` and `B` are square matrices in Schur forms, and `A` and `-B` must not have
common reciprocal eigenvalues. `C` contains on output the solution `X`.
"""
function sylvds!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::AbstractMatrix{T1}, W::AbstractMatrix{T1} = similar(A,size(A,1),2); adjA::Bool = false, adjB::Bool = false) where  T1<:BlasReal
   """
   An extension of the Bartels-Stewart Schur form based approach is employed.

   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   m, n = LinearAlgebra.checksquare(A,B)
   (size(C,1) == m && size(C,2) == n ) || throw(DimensionMismatch("C must be an $m x $n matrix"))
   (m, 2) == size(W) || throw(DimensionMismatch("W must be an $m x 2 matrix"))
   ONE = one(T1)

   # determine the structure of the real Schur form of A
   ba, pa = sfstruct(A)
   bb, pb = sfstruct(B)

   G = Matrix{T1}(undef,2,2)
   WA = Matrix{T1}(undef,2,2)

   Xw = Matrix{T1}(undef,4,4)
   if !adjA && !adjB
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-left corner column by column by

      #            A(K,K)*X(K,L)*B(L,L) + X(K,L) = C(K,L) - R(K,L)

      # where
      #                        M
      #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L) +
      #                      J=K+1
      #                        M             L-1
      #                       SUM { A(K,J) * SUM [X(J,I)*B(I,L)] }.
      #                       J=K            I=1
      # """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = m
          for kk = pa:-1:1
              dk = ba[kk]
              dkk = 1:dk
              i1 = i-dk+1
              k = i1:i
              Ckl = view(C,k,l)
              y = view(G,1:dk,1:dl)
              copyto!(y,Ckl)
              if kk < pa
                 ir = i+1:m
                 W1 = view(WA,dkk,dll)
                 mul!(W1,view(A,k,ir),view(C,ir,l))
                 mul!(y,W1,view(B,l,l),-ONE,ONE)
              end
              if ll > 1
                 ic = i1:m
                 mul!(view(W,k,dll),view(C,k,il1),view(B,il1,l))
                 mul!(y,view(A,k,ic),view(W,ic,dll),-ONE,ONE)
              end
              sylvd2!(adjA,adjB,y,dk,dl,view(A,k,k),view(B,l,l),Xw)
              copyto!(Ckl,y)
           i -= dk
          end
          j += dl
      end
   elseif !adjA && adjB
         # """
         # The (K,L)th block of X is determined starting from
         # bottom-right corner column by column by

         #             A(K,K)*X(K,L)*B(L,L)' + X(K,L) = C(K,L) - R(K,L)

         # where
         #                        M
         #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L)' +
         #                      J=K+1
         #                        M              N
         #                       SUM { A(K,J) * SUM [X(J,I)*B(L,I)'] }.
         #                       J=K           I=L+1
         # """
         j = n
         for ll = pb:-1:1
             dl = bb[ll]
             dll = 1:dl
             il1 = j+1:n
             l = j-dl+1:j
             i = m
             for kk = pa:-1:1
                 dk = ba[kk]
                 dkk = 1:dk
                 i1 = i-dk+1
                 k = i1:i
                 Ckl = view(C,k,l)
                 y = view(G,1:dk,1:dl)
                 copyto!(y,Ckl)
                 if kk < pa
                    ir = i+1:m
                    W1 = view(WA,dkk,dll)
                    mul!(W1,view(A,k,ir),view(C,ir,l))
                    mul!(y,W1,adjoint(view(B,l,l)),-ONE,ONE)
                 end
                 if ll < pb
                    ic = i1:m
                    mul!(view(W,k,dll),view(C,k,il1),adjoint(view(B,l,il1)))
                    mul!(y,view(A,k,ic),view(W,ic,dll),-ONE,ONE)
                 end
                 sylvd2!(adjA,adjB,y,dk,dl,view(A,k,k),view(B,l,l),Xw)
                 copyto!(Ckl,y)
                 i -= dk
             end
             j -= dl
         end
   elseif adjA && !adjB
      # """
      # The (K,L)th block of X is determined starting from the
      # upper-left corner column by column by

      # A(K,K)'*X(K,L)*B(L,L) + X(K,L) = C(K,L) - R(K,L),

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(I,L)] }.
      #                       J=1             I=1
      # """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              Ckl = view(C,k,l)
              y = view(G,1:dk,1:dl)
              copyto!(y,Ckl)
              if kk > 1
                 ir = 1:i-1
                 W1 = view(WA,dkk,dll)
                 mul!(W1,adjoint(view(A,ir,k)),view(C,ir,l))
                 mul!(y,W1,view(B,l,l),-ONE,ONE)
              end
              if ll > 1
                 ic = 1:i1
                 mul!(view(W,k,dll),view(C,k,il1),view(B,il1,l))
                 mul!(y,adjoint(view(A,ic,k)),view(W,ic,dll),-ONE,ONE)
              end
              sylvd2!(adjA,adjB,y,dk,dl,view(A,k,k),view(B,l,l),Xw)
              copyto!(Ckl,y)
           i += dk
          end
          j += dl
      end
   elseif adjA && adjB
      # """
      # The (K,L)th block of X is determined starting from the
      # lower-left corner column by column by

      #            A(K,K)'*X(K,L)*B(L,L)' + X(K,L) = C(K,L) - R(K,L)

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(L,I)'] }.
      #                       J=1            I=L+1
      # """
      j = n
      for ll = pb:-1:1
          dl = bb[ll]
          dll = 1:dl
          il1 = j+1:n
          l = j-dl+1:j
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              Ckl = view(C,k,l)
              y = view(G,1:dk,1:dl)
              copyto!(y,Ckl)
              if kk > 1
                 ir = 1:i-1
                 W1 = view(WA,dkk,dll)
                 mul!(W1,adjoint(view(A,ir,k)),view(C,ir,l))
                 mul!(y,W1,adjoint(view(B,l,l)),-ONE,ONE)
              end
              if ll < pb
                 ic = 1:i1
                 mul!(view(W,k,dll),view(C,k,il1),adjoint(view(B,l,il1)))
                 mul!(y,adjoint(view(A,ic,k)),view(W,ic,dll),-ONE,ONE)
              end
              sylvd2!(adjA,adjB,y,dk,dl,view(A,k,k),view(B,l,l),Xw)
              copyto!(Ckl,y)
           i += dk
          end
          j -= dl
      end
   end
   return C
end
function sylvds!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::AbstractMatrix{T1}, W::AbstractVector{T1} = similar(A,size(A,1)); adjA::Bool = false, adjB::Bool = false) where  T1<:BlasComplex
   """
   An extension of the Bartels-Stewart Schur form based approach is employed.

   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   m, n = LinearAlgebra.checksquare(A,B)
   (size(C,1) == m && size(C,2) == n ) || throw(DimensionMismatch("C must be an $m x $n matrix"))

   ONE = one(T1)
   ZERO = zero(T1)
   if !adjA && !adjB
      # """
      # The (K,L)th element of X is determined starting from
      # bottom-left corner column by column by

      #            A(K,K)*X(K,L)*B(L,L) + X(K,L) = C(K,L) - R(K,L)

      # where
      #                        M
      #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L) +
      #                      J=K+1
      #                        M             L-1
      #                       SUM { A(K,J) * SUM [X(J,I)*B(I,L)] }.
      #                       J=K            I=1
      # """
      for l = 1:n
          il1 = 1:l-1
          for k = m:-1:1
              y = C[k,l]
              if k < m
                 ta = ZERO
                 for ir = k+1:m
                     ta += A[k,ir]*C[ir,l]
                 end
                 y -= ta*B[l,l]
              end
              if l > 1
                 tz = ZERO
                 for ir = il1
                     tz += C[k,ir]*B[ir,l]
                 end
                 W[k] = tz
                 for ic = k:m
                     y -= A[k,ic]*W[ic]
                 end
              end
              C[k,l] = y/(B[l,l]*A[k,k]+ONE)
              isfinite(C[k,l]) || throw("ME:SingularException: A and -B have common or close reciprocal eigenvalues")
             end
      end
   elseif !adjA && adjB
         # """
         # The (K,L)th element of X is determined starting from
         # bottom-right corner column by column by

         #          A(K,K)*X(K,L)*B(L,L)' + X(K,L) = C(K,L) - R(K,L)

         # where
         #                        M
         #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L)' +
         #                      J=K+1
         #                        M              N
         #                       SUM { A(K,J) * SUM [X(J,I)*B(L,I)'] }.
         #                       J=K           I=L+1
         # """
         for l = n:-1:1
             il1 = l+1:n
             for k = m:-1:1
                 y = C[k,l]
                 if k < m
                    ta = ZERO
                    for ir = k+1:m
                       ta += A[k,ir]*C[ir,l]
                    end
                    y -= ta*B[l,l]'
                 end
                 if l < n
                    tz = ZERO
                    for ir = il1
                        tz += C[k,ir]*B[l,ir]'
                    end
                    W[k] = tz
                    for ic = k:m
                        y -= A[k,ic]*W[ic]
                    end
                 end
                 C[k,l] = y/(B[l,l]'*A[k,k]+ONE)
                 isfinite(C[k,l]) || throw("ME:SingularException: A and -B have common or close reciprocal eigenvalues")
             end
         end
   elseif adjA && !adjB
      # """
      # The (K,L)th element of X is determined starting from the
      # upper-left corner column by column by

      #          A(K,K)'*X(K,L)*B(L,L) + X(K,L) = C(K,L) - R(K,L),

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(I,L)] }.
      #                       J=1             I=1
      # """
      for l = 1:n
          il1 = 1:l-1
          for k = 1:m
              y = C[k,l]
              if k > 1
                 ta = ZERO
                 for ir = 1:k-1
                     ta += A[ir,k]'*C[ir,l]
                 end
                 y -= ta*B[l,l]
              end
              if l > 1
                 ic = 1:m
                 tz = ZERO
                 for ir = il1
                     tz += C[k,ir]*B[ir,l]
                 end
                 W[k] = tz
                 for ic = 1:m
                     y -= A[ic,k]'*W[ic]
                 end
              end
              C[k,l] = y/(B[l,l]*A[k,k]'+ONE)
              isfinite(C[k,l]) || throw("ME:SingularException: A and -B have common or close reciprocal eigenvalues")
         end
      end
   elseif adjA && adjB
      # """
      # The (K,L)th element of X is determined starting from the
      # upper-right corner column by column by

      #         A(K,K)'*X(K,L)*B(L,L)' + X(K,L) = C(K,L) - R(K,L)

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(L,I)'] }.
      #                       J=1            I=L+1
      # """
      for l = n:-1:1
          il1 = l+1:n
          for k = 1:m
              y = C[k,l]
              if k > 1
                 ta = ZERO
                 for ir = 1:k-1
                     ta += A[ir,k]'*C[ir,l]
                 end
                 y -= ta*B[l,l]'
              end
              if l < n
                 ic = 1:m
                 tz = ZERO
                 for ir = il1
                     tz += C[k,ir]*B[l,ir]'
                 end
                 W[k] = tz
                 for ic = 1:m
                     y -= A[ic,k]'*W[ic]
                 end
              end
              C[k,l] = y/(B[l,l]'*A[k,k]'+ONE)
              isfinite(C[k,l]) || throw("ME:SingularException: A and -B have common or close reciprocal eigenvalues")
          end
      end
   end
   return C
end
"""
    X = gsylvs!(A,B,C,D,E; adjAC=false, adjBD=false, CASchur = false, DBSchur = false)

Solve the generalized Sylvester matrix equation

                op1(A)Xop2(B) + op1(C)Xop2(D) = E,

where `A`, `B`, `C` and `D` are square matrices, and

`op1(A) = A` and `op1(C) = C` if `adjAC = false`;

`op1(A) = A'` and `op1(C) = C'` if `adjAC = true`;

`op2(B) = B` and `op2(D) = D` if `adjBD = false`;

`op2(B) = B'` and `op2(D) = D'` if `adjBD = true`.

The matrix pair `(A,C)` is in a generalized real or complex Schur form.
The matrix pair `(B,D)` is in a generalized real or complex Schur form if `DBSchur = false`
or the matrix pair `(D,B)` is in a generalized real or complex Schur form if `DBSchur = true`.
The pencils `A-λC` and `D+λB` must be regular and must not have common eigenvalues.
"""
function gsylvs!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::AbstractMatrix{T1}, D::AbstractMatrix{T1}, E::AbstractMatrix{T1},
                 WB::AbstractMatrix{T1} = similar(A,size(A,1),2), WD::AbstractMatrix{T1} = similar(A,size(A,1),2);
                 adjAC::Bool = false, adjBD::Bool = false, CASchur::Bool = false, DBSchur::Bool = false) where T1<:BlasReal
   """
   An extension proposed in [1] of the Bartels-Stewart Schur form based approach [2] is employed.

   References:
   [1] K.-W. E. Chu. The solution of the matrix equation AXB – CXD = E and
       (YA – DZ, YC– BZ) = (E, F). Lin. Alg. Appl., 93:93-105, 1987.
   [2] R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
       Comm. ACM, 15:820–826, 1972.
   """
   m, n = size(E);
   [m; n; m; n] == LinearAlgebra.checksquare(A,B,C,D) ||
      throw(DimensionMismatch("A, B, C, D and E have incompatible dimensions"))

   (m, 2) == size(WB) || throw(DimensionMismatch("WB must be an $m x 2 matrix"))
   (m, 2) == size(WD) || throw(DimensionMismatch("WD must be an $m x 2 matrix"))

   ONE = one(T1)

   # determine the structure of the generalized real Schur form of (A,C)
   CASchur ? ((ba, pa) = sfstruct(C)) : ((ba, pa) = sfstruct(A))
   DBSchur ? ((bb, pb) = sfstruct(D)) : ((bb, pb) = sfstruct(B))

   # WB = zeros(T1,m,2)
   # WD = zeros(T1,m,2)
   Xw = Matrix{T1}(undef,4,4)
   if !adjAC && !adjBD
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-left corner column by column by

      #       A(K,K)*X(K,L)*B(L,L) + C(K,K)*X(K,L)*D(L,L) = E(K,L) - R(K,L)

      # where
      #                        M
      #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L) +
      #                      J=K+1
      #                        M             L-1
      #                       SUM { A(K,J) * SUM [X(J,I)*B(I,L)] } +
      #                       J=K            I=1

      #                        M
      #                     { SUM [C(K,J)*X(J,L)] } * D(L,L) +
      #                      J=K+1
      #                        M             L-1
      #                       SUM { C(K,J) * SUM [X(J,I)*D(I,L)] }.
      #                       J=K            I=1
      # """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = m
          for kk = pa:-1:1
              dk = ba[kk]
              dkk = 1:dk
              i1 = i-dk+1
              k = i1:i
              y = view(E,k,l)
              W1 = view(Xw,dkk,dll)
              if kk < pa
                 ir = i+1:m
                 #W1 = A[k,ir]*E[ir,l]
                 mul!(W1,view(A,k,ir),view(E,ir,l))
                 #y -= W1*B[l,l]
                 mul!(y,W1,view(B,l,l),-ONE,ONE)
                 #W1 = C[k,ir]*E[ir,l]
                 mul!(W1,view(C,k,ir),view(E,ir,l))
                 #y -= W1*D[l,l]
                 mul!(y,W1,view(D,l,l),-ONE,ONE)
              end
              if ll > 1
                 ic = i1:m
                 # WB[k,dll] = E[k,il1]*B[il1,l]
                 mul!(view(WB,k,dll),view(E,k,il1),view(B,il1,l))
                 # WD[k,dll] = E[k,il1]*D[il1,l]
                 mul!(view(WD,k,dll),view(E,k,il1),view(D,il1,l))
                 # y -= (A[k,ic]*WB[ic,dll] + C[k,ic]*WD[ic,dll])
                 mul!(y,view(A,k,ic),view(WB,ic,dll),-ONE,ONE)
                 mul!(y,view(C,k,ic),view(WD,ic,dll),-ONE,ONE)
              end
              gsylv2!(adjAC,adjBD,y,dk,dl,view(A,k,k),view(B,l,l),view(C,k,k),view(D,l,l),Xw)
              i -= dk
          end
          j += dl
      end
   elseif !adjAC && adjBD
         # """
         #  The (K,L)th block of X is determined starting from
         #  bottom-right corner column by column by

         #       A(K,K)*X(K,L)*B(L,L)' + C(K,K)*X(K,L)*D(L,L)' = E(K,L) - R(K,L)

         #  where
         #                        M
         #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L)' +
         #                      J=K+1
         #                        M              N
         #                       SUM { A(K,J) * SUM [X(J,I)*B(L,I)'] } +
         #                       J=K           I=L+1

         #                       M
         #                    { SUM [C(K,J)*X(J,L)] } * D(L,L)' +
         #                     J=K+1
         #                       M              N
         #                      SUM { C(K,J) * SUM [X(J,I)*D(L,I)'] }.
         #                      J=K           I=L+1
         # """
         j = n
         for ll = pb:-1:1
             dl = bb[ll]
             dll = 1:dl
             il1 = j+1:n
             l = j-dl+1:j
             i = m
             for kk = pa:-1:1
                 dk = ba[kk]
                 dkk = 1:dk
                 i1 = i-dk+1
                 k = i1:i
                 y = view(E,k,l)
                 W1 = view(Xw,dkk,dll)
                 if kk < pa
                    ir = i+1:m
                    #W1 = A[k,ir]*E[ir,l]
                    mul!(W1,view(A,k,ir),view(E,ir,l))
                    #y -= W1*B[l,l]'
                    mul!(y,W1,transpose(view(B,l,l)),-ONE,ONE)
                    #W2 = C[k,ir]*E[ir,l]
                    mul!(W1,view(C,k,ir),view(E,ir,l))
                    # y -= W1*D[l,l]'
                    mul!(y,W1,transpose(view(D,l,l)),-ONE,ONE)
                  end
                 if ll < pb
                    ic = i1:m
                    #WB[k,dll] = E[k,il1]*B[l,il1]'
                    mul!(view(WB,k,dll),view(E,k,il1),transpose(view(B,l,il1)))
                    #WD[k,dll] = E[k,il1]*D[l,il1]'
                    mul!(view(WD,k,dll),view(E,k,il1),transpose(view(D,l,il1)))
                    #y -= (A[k,ic]*WB[ic,dll]+C[k,ic]*WD[ic,dll])
                    mul!(y,view(A,k,ic),view(WB,ic,dll),-ONE,ONE)
                    mul!(y,view(C,k,ic),view(WD,ic,dll),-ONE,ONE)
                 end
                 gsylv2!(adjAC,adjBD,y,dk,dl,view(A,k,k),view(B,l,l),view(C,k,k),view(D,l,l),Xw)
                 i -= dk
             end
             j -= dl
         end
   elseif adjAC && !adjBD
      # """
      # The (K,L)th block of X is determined starting from the
      # upper-left corner column by column by

      # A(K,K)'*X(K,L)*B(L,L) + C(K,K)'*X(K,L)*D(L,L) = E(K,L) - R(K,L),

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(I,L)] } +
      #                       J=1             I=1

      #                       K-1
      #                     { SUM [C(J,K)'*X(J,L)] } * D(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM C(J,K)' * { SUM [X(J,I)*D(I,L)] }.
      #                       J=1             I=1
      # """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = view(E,k,l)
              W1 = view(Xw,dkk,dll)
              if kk > 1
                 ir = 1:i-1
                 # W1 = A[ir,k]'*E[ir,l]
                 mul!(W1,transpose(view(A,ir,k)),view(E,ir,l))
                 #y -= W1*B[l,l]
                 mul!(y,W1,view(B,l,l),-ONE,ONE)
                 mul!(W1,transpose(view(C,ir,k)),view(E,ir,l))
                 #y -= W1*D[l,l]
                 mul!(y,W1,view(D,l,l),-ONE,ONE)
              end
              if ll > 1
                 ic = 1:i1
                 #WB[k,dll] = E[k,il1]*B[il1,l]
                 mul!(view(WB,k,dll),view(E,k,il1),view(B,il1,l))
                 #y -= A[ic,k]'*WB[ic,dll]
                 mul!(y,transpose(view(A,ic,k)),view(WB,ic,dll),-ONE,ONE)
                 #WD[k,dll] = E[k,il1]*D[il1,l]
                 mul!(view(WD,k,dll),view(E,k,il1),view(D,il1,l))
                 #y -= C[ic,k]'*WD[ic,dll]
                 mul!(y,transpose(view(C,ic,k)),view(WD,ic,dll),-ONE,ONE)
              end
              gsylv2!(adjAC,adjBD,y,dk,dl,view(A,k,k),view(B,l,l),view(C,k,k),view(D,l,l),Xw)
              i += dk
          end
          j += dl
      end
   elseif adjAC && adjBD
      # """
      # The (K,L)th block of X is determined starting from
      # upper-right corner column by column by

      #            A(K,K)'*X(K,L)*B(L,L)' + C(K,K)'*X(K,L)*D(L,L)' = E(K,L) - R(K,L)

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(L,I)'] }+
      #                       J=1            I=L+1

      #                       K-1
      #                     { SUM [C(J,K)'*X(J,L)] } * D(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM C(J,K)' * { SUM [X(J,I)*D(L,I)'] }.
      #                       J=1            I=L+1
      # """
      j = n
      for ll = pb:-1:1
          dl = bb[ll]
          dll = 1:dl
          il1 = j+1:n
          l = j-dl+1:j
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = view(E,k,l)
              W1 = view(Xw,dkk,dll)
              if kk > 1
                 ir = 1:i-1
                 #W1 = A[ir,k]'*E[ir,l]
                 mul!(W1,transpose(view(A,ir,k)),view(E,ir,l))
                 #y -= W1*B[l,l]'
                 mul!(y,W1,transpose(view(B,l,l)),-ONE,ONE)
                 #   W2 = C[ir,k]'*E[ir,l]
               #   y -= W2*D[l,l]'
                 mul!(W1,transpose(view(C,ir,k)),view(E,ir,l))
                 #y -= W1*D[l,l]'
                 mul!(y,W1,transpose(view(D,l,l)),-ONE,ONE)
               end
              if ll < pb
                 ic = 1:i1
                 #WB[k,dll] = E[k,il1]*B[l,il1]'
                 mul!(view(WB,k,dll),view(E,k,il1),transpose(view(B,l,il1)))
                 #WD[k,dll] = E[k,il1]*D[l,il1]'
                 mul!(view(WD,k,dll),view(E,k,il1),transpose(view(D,l,il1)))
                 #y -= (A[ic,k]'*WB[ic,dll] + C[ic,k]'*WD[ic,dll])
                 mul!(y,transpose(view(A,ic,k)),view(WB,ic,dll),-ONE,ONE)
                 mul!(y,transpose(view(C,ic,k)),view(WD,ic,dll),-ONE,ONE)
              end
              gsylv2!(adjAC,adjBD,y,dk,dl,view(A,k,k),view(B,l,l),view(C,k,k),view(D,l,l),Xw)
              i += dk
          end
          j -= dl
      end
   end
   return E
end
@inline function gsylv2!(adjAC::Bool,adjBD::Bool,E::StridedMatrix{T},na::Int,nb::Int,A::AbstractMatrix{T},B::AbstractMatrix{T},C::AbstractMatrix{T},D::AbstractMatrix{T},Xw::StridedMatrix{T}) where T <:BlasReal
   # speed and reduced allocation oriented implementation of a solver for 1x1 and 2x2 generalized Sylvester equations:
   #      A*X*B + C*X*D = E     if adjAC = false and adjBD = false -> R = kron(B',A)  + kron(D',C)
   #      A'*X*B + C'*X*D = E   if adjAC = true and adjBD = false  -> R = kron(B',A') + kron(D',C')
   #      A*X*B' + C*X*D' = E   if adjAC = false and adjBD = true  -> R = kron(B,A)   + kron(D,C)
   #      A'*X*B' + C'*X*D' = E if adjAC = true and adjBD = true   -> R = kron(B,A')  + kron(D,C')
   if na == 1 && nb == 1
      temp = A[1,1]*B[1,1] + C[1,1]*D[1,1]
      rmul!(E,inv(temp))
      any(!isfinite, E) &&  throw("ME:SingularException: `A-λC` and `D+λB` have common eigenvalues")
      return E
   end
   n = na*nb
   i1 = 1:n
   R = view(Xw,i1,i1)
   Y = reshape(E, n)
   if !adjAC && !adjBD
      if na == 1
         # R12 =
         # [ a11*b11 + c11*d11, a11*b21 + c11*d21]
         # [ a11*b12 + c11*d12, a11*b22 + c11*d22]
         # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,1]*B[2,1]+C[1,1]*D[2,1];
         #                 A[1,1]*B[1,2]+C[1,1]*D[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
         @inbounds  R[1,2] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
         @inbounds  R[2,1] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
         # @inbounds  Y[1] = E[1,1]
         # @inbounds  Y[2] = E[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + c11*d11, a12*b11 + c12*d11]
            # [ a21*b11 + c21*d11, a22*b11 + c22*d11]
            # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]+C[1,2]*D[1,1];
            #                 A[2,1]*B[1,1]+C[2,1]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1] ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
         else
            # R =
            # [ a11*b11 + c11*d11, a12*b11 + c12*d11, a11*b21 + c11*d21, a12*b21 + c12*d21]
            # [ a21*b11 + c21*d11, a22*b11 + c22*d11, a21*b21 + c21*d21, a22*b21 + c22*d21]
            # [ a11*b12 + c11*d12, a12*b12 + c12*d12, a11*b22 + c11*d22, a12*b22 + c12*d22]
            # [ a21*b12 + c21*d12, a22*b12 + c22*d12, a21*b22 + c21*d22, a22*b22 + c22*d22]
            # (iszero(C[2,1]) && iszero(D[2,1]) && iszero(C[1,2]) && iszero(D[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]      A[1,1]*B[2,1]      A[1,2]*B[2,1];
            # A[2,1]*B[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[2,1]*B[2,1]      A[2,2]*B[2,1];
            # A[1,1]*B[1,2]      A[1,2]*B[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[1,2]*B[2,2];
            # A[2,1]*B[1,2]      A[2,2]*B[1,2]      A[2,1]*B[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]+C[1,2]*D[1,1]      A[1,1]*B[2,1]+C[1,1]*D[2,1]      A[1,2]*B[2,1]+C[1,2]*D[2,1];
            # A[2,1]*B[1,1]+C[2,1]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[2,1]*B[2,1]+C[2,1]*D[2,1]      A[2,2]*B[2,1]+C[2,2]*D[2,1];
            # A[1,1]*B[1,2]+C[1,1]*D[1,2]      A[1,2]*B[1,2]+C[1,2]*D[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[1,2]*B[2,2]+C[1,2]*D[2,2];
            # A[2,1]*B[1,2]+C[2,1]*D[1,2]      A[2,2]*B[1,2]+C[2,2]*D[1,2]      A[2,1]*B[2,2]+C[2,1]*D[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[1,3] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
            @inbounds  R[1,4] = A[1,2]*B[2,1]+C[1,2]*D[2,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            @inbounds  R[2,3] = A[2,1]*B[2,1]+C[2,1]*D[2,1]
            @inbounds  R[2,4] = A[2,2]*B[2,1]+C[2,2]*D[2,1]
            @inbounds  R[3,1] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
            @inbounds  R[3,2] = A[1,2]*B[1,2]+C[1,2]*D[1,2]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
            @inbounds  R[3,4] = A[1,2]*B[2,2]+C[1,2]*D[2,2]
            @inbounds  R[4,1] = A[2,1]*B[1,2]+C[2,1]*D[1,2]
            @inbounds  R[4,2] = A[2,2]*B[1,2]+C[2,2]*D[1,2]
            @inbounds  R[4,3] = A[2,1]*B[2,2]+C[2,1]*D[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+C[2,2]*D[2,2]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
            # @inbounds  Y[3] = E[1,2]
            # @inbounds  Y[4] = E[2,2]
         end
      end
      #R = kron(transpose(B),A) + kron(transpose(D),C)
   elseif adjAC && !adjBD
      if na == 1
         # R12 =
         # [ a11*b11 + c11*d11, a11*b21 + c11*d21]
         # [ a11*b12 + c11*d12, a11*b22 + c11*d22]
         # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,1]*B[2,1]+C[1,1]*D[2,1];
         #                 A[1,1]*B[1,2]+C[1,1]*D[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
         @inbounds  R[1,2] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
         @inbounds  R[2,1] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
         # @inbounds  Y[1] = E[1,1]
         # @inbounds  Y[2] = E[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + c11*d11, a21*b11 + c21*d11]
            # [ a12*b11 + c12*d11, a22*b11 + c22*d11]
            # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]+C[2,1]*D[1,1];
            #                 A[1,2]*B[1,1]+C[1,2]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1] ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
         else
            # R =
            # [ a11*b11 + c11*d11, a21*b11 + c21*d11, a11*b21 + c11*d21, a21*b21 + c21*d21]
            # [ a12*b11 + c12*d11, a22*b11 + c22*d11, a12*b21 + c12*d21, a22*b21 + c22*d21]
            # [ a11*b12 + c11*d12, a21*b12 + c21*d12, a11*b22 + c11*d22, a21*b22 + c21*d22]
            # [ a12*b12 + c12*d12, a22*b12 + c22*d12, a12*b22 + c12*d22, a22*b22 + c22*d22]
            # (iszero(C[2,1]) && iszero(D[2,1]) && iszero(C[1,2]) && iszero(D[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]      A[1,1]*B[2,1]      A[2,1]*B[2,1];
            # A[1,2]*B[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[1,2]*B[2,1]      A[2,2]*B[2,1];
            # A[1,1]*B[1,2]      A[2,1]*B[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[2,1]*B[2,2];
            # A[1,2]*B[1,2]      A[2,2]*B[1,2]      A[1,2]*B[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]+C[2,1]*D[1,1]      A[1,1]*B[2,1]+C[1,1]*D[2,1]      A[2,1]*B[2,1]+C[2,1]*D[2,1];
            # A[1,2]*B[1,1]+C[1,2]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[1,2]*B[2,1]+C[1,2]*D[2,1]      A[2,2]*B[2,1]+C[2,2]*D[2,1];
            # A[1,1]*B[1,2]+C[1,1]*D[1,2]      A[2,1]*B[1,2]+C[2,1]*D[1,2]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[2,1]*B[2,2]+C[2,1]*D[2,2];
            # A[1,2]*B[1,2]+C[1,2]*D[1,2]      A[2,2]*B[1,2]+C[2,2]*D[1,2]      A[1,2]*B[2,2]+C[1,2]*D[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[1,3] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
            @inbounds  R[1,4] = A[2,1]*B[2,1]+C[2,1]*D[2,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            @inbounds  R[2,3] = A[1,2]*B[2,1]+C[1,2]*D[2,1]
            @inbounds  R[2,4] = A[2,2]*B[2,1]+C[2,2]*D[2,1]
            @inbounds  R[3,1] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
            @inbounds  R[3,2] = A[2,1]*B[1,2]+C[2,1]*D[1,2]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
            @inbounds  R[3,4] = A[2,1]*B[2,2]+C[2,1]*D[2,2]
            @inbounds  R[4,1] = A[1,2]*B[1,2]+C[1,2]*D[1,2]
            @inbounds  R[4,2] = A[2,2]*B[1,2]+C[2,2]*D[1,2]
            @inbounds  R[4,3] = A[1,2]*B[2,2]+C[1,2]*D[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+C[2,2]*D[2,2]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
            # @inbounds  Y[3] = E[1,2]
            # @inbounds  Y[4] = E[2,2]
         end
      end
      #R = kron(transpose(B),transpose(A)) + kron(transpose(D),transpose(C))
   elseif !adjAC && adjBD
      if na == 1
         # R12 =
         # [ a11*b11 + c11*d11, a11*b12 + c11*d12]
         # [ a11*b21 + c11*d21, a11*b22 + c11*d22]
         # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,1]*B[1,2]+C[1,1]*D[1,2];
         #                 A[1,1]*B[2,1]+C[1,1]*D[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
         @inbounds  R[1,2] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
         @inbounds  R[2,1] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
         # @inbounds  Y[1] = E[1,1]
         # @inbounds  Y[2] = E[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + c11*d11, a12*b11 + c12*d11]
            # [ a21*b11 + c21*d11, a22*b11 + c22*d11]
            # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]+C[1,2]*D[1,1];
            #                 A[2,1]*B[1,1]+C[2,1]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1] ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
         else
            # R =
            # [ a11*b11 + c11*d11, a12*b11 + c12*d11, a11*b12 + c11*d12, a12*b12 + c12*d12]
            # [ a21*b11 + c21*d11, a22*b11 + c22*d11, a21*b12 + c21*d12, a22*b12 + c22*d12]
            # [ a11*b21 + c11*d21, a12*b21 + c12*d21, a11*b22 + c11*d22, a12*b22 + c12*d22]
            # [ a21*b21 + c21*d21, a22*b21 + c22*d21, a21*b22 + c21*d22, a22*b22 + c22*d22]
            # (iszero(C[2,1]) && iszero(D[2,1]) && iszero(C[1,2]) && iszero(D[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]      A[1,1]*B[1,2]      A[1,2]*B[1,2];
            # A[2,1]*B[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[2,1]*B[1,2]      A[2,2]*B[1,2];
            # A[1,1]*B[2,1]      A[1,2]*B[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[1,2]*B[2,2]+C[1,2]*D[2,2];
            # A[2,1]*B[2,1]      A[2,2]*B[2,1]      A[2,1]*B[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,2]*B[1,1]+C[1,2]*D[1,1]      A[1,1]*B[1,2]+C[1,1]*D[1,2]      A[1,2]*B[1,2]+C[1,2]*D[1,2];
            # A[2,1]*B[1,1]+C[2,1]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[2,1]*B[1,2]+C[2,1]*D[1,2]      A[2,2]*B[1,2]+C[2,2]*D[1,2];
            # A[1,1]*B[2,1]+C[1,1]*D[2,1]      A[1,2]*B[2,1]+C[1,2]*D[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[1,2]*B[2,2]+C[1,2]*D[2,2];
            # A[2,1]*B[2,1]+C[2,1]*D[2,1]      A[2,2]*B[2,1]+C[2,2]*D[2,1]      A[2,1]*B[2,2]+C[2,1]*D[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[1,3] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
            @inbounds  R[1,4] = A[1,2]*B[1,2]+C[1,2]*D[1,2]
            @inbounds  R[2,1] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            @inbounds  R[2,3] = A[2,1]*B[1,2]+C[2,1]*D[1,2]
            @inbounds  R[2,4] = A[2,2]*B[1,2]+C[2,2]*D[1,2]
            @inbounds  R[3,1] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
            @inbounds  R[3,2] = A[1,2]*B[2,1]+C[1,2]*D[2,1]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
            @inbounds  R[3,4] = A[1,2]*B[2,2]+C[1,2]*D[2,2]
            @inbounds  R[4,1] = A[2,1]*B[2,1]+C[2,1]*D[2,1]
            @inbounds  R[4,2] = A[2,2]*B[2,1]+C[2,2]*D[2,1]
            @inbounds  R[4,3] = A[2,1]*B[2,2]+C[2,1]*D[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+C[2,2]*D[2,2]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
            # @inbounds  Y[3] = E[1,2]
            # @inbounds  Y[4] = E[2,2]
         end
      end
      #R = kron(B,A) + kron(D,C)
   else
      if na == 1
         # R12 =
         # [ a11*b11 + c11*d11, a11*b12 + c11*d12]
         # [ a11*b21 + c11*d21, a11*b22 + c11*d22]
         # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[1,1]*B[1,2]+C[1,1]*D[1,2];
         #                 A[1,1]*B[2,1]+C[1,1]*D[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
         @inbounds  R[1,2] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
         @inbounds  R[2,1] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
         @inbounds  R[2,2] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
         # @inbounds  Y[1] = E[1,1]
         # @inbounds  Y[2] = E[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 + c11*d11, a21*b11 + c21*d11]
            # [ a12*b11 + c12*d11, a22*b11 + c22*d11]
            # @inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]+C[2,1]*D[1,1];
            #                 A[1,2]*B[1,1]+C[1,2]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1] ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
         else
            # R =
            # [ a11*b11 + c11*d11, a21*b11 + c21*d11, a11*b12 + c11*d12, a21*b12 + c21*d12]
            # [ a12*b11 + c12*d11, a22*b11 + c22*d11, a12*b12 + c12*d12, a22*b12 + c22*d12]
            # [ a11*b21 + c11*d21, a21*b21 + c21*d21, a11*b22 + c11*d22, a21*b22 + c21*d22]
            # [ a12*b21 + c12*d21, a22*b21 + c22*d21, a12*b22 + c12*d22, a22*b22 + c22*d22]
            # (iszero(C[2,1]) && iszero(D[2,1]) && iszero(C[1,2]) && iszero(D[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]      A[1,1]*B[1,2]      A[2,1]*B[1,2];
            # A[1,2]*B[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[1,2]*B[1,2]      A[2,2]*B[1,2];
            # A[1,1]*B[2,1]      A[2,1]*B[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[2,1]*B[2,2];
            # A[1,2]*B[2,1]      A[2,2]*B[2,1]      A[1,2]*B[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]+C[1,1]*D[1,1]      A[2,1]*B[1,1]+C[2,1]*D[1,1]      A[1,1]*B[1,2]+C[1,1]*D[1,2]      A[2,1]*B[1,2]+C[2,1]*D[1,2];
            # A[1,2]*B[1,1]+C[1,2]*D[1,1]  A[2,2]*B[1,1]+C[2,2]*D[1,1]      A[1,2]*B[1,2]+C[1,2]*D[1,2]      A[2,2]*B[1,2]+C[2,2]*D[1,2];
            # A[1,1]*B[2,1]+C[1,1]*D[2,1]      A[2,1]*B[2,1]+C[2,1]*D[2,1]  A[1,1]*B[2,2]+C[1,1]*D[2,2]      A[2,1]*B[2,2]+C[2,1]*D[2,2];
            # A[1,2]*B[2,1]+C[1,2]*D[2,1]      A[2,2]*B[2,1]+C[2,2]*D[2,1]      A[1,2]*B[2,2]+C[1,2]*D[2,2]  A[2,2]*B[2,2]+C[2,2]*D[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]+C[1,1]*D[1,1]
            @inbounds  R[1,2] = A[2,1]*B[1,1]+C[2,1]*D[1,1]
            @inbounds  R[1,3] = A[1,1]*B[1,2]+C[1,1]*D[1,2]
            @inbounds  R[1,4] = A[2,1]*B[1,2]+C[2,1]*D[1,2]
            @inbounds  R[2,1] = A[1,2]*B[1,1]+C[1,2]*D[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]+C[2,2]*D[1,1]
            @inbounds  R[2,3] = A[1,2]*B[1,2]+C[1,2]*D[1,2]
            @inbounds  R[2,4] = A[2,2]*B[1,2]+C[2,2]*D[1,2]
            @inbounds  R[3,1] = A[1,1]*B[2,1]+C[1,1]*D[2,1]
            @inbounds  R[3,2] = A[2,1]*B[2,1]+C[2,1]*D[2,1]
            @inbounds  R[3,3] = A[1,1]*B[2,2]+C[1,1]*D[2,2]
            @inbounds  R[3,4] = A[2,1]*B[2,2]+C[2,1]*D[2,2]
            @inbounds  R[4,1] = A[1,2]*B[2,1]+C[1,2]*D[2,1]
            @inbounds  R[4,2] = A[2,2]*B[2,1]+C[2,2]*D[2,1]
            @inbounds  R[4,3] = A[1,2]*B[2,2]+C[1,2]*D[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]+C[2,2]*D[2,2]
            # @inbounds  Y[1] = E[1,1]
            # @inbounds  Y[2] = E[2,1]
            # @inbounds  Y[3] = E[1,2]
            # @inbounds  Y[4] = E[2,2]
         end
      end
      #R = kron(B,transpose(A)) + kron(D,transpose(C))
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalue(s) α and B has eingenvalu(s) β such that αβ = -1")
   # E[:,:] = Y
   return E
end

function gsylvs!(A::AbstractMatrix{T1}, B::AbstractMatrix{T1}, C::AbstractMatrix{T1}, D::AbstractMatrix{T1}, E::AbstractMatrix{T1},
                 WB::AbstractVector{T1} = similar(A,size(A,1)), WD::AbstractVector{T1} = similar(A,size(A,1));
                 adjAC::Bool = false, adjBD::Bool = false, CASchur::Bool = false, DBSchur::Bool = false) where T1<:BlasComplex
   """
   An extension proposed in [1] of the Bartels-Stewart Schur form based approach [2] is employed.

   References:
   [1] K.-W. E. Chu. The solution of the matrix equation AXB – CXD = E and
       (YA – DZ, YC– BZ) = (E, F). Lin. Alg. Appl., 93:93-105, 1987.
   [2] R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
       Comm. ACM, 15:820–826, 1972.
   """
   m, n = size(E);
   [m; n; m; n] == LinearAlgebra.checksquare(A,B,C,D) ||
      throw(DimensionMismatch("A, B, C, D and E have incompatible dimensions"))
   m == length(WB) || throw(DimensionMismatch("WB must be an $m - dimensional vector"))
   m == length(WD) || throw(DimensionMismatch("WD must be an $m - dimensional vector"))

   WB = zeros(T1,m)
   WD = zeros(T1,m)
   ZERO = zero(T1)
   if !adjAC && !adjBD
      # """
      # The (K,L)th element of X is determined starting from
      # bottom-left corner column by column by

      #       A(K,K)*X(K,L)*B(L,L) +C(K,K)*X(K,L)*D(L,L) = E(K,L) - R(K,L)

      # where
      #                        M
      #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L) +
      #                      J=K+1
      #                        M             L-1
      #                       SUM { A(K,J) * SUM [X(J,I)*B(I,L)] } +
      #                       J=K            I=1

      #                       M
      #                    { SUM [C(K,J)*X(J,L)] } * D(L,L) +
      #                     J=K+1
      #                       M             L-1
      #                      SUM { C(K,J) * SUM [X(J,I)*D(I,L)] } +
      #                      J=K            I=1
      # """
      for l = 1:n
          il1 = 1:l-1
          for k = m:-1:1
              y = E[k,l]
              if k < m
                 ta = ZERO
                 tc = ZERO
                 for ir = k+1:m
                     ta += A[k,ir]*E[ir,l]
                     tc += C[k,ir]*E[ir,l]
                 end
                 y -= (ta*B[l,l]+tc*D[l,l])
              end
              if l > 1
                 ta = ZERO
                 tc = ZERO
                 for ir = il1
                     ta += E[k,ir]*B[ir,l]
                     tc += E[k,ir]*D[ir,l]
                 end
                 WB[k] = ta
                 WD[k] = tc
                 for ic = k:m
                     y -= (A[k,ic]*WB[ic]+C[k,ic]*WD[ic])
                 end
              end
              E[k,l] = y/(B[l,l]*A[k,k]+D[l,l]*C[k,k])
              isfinite(E[k,l]) || throw("ME:SingularException: A-λC and D+λB have common or close eigenvalues")
          end
      end
   elseif !adjAC && adjBD
         # """
         #  The (K,L)th element of X is determined starting from
         #  bottom-right corner column by column by

         #       A(K,K)*X(K,L)*B(L,L)' + C(K,K)*X(K,L)*D(L,L)' = E(K,L) - R(K,L)

         #  where
         #                        M
         #            R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L)' +
         #                      J=K+1
         #                        M              N
         #                       SUM { A(K,J) * SUM [X(J,I)*B(L,I)'] } +
         #                       J=K           I=L+1

         #                       M
         #                    { SUM [C(K,J)*X(J,L)] } * D(L,L)' +
         #                     J=K+1
         #                       M              N
         #                      SUM { C(K,J) * SUM [X(J,I)*D(L,I)'] }.
         #                      J=K           I=L+1
         # """
         for l = n:-1:1
             il1 = l+1:n
             for k = m:-1:1
                 y = E[k,l]
                 if k < m
                    ta = ZERO
                    tc = ZERO
                    for ir = k+1:m
                        ta += A[k,ir]*E[ir,l]
                        tc += C[k,ir]*E[ir,l]
                    end
                    y -= (ta*B[l,l]'+tc*D[l,l]')
                 end
                 if l < n
                    ta = ZERO
                    tc = ZERO
                    for ir = il1
                        ta += E[k,ir]*B[l,ir]'
                        tc += E[k,ir]*D[l,ir]'
                    end
                    WB[k] = ta
                    WD[k] = tc
                    for ic = k:m
                        y -= (A[k,ic]*WB[ic]+C[k,ic]*WD[ic])
                    end
                 end
                 E[k,l] = y/(B[l,l]'*A[k,k]+D[l,l]'*C[k,k])
                 isfinite(E[k,l]) || throw("ME:SingularException: A-λC and D'+λB' have common or close eigenvalues")
             end
         end
   elseif adjAC && !adjBD
      # """
      # The (K,L)th element of X is determined starting from the
      # upper-left corner column by column by

      # A(K,K)'*X(K,L)*B(L,L) + C(K,K)'*X(K,L)*D(L,L) = E(K,L) - R(K,L),

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(I,L)] } +
      #                       J=1             I=1

      #                       K-1
      #                     { SUM [C(J,K)'*X(J,L)] } * D(L,L) +
      #                       J=1
      #                        K              L-1
      #                       SUM C(J,K)' * { SUM [X(J,I)*D(I,L)] }.
      #                       J=1             I=1
      # """
      for l = 1:n
          il1 = 1:l-1
          for k = 1:m
              y = E[k,l]
              if k > 1
                 ta = ZERO
                 tc = ZERO
                 for ir = 1:k-1
                     ta += A[ir,k]'*E[ir,l]
                     tc += C[ir,k]'*E[ir,l]
                 end
                 y -= (ta*B[l,l]+tc*D[l,l])
              end
              if l > 1
                 ta = ZERO
                 tc = ZERO
                 for ir = il1
                     ta += E[k,ir]*B[ir,l]
                     tc += E[k,ir]*D[ir,l]
                 end
                 WB[k] = ta
                 WD[k] = tc
                 for ic = 1:m
                     y -= (A[ic,k]'*WB[ic]+C[ic,k]'*WD[ic])
                 end
              end
              E[k,l] = y/(B[l,l]*A[k,k]'+D[l,l]*C[k,k]')
              isfinite(E[k,l]) || throw("ME:SingularException: A'-λC' and D+λB have common or close eigenvalues")
          end
      end
   elseif adjAC && adjBD
      # """
      # The (K,L)th element of X is determined starting from
      # upper-rght corner column by column by

      #       A(K,K)'*X(K,L)*B(L,L)' + C(K,K)'*X(K,L)*D(L,L)' = E(K,L) - R(K,L)

      # where
      #                       K-1
      #            R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM A(J,K)' * { SUM [X(J,I)*B(L,I)'] }+
      #                       J=1            I=L+1

      #                       K-1
      #                     { SUM [C(J,K)'*X(J,L)] } * D(L,L)' +
      #                       J=1
      #                        K               N
      #                       SUM C(J,K)' * { SUM [X(J,I)*D(L,I)'] }.
      #                       J=1            I=L+1
      # """
      for l = n:-1:1
          il1 = l+1:n
          for k = 1:m
              y = E[k,l]
              if k > 1
                 ta = ZERO
                 tc = ZERO
                 for ir = 1:k-1
                     ta += A[ir,k]'*E[ir,l]
                     tc += C[ir,k]'*E[ir,l]
                 end
                 y -= (ta*B[l,l]'+tc*D[l,l]')
               end
               if l < n
                 ic = 1:m
                 ta = ZERO
                 tc = ZERO
                 for ir = il1
                     ta += E[k,ir]*B[l,ir]'
                     tc += E[k,ir]*D[l,ir]'
                 end
                 WB[k] = ta
                 WD[k] = tc
                 for ic = 1:m
                     y -= (A[ic,k]'*WB[ic]+C[ic,k]'*WD[ic])
                 end
               end
              E[k,l] = y/(B[l,l]'*A[k,k]'+D[l,l]'*C[k,k]')
              isfinite(E[k,l]) || throw("ME:SingularException: A'-λC' and D'+λB' have common or close eigenvalues")
          end
      end
   end
   return E
end
"""
    (X,Y) = sylvsyss!(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F,

where `(A,D)`, `(B,E)` are pairs of square matrices of the same size in generalized Schur forms.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues. The computed
solution `(X,Y)` is contained in `(C,F)`.

_Note:_ This is an enhanced interface to the `LAPACK.tgsyl!` function to also cover the case when
`A`, `B`, `D` and `E` are real matrices and `C` and `F` are complex matrices.
"""
function sylvsyss!(A::T1, B::T1, C::T1, D::T1, E::T1, F::T1) where {T<:BlasFloat,T1<:Matrix{T}}
   """
   This is a wrapper to the LAPACK.tgsyl! function with `trans = 'N'`.
   """
   C, F, scale =  tgsyl!('N',A,B,C,D,E,F)
   return rmul!(C,inv(scale)), rmul!(F,inv(-scale))
end
"""
    (X,Y) = dsylvsyss!(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

    A'X + D'Y = C
    XB' + YE' = F,

where `(A,D)`, `(B,E)` are pairs of square matrices of the same size in generalized Schur forms.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues. The computed
solution `(X,Y)` is contained in `(C,F)`.
"""
function dsylvsyss!(A::T1, B::T1, C::T1, D::T1, E::T1, F::T1) where {T<:BlasFloat,T1<:Matrix{T}}
   """
   This is an interface to the LAPACK.tgsyl! function with `trans = 'T' or `trans = 'C'`.
   """
   # MF = -F
   # E, F, scale =  tgsyl!(T <: Complex ? 'C' : 'T', A, B, C, D, E, MF)
   # F = MF
   # return rmul!(C[:,:],inv(scale)), rmul!(F[:,:],inv(scale))
   F = -F
   C, F, scale =  tgsyl!(T <: Complex ? 'C' : 'T', A, B, C, D, E, F)
   return rmul!(C[:,:],inv(scale)), rmul!(F[:,:],inv(scale))
end


text:
abstract type AbstractTriangulation end

abstract type AbstractTriangulationPartition{D, T} <: AbstractTriangulation end
abstract type AbstractMutableTriangulationPartition{D, T} <: AbstractTriangulationPartition{D, T} end

abstract type AbstractTriangulationPartitionFull{D, T} <: AbstractTriangulationPartition{D, T} end
abstract type AbstractMutableTriangulationPartitionFull{D, T} <: AbstractTriangulationPartitionFull{D, T} end

Base.eltype(::AbstractTriangulationPartition{D,T}) where {D,T} = T

####################
# Various
####################
data(tp::AbstractTriangulationPartition) = tp.data
simplexindices(tp::AbstractTriangulationPartition) = tp.simplexindices
getsimplex(tp::AbstractTriangulationPartition, i) = tp.simplexindices[i]
nsimplices(tp::AbstractTriangulationPartition) = length(tp.simplexindices)

dimension(tp::AbstractTriangulationPartition) = dimension(tp.simplexindices)

####################
# Pretty printing
####################
function summarise(tp::AbstractTriangulationPartition)
    _type = typeof(tp)
    n_simplices = nsimplices(tp)
    D = dimension(tp)
    summary = "$_type with $n_simplices simplices"
end

Base.show(io::IO, tp::AbstractTriangulationPartition) = println(io, summarise(tp))

####################
# Getting simplices
####################

"""
Materialize the simplices of a triangulation.
"""
function simplices(tpp::AbstractTriangulationPartition)

end

function getpoints(tpp::AbstractTriangulationPartition) end
function getsimplices(tpp::AbstractTriangulationPartition) end

export
AbstractTriangulationPartition,
AbstractMutableTriangulationPartition,
AbstractTriangulationPartitionFull,
AbstractMutableTriangulationPartitionFull,
data,
simplexindices,
getsimplex,
nsimplices,
dimension,
simplices,

getpoints,
getsimplices


text:
@testset "KNN" begin
    @testset "Default" begin
        tester = ImputorTester(KNN)
        test_hashing(tester)
        test_equality(tester)
        test_matrix(tester)
        test_axisarray(tester)
        test_nameddimsarray(tester)
        test_keyedarray(tester)
    end
    @testset "Iris" begin
        # Reference
        # P. Schimitt, et. al
        # A comparison of six methods for missing data imputation
        iris = Impute.dataset("test/table/iris") |> DataFrame
        data = Array(iris[:, [:SepalLength, :SepalWidth, :PetalLength, :PetalWidth]])

        @testset "MCAR $r" for r in (0.15, 0.25, 0.35)
            # Original authors said they ran the test 1000 times
            results = map(1:1000) do i
                X = add_missings(data, r)
                knn_imputed = impute(copy(X), Impute.KNN(; k=3); dims=:cols)
                sub_imputed = impute(copy(X), Substitute(); dims=:cols)

                return nrmsd(data, knn_imputed), nrmsd(data, sub_imputed)
            end

            knn_nrmsd = first.(results)
            sub_nrmsd = last.(results)

            # Test that knn is better on average
            @test mean(knn_nrmsd) < mean(sub_nrmsd)
            # Test that the difference is "significant".
            @test pvalue(UnequalVarianceTTest(knn_nrmsd, sub_nrmsd)) < 0.05

            # test type stability
            @test typeof(impute(add_missings(data, r), Impute.KNN(; k=3); dims=:cols)) == Matrix{Union{Float64, Missing}}
            @test typeof(impute(add_missings(data, r), Substitute(); dims=:cols)) == Matrix{Union{Float64, Missing}}
        end
    end

    # Test a case where we expect kNN to perform well (e.g., many variables, )
    @testset "Data match" begin
        data = mapreduce(hcat, 1:1000) do i
            seeds = [sin(i), cos(i), tan(i), atan(i)]
            mapreduce(vcat, combinations(seeds)) do args
                [
                    +(args...),
                    *(args...),
                    +(args...) * 100,
                    +(abs.(args)...),
                    (+(args...) * 10) ^ 2,
                    (+(abs.(args)...) * 10) ^ 2,
                    log(+(abs.(args)...) * 100),
                    +(args...) * 100 + rand(-10:0.1:10),
                ]
            end
        end

        X = add_missings(data')
        num_tests = 100

        knn_nrmsd, mean_nrmsd = 0.0, 0.0

        for i = 1:num_tests
            knn_imputed = impute(copy(X), Impute.KNN(; k=4); dims=:cols)
            mean_imputed = impute(copy(X), Substitute(); dims=:cols)

            knn_nrmsd = ((i - 1) * knn_nrmsd + nrmsd(data', knn_imputed)) / i
            mean_nrmsd = ((i - 1) * mean_nrmsd + nrmsd(data', mean_imputed)) / i
        end

        # @show knn_nrmsd mean_nrmsd
        @test knn_nrmsd < mean_nrmsd
        # test type stability
        @test typeof(X) == typeof(impute(copy(X), Impute.KNN(; k=4); dims=:cols))
        @test typeof(X) == typeof(impute(copy(X), Substitute(); dims=:cols))
    end
end