text:
abstract type Register end

mutable struct CircuitReference
    reference::Ref{Circuit}
    function CircuitReference()
        new()
    end
end

"""
    QRegister{N}

Quantum Register of `N` qubits.

"""
struct QRegister <: Register
    # TODO: Create function to read QRegister using ind
    "size of register"
    N::Int64
    "vector of quantum wire numbers in this quantum register"
    ind::Vector{Int}
    "ref to circuit"
    circuit::CircuitReference

    # function QRegister(N::Int64, ind::Vector{Int})
    #     new(N, ind)
    # end
end

function (qreg::QRegister)(indices::Int...)
    qreg, qreg[[indices...]]
end

function (qreg::QRegister)(indices::Array{Int,1})
    qreg, qreg[indices]
end

# TODO: Reimplement CRegister smartly
# """
#     CRegister{N}

# Classical register of `N` bits.
# """
# struct CRegister <: Register
#     # TODO: Create function to read QRegister using ind
#     "size of register"
#     N::Int
#     "vector of quantum wire numbers in this quantum register"
#     ind::Vector{Int}
#     "bit array of bits in this classical register"
#     n::BitArray{1}
#     "ref to circuit"
#     circuit::CircuitReference

#     # function CRegister(N::Int64, ind::Vector{Int}, n::BitArray{1})
#     #     new(N, ind, n)
#     # end
# end

"""
    qreg(n::Int)

constructs quantum register of size `n`
"""
function qreg(n::Int)
    QRegister(n, zeros(Int, n), CircuitReference())
end

# """
#     creg(n::Int)

# constructs classical register of size `n`
# """
# function creg(n::Int)
#     CRegister(n, zeros(Int, n), BitArray(undef, (n)), CircuitReference())
# end

# make Registers indexable
function Base.getindex(r::Register, i::Int)
    !any(r.ind .== 0) || error("Register object has yet to be used in a Circuit object")
    minimum(i) >= 1 || throw(BoundsError(r.ind, minimum(i)))
    maximum(i) <= length(r.ind) || throw(BoundsError(r.ind, maximum(i)))

    return r.ind[i]
end

function Base.getindex(r::Register, i::Array{Int,1})
    !any(r.ind .== 0) || error("Register object has yet to be used in a Circuit object")
    minimum(i) >= 1 || throw(BoundsError(r.ind, minimum(i)))
    maximum(i) <= length(r.ind) || throw(BoundsError(r.ind, maximum(i)))

    return r.ind[i]
end

function Base.getindex(r::Register, ::Colon)
    !any(r.ind .== 0) || error("Register object has yet to be used in a Circuit object")

    return r.ind[:]
end

function Base.getindex(r::Register, ur::UnitRange{Int64})
    !any(r.ind .== 0) || error("Register object has yet to be used in a Circuit object")

    return r.ind[ur]
end

# """
#     reg_check(c::Qaintessent.CRegister, val::Int)

# create expression to verify that classical register is equal to integer value `val`
# """
# function reg_check(c::Qaintessent.CRegister, val::Int)
#     r = Ref(c.n)
#     return :(Qaintessent.binary_to_int($r[]) == $val)
# end

# """
#     set_creg!(reg::CRegister, index::Int, value::Bool) where {N}

# Set the classical register `reg` bit `index` to value `value`
# """
# function set_creg!(reg::CRegister, index::Int, value::Bool)
#     index = length(reg.n) - index + 1
#     reg.n[index] = value
# end

# """
#     set_creg!(reg::CRegister, value::Int) where {N}

# Set the classical register `reg` to value `value`
# """
# function set_creg!(reg::CRegister, value::Int)
#     l = length(reg.n)
#     log2(value) <= l || error("Unable to store integer value " * string(value) * " in BitArray of size " * string(l))
#     reg.n .= binary_digits(l, value)
# end

# num_wires(creg::CRegister) = creg.N
num_wires(qreg::QRegister) = qreg.N

"""
    Circuit(regs::Register...)

helper function, creates Circuit from multiple quantum or classical registers
"""
function Circuit(regs::Register...)
    N = 0
    j = 1
    for i in 1:length(regs)
        !isdefined(regs[i].circuit, :reference) || error("Input Register $i is already used in another circuit.")
        N += regs[i].N
        regs[i].ind .= j:regs[i].N+j-1
        j = regs[i].N+j
    end
    c = Circuit{N}()
    for i in 1:length(regs)
        regs[i].circuit.reference = c
    end
    return c
end

function add!(g::AbstractGate, q::QRegister)
    isdefined(q.circuit, :reference) || error("Register object has yet to be used in a Circuit object")
    if num_wires(g) == num_wires(q)
        append!(q.circuit.reference[], circuit_gate(Tuple(q[:]), g))
        return
    end
    num_wires(g) == 1 || error("AbstractGate $g is applied to $(num_wires(g)) wires, quantum register of size $(num_wires(q)) provided. Unable to determine gates to be applied.")

    for i in 1:num_wires(q)
        append!(q.circuit.reference[], circuit_gate(q[i], g))
    end
end

function add!(g::AbstractGate, target_wires::Tuple{QRegister, Array{Int}}...)
    circuit_length = num_wires(g)
    circuit_references = getfield.(target_wires, 1)

    wires = reduce(hcat, getfield.(target_wires, 2))
    length(wires) == circuit_length || error("AbstractGate $g is applied to $circuit_length wires. However, $(length(wires)) wires provided")
    for ref in circuit_references[2:end]
        isdefined(ref.circuit, :reference) || error("Register object has yet to be used in a Circuit object")
        ref.circuit.reference[] === circuit_references[1].circuit.reference[] || error("Targetted Registers belong to different Circuit objects")
    end
    append!(circuit_references[1].circuit.reference[], CircuitGate(Tuple(wires), g))
    return
end

function add_control!(g::AbstractGate, target_q::QRegister, control_q::QRegister)
    num_wires(g) == 1 || error("Helper function for controlled gates across entire registers only supports 1 qubit gates")
    isdefined(target_q.circuit, :reference) || error("Target quantum register is not included in any circuit")
    isdefined(control_q.circuit, :reference) || error("Control quantum register is not included in any circuit")
    target_q.circuit.reference[] === control_q.circuit.reference[] || error("Target quantum register and control quantum register are used in different circuits")

    num_wires(target_q) == num_wires(control_q) || error("Control register and Target register have differing number of qubits, unable to determine gates to be applied")
    for i in 1:num_wires(control_q)
        append!(target_q.circuit.reference[], circuit_gate(target_q.ind[i], g, control_q.ind[i]))
    end
end

function add_control!(g::AbstractGate, target_wires::Tuple{QRegister, Array{Int}}, control_wires::Tuple{QRegister, Array{Int}})
    isdefined(target_wires[1].circuit, :reference) || error("Target quantum register is not included in any circuit")
    isdefined(control_wires[1].circuit, :reference) || error("Control quantum register is not included in any circuit")
    target_wires[1].circuit.reference[] === control_wires[1].circuit.reference[] || error("Target quantum register and control quantum register are used in different circuits")

    circuit_length = num_wires(g)
    num_wires(g) == length(target_wires[2]) || error("Gate affecting $(num_wires(g)) qubits applied to $(length(target_wires[2])) qubits")
    append!(target_wires[1].circuit.reference[], circuit_gate((target_wires[2]...,), g, control_wires[2]...))
end


text:
function doubling{T<:AbstractFloat,S<:Int}(a::Array{T,2},b::Array{T,2},g::Array{T,2},tol::T,maxiters::S)

  retcode = 0
  j = 1
  len = Inf

  gn = zeros(g)

  while len > tol && j <= maxiters
    gn = g+a*g*b
    len = maxabs(gn-g)
    a = a*a
    b = b*b
    g = copy(gn)
    j += 1
  end

  if j > maxiters
    retcode = 1
  end

  return gn,retcode

end

ext: jl


text:
using Interpolations
using BenchmarkTools

function interpolate!(a_hi, a_lo, x_hi, x_lo)
    interp_a = interpolate((x_lo, x_lo, x_lo), a_lo, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    a_hi[:, :, :] .= interp_a(x_hi, x_hi, x_hi)
end

n_hi = 256
n_lo = 128

dx_hi = 1//n_hi
dx_lo = 1//n_lo

x_hi = 1//2*dx_hi:dx_hi:1
x_lo = -1//2*dx_lo:dx_lo:1+1//2*dx_lo

a0_lo = rand(n_lo+2, n_lo+2, n_lo+2)
a0_hi = zeros(n_hi, n_hi, n_hi)

a1_lo = rand(n_lo+2, n_lo+2, n_lo+2)
a1_hi = zeros(n_hi, n_hi, n_hi)

a2_lo = rand(n_lo+2, n_lo+2, n_lo+2)
a2_hi = zeros(n_hi, n_hi, n_hi)

a3_lo = rand(n_lo+2, n_lo+2, n_lo+2)
a3_hi = zeros(n_hi, n_hi, n_hi)

@btime begin
    @sync begin
        Threads.@spawn interpolate!(a0_hi, a0_lo, x_hi, x_lo)
        Threads.@spawn interpolate!(a1_hi, a1_lo, x_hi, x_lo)
        Threads.@spawn interpolate!(a2_hi, a2_lo, x_hi, x_lo)
        Threads.@spawn interpolate!(a3_hi, a3_lo, x_hi, x_lo)
    end
end


text:
# Make sure we're running from the examples directory.
cd(joinpath(@__DIR__(), ".."))

using Eratosthenes # Bring in the top-level sim environment.
using Plots # For plotting (Pkg.add("Plots"); Pkg.add("PyPlot");)
import HDF5 # For loading the log file (Pkg.add("HDF5");)

# Run the sim on the simplest little scenario.
scenario = simulate("scenarios/min-scenario.yaml")

# Open the logs and plot some things.
if !isempty(scenario.sim.log_file)

    # Read some data from the logs. This do block takes care of
    # error-handling, so we can just read whatever we want and return it.
    t, q_BI, ω_BI_B = HDF5.h5open(scenario.sim.log_file, "r") do logs
        (read(logs, "/cube1/body/state/time"),
         read(logs, "/cube1/body/state/data/q_BI"),
         read(logs, "/cube1/body/state/data/ω_BI_B"))
    end

    # Scalars are logged as 1-by-n (don't know why). Convert to just n for
    # the sake of plotting.
    t = dropdims(t, dims=1)

    # Choose a friendly plotting package. PyPlot behaves nicely.
    pyplot()

    # Define each plot that we'll need.
    p1 = plot(t, transpose(q_BI),
              label  = ["q1" "q2" "q3" "q4"],
              ylabel = "Attitude Quaternion")
    p2 = plot(t, 180/π * transpose(ω_BI_B),
              label  = ["ω1" "ω2" "ω3"],
              ylabel = "Rotation Rate (deg/s)")

    # Save the individual plots before showing them together.
    savefig(p1, "out/sim-results-q.png")
    savefig(p2, "out/sim-results-ω.png")

    # Display those plots as subplots using the default "2" layout. Put any
    # common descriptors here.
    display(plot(p1, p2,
                 layout = 2,
                 xlabel = "Time (s)"))

end # plotting


text:
# Continuous Lyapunov equations
"""
    X = lyapc(A, C)

Compute `X`, the solution of the continuous Lyapunov equation

      AX + XA' + C = 0,

where `A` is a square real or complex matrix and `C` is a square matrix.
`A` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
The solution `X` is symmetric or hermitian if `C` is a symmetric or hermitian.

The following particular cases are also adressed:

    X = lyapc(α*I,C)  or  X = lyapc(α,C)

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

    x = lyapc(α,γ)

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

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

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

julia> X = lyapc(A, C)
2×2 Array{Float64,2}:
  0.5  -0.5
 -0.5   0.25

julia> A*X + X*A' + C
2×2 Array{Float64,2}:
 -8.88178e-16   2.22045e-16
  2.22045e-16  -4.44089e-16
```
"""
function lyapc(A::AbstractMatrix, C::AbstractMatrix)
   """
   The Bartels-Steward Schur form based method 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.
   """
   n = LinearAlgebra.checksquare(A)
   LinearAlgebra.checksquare(C) == n ||
      throw(DimensionMismatch("C must be a square matrix of dimension $n"))

   adj = isa(A,Adjoint)
   her = ishermitian(C)

   T2 = promote_type(eltype(A), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce A to Schur form and transform C
   if adj
      AS, Q = schur(A.parent)
   else
      AS, Q = schur(A)
   end

   #X = Q'*C*Q
   if her
      X = utqu(C,Q)
      lyapcs!(AS, X, adj = adj)
      #X <- Q*X*Q'
      utqu!(X,Q')
      return X
   else
      X = Q' * C * Q
      adj ? (sylvcs!(AS, AS, X, adjA = true)) : (sylvcs!(AS, AS, X, adjB = true))
      return rmul!(Q * X * Q',-1)
   end
end
# (α+α')X  + C = 0
lyapc(A::UniformScaling, C::AbstractMatrix) = -C/(A+A')
lyapc(A::Union{Real,Complex}, C::AbstractMatrix) = real(A) == 0 ? throw(SingularException(1)) : -C/(A+A')
# (α+α')x + γ = 0
lyapc(A::Union{Real,Complex}, C::Union{Real,Complex})  = real(A) == 0 ? throw(SingularException(1)) : -C/(A+A')
"""
    X = lyapc(A, E, C)

Compute `X`, the solution of the generalized continuous Lyapunov equation

     AXE' + EXA' + C = 0,

where `A` and `E` are square real or complex matrices and `C` is a square matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
The solution `X` is symmetric or hermitian if `C` is a symmetric or hermitian.

The following particular cases are also adressed:

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

Solve the matrix equation `AXβ' + βXA' + C = 0`.

    X = lyapc(α*I,E,C)  or  X = lyapc(α,E,C)

Solve the matrix equation `αXE' + EXα' + C = 0`.

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

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

    x = lyapc(α,β,γ)

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

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

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

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

julia> X = lyapc(A, E, C)
2×2 Array{Float64,2}:
 -2.5   2.5
  2.5  -2.25

julia> A*X*E' + E*X*A' + C
2×2 Array{Float64,2}:
 -5.32907e-15  -2.66454e-15
 -4.44089e-15   0.0
```
"""
function lyapc(A::AbstractMatrix, E::AbstractMatrix, C::AbstractMatrix)
   """
   The extension of the Bartels-Steward method based on the generalized Schur form
   is employed.

   Reference:
   T. Penzl. Numerical solution of generalized Lyapunov equations.
   Adv. Comput. Math., 8:33–48, 1998.
   """
   n = LinearAlgebra.checksquare(A)
   LinearAlgebra.checksquare(C) == n ||
      throw(DimensionMismatch("C must be a square matrix of dimension $n"))
   isequal(E,I) && size(E,1) == n && (return lyapc(A, C))
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a square matrix of dimension $n"))

   adjA = isa(A,Adjoint)
   adjE = isa(E,Adjoint)
   if adjA && !adjE
      A = copy(A)
      adjA = false
   elseif !adjA && adjE
      E = copy(E)
      adjE = false
   end

   adj = adjA & adjE
   her = ishermitian(C)

   T2 = promote_type(eltype(A), eltype(E), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(E) == T2 || (adj ? E = convert(Matrix{T2},E.parent)' : E = convert(Matrix{T2},E))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce (A,E) to generalized Schur form and transform C
   # (as,es) = (q'*A*z, q'*E*z)
   if adj
      as, es, z, q = schur(A.parent,E.parent)
   else
      as, es, q, z = schur(A,E)
   end
   if her
      #x = q'*C*q
      x = utqu(C,q)
      lyapcs!(as,es,x, adj = adj)
      #x = z*x*z'
      utqu!(x,z')
      return x
   else
      x = q' * C * q
      adj ? (gsylvs!(as,es,es,as,x,adjAC = true,DBSchur = true)) : (gsylvs!(as,es,es,as,x,adjBD = true,DBSchur = true))
      return rmul!(z*x*z', -1)
   end
end
# Aβ'X + XA'β + C = 0
lyapc(A::AbstractMatrix, E::Union{Real,Complex,UniformScaling}, C::AbstractMatrix) =
      isa(A,Adjoint) ? lyapc((A.parent*E')',C) : lyapc(A*E',C)
# αXE' + EXα' + C = 0
lyapc(A::Union{Real,Complex,UniformScaling}, E::AbstractMatrix, C::AbstractMatrix) =
      isa(E,Adjoint) ? lyapc((E.parent*A')',C) : lyapc(A'*E,C)
# (αβ'+α'β)X + C = 0
lyapc(A::UniformScaling, E::UniformScaling, C::AbstractMatrix) = -(A*E'+A'*E)\C
lyapc(A::Union{Real,Complex}, E::Union{Real,Complex}, C::AbstractMatrix) = real(A*E') == 0 ? throw(SingularException(1)) : -C/(A*E'+A'*E)
# (αβ'+α'β)X  + γ = 0
lyapc(A::Union{Real,Complex}, E::Union{Real,Complex}, C::Union{Real,Complex}) = real(A*E') == 0 ? throw(SingularException(1)) : -C/(A*E'+A'*E)
# Discrete Lyapunov equations
"""
    X = lyapd(A, C)

Compute `X`, the solution of the discrete Lyapunov equation

       AXA' - X + C = 0,

where `A` is a square real or complex matrix and `C` is a square matrix.
`A` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The solution `X` is symmetric or hermitian if `C` is a symmetric or hermitian.

The following particular cases are also adressed:

    X = lyapd(α*I,C)  or  X = lyapd(α,C)

Solve the matrix equation `(αα'-1)X + C = 0`.

    x = lyapd(α,γ)

Solve the equation `(αα'-1)x + γ = 0`.

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

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

julia> X = lyapd(A, C)
2×2 Array{Float64,2}:
  0.2375  -0.2125
 -0.2125   0.1375

julia> A*X*A' - X + C
2×2 Array{Float64,2}:
 5.55112e-16  6.66134e-16
 2.22045e-16  4.44089e-16
```
"""
function lyapd(A::AbstractMatrix, C::AbstractMatrix)
   """
   The discrete analog of the Bartels-Steward method based on the Schur form
   is employed.

   Reference:
   G. Kitagawa. An Algorithm for solving the matrix equation X = F X F' + S,
   International Journal of Control, 25:745-753, 1977.
   """
   n = LinearAlgebra.checksquare(A)
   LinearAlgebra.checksquare(C) == n ||
       throw(DimensionMismatch("C must be a square matrix of dimension $n"))

   adj = isa(A,Adjoint)
   her = ishermitian(C)

   T2 = promote_type(eltype(A), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce A to Schur form and transform C
   if adj
      AS, Q = schur(A.parent)
   else
      AS, Q = schur(A)
   end
   #X = Q'*C*Q
   if her
      X = utqu(C,Q)
      lyapds!(AS, X, adj = adj)
      #X <- Q*X*Q'
      utqu!(X,Q')
      return X
   else
      #X = rmul!( Q' * C * Q, -1)
      X = Q' * C * Q
      adj ? (sylvds!(-AS, AS, X, adjA = true)) : (sylvds!(-AS, AS, X, adjB = true))
      return Q * X * Q'
   end
end
# (αα'-1)X + C = 0
lyapd(A::UniformScaling, C::AbstractMatrix) = (I-A'*A)\C
lyapd(A::Union{Real,Complex}, C::AbstractMatrix) = A*A' == 1 ? throw(SingularException(1)) : C/(1-A'*A)
# (αα'-1)x + γ = 0
lyapd(A::Union{Real,Complex}, C::Union{Real,Complex}) = A*A' == 1 ? throw(SingularException(1)) : C/(one(C)-A'*A)
"""
    X = lyapd(A, E, C)

Compute `X`, the solution of the generalized discrete Lyapunov equation

         AXA' - EXE' + C = 0,

where `A` and `E` are square real or complex matrices and `C` is a square matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The solution `X` is symmetric or hermitian if `C` is a symmetric or hermitian.

The following particular cases are also adressed:

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

Solve the matrix equation `AXA' - βXβ' + C = 0`.

    X = lyapd(α*I,E,C)  or  X = lyapd(α,E,C)

Solve the matrix equation `αXα' - EXE' + C = 0`.

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

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

    x = lyapd(α,β,γ)

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

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

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

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

julia> X = lyapd(A, E, C)
2×2 Array{Float64,2}:
  1.775  -1.225
 -1.225   0.775

julia> A*X*A' - E*X*E' + C
2×2 Array{Float64,2}:
 -2.22045e-16  -4.44089e-16
 -1.33227e-15   1.11022e-15
```
"""
function lyapd(A::AbstractMatrix, E::AbstractMatrix, C::AbstractMatrix)
   """
   The extension of the Bartels-Steward method based on the generalized Schur form
   is employed.

   Reference:
   T. Penzl. Numerical solution of generalized Lyapunov equations.
   Adv. Comput. Math., 8:33–48, 1998.
   """
   n = LinearAlgebra.checksquare(A)
   LinearAlgebra.checksquare(C) == n ||
      throw(DimensionMismatch("C must be a square matrix of dimension $n"))
   isequal(E,I) && size(E,1) == n && (return lyapd(A, C))
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a square matrix of dimension $n"))

   adjA = isa(A,Adjoint)
   adjE = isa(E,Adjoint)
   if adjA && !adjE
      A = copy(A)
      adjA = false
   elseif !adjA && adjE
      E = copy(E)
      adjE = false
   end

   adj = adjA & adjE
   her = ishermitian(C)

   T2 = promote_type(eltype(A), eltype(E), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(E) == T2 || (adj ? E = convert(Matrix{T2},E.parent)' : E = convert(Matrix{T2},E))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce (A,E) to generalized Schur form and transform C
   # (as,es) = (q'*A*z, q'*E*z)
   if adj
      as, es, z, q = schur(A.parent,E.parent)
   else
      as, es, q, z = schur(A,E)
   end
   if her
      #x = q'*C*q
      x = utqu(C,q)
      lyapds!(as,es,x, adj = adj)
      #x = z*x*z'
      utqu!(x,z')
      return x
   else
      x = q' * C * q
      gsylvs!(as, as, -es, es, x, adjAC = adj, adjBD = !adj)
      return rmul!(z*x*z', -1)
   end
end
# AXA' - Xββ' + C = 0
lyapd(A::AbstractMatrix, E::Union{Real,Complex,UniformScaling}, C::AbstractMatrix) =
      isa(A,Adjoint) ? lyapd((A.parent/E)',C/(E*E')) : lyapd(A/E,C/(E*E'))
# αXα' - EXE' + C = 0
lyapd(A::Union{Real,Complex,UniformScaling}, E::AbstractMatrix, C::AbstractMatrix) =
      isa(E,Adjoint) ? lyapd((E.parent/A)',C/(-A*A')) : lyapd(E/A,C/(-A*A'))
# (α'α-β'β)X + C = 0
lyapd(A::UniformScaling, E::UniformScaling, C::AbstractMatrix) = C/(E'*E-A'*A)
lyapd(A::Union{Real,Complex}, E::Union{Real,Complex}, C::AbstractMatrix) = A*A'== E*E' ? throw(SingularException(1)) : C/(E'*E-A'*A)
# (α'α-β'β)X + γ  = 0
lyapd(A::Union{Real,Complex}, E::Union{Real,Complex}, C::Union{Real,Complex}) = A*A'== E*E' ? throw(SingularException(1)) : C/(E'*E-A'*A)
"""
    lyapcs!(A,C;adj = false)

Solve the continuous Lyapunov matrix equation

                op(A)X + Xop(A)' + C = 0,

where `op(A) = A` if `adj = false` and `op(A) = A'` if `adj = true`.
`A` is a square real matrix in a real Schur form, or a square complex matrix in a
complex Schur form and `C` is a symmetric or hermitian matrix.
`A` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
`C` contains on output the solution `X`.
"""
function lyapcs!(A::AbstractMatrix{T1},C::AbstractMatrix{T1}; adj = false) where {T1<:BlasReal}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && issymmetric(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric matrix"))

   ONE = one(T1)
   ZERO = zero(T1)

   # determine the structure of the real Schur form
   ba, p = sfstruct(A)

   #W = Array{T1,2}(I,2,2)
   Xw = Matrix{T1}(undef,4,4)
   Yw = Vector{T1}(undef,4)
   if adj
      # """
      # The (K,L)th block of X is determined starting from
      # upper-left corner column by column by

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

      # where
      #          K-1                    L-1
      # R(K,L) = SUM [A(I,K)'*X(I,L)] + SUM [X(K,J)*A(J,L)].
      #          I=1                    J=1
      # """
      j = 1
      for ll = 1:p
          dl = ba[ll]
          l = j:j+dl-1
          i = j
          for kk = ll:p
              dk = ba[kk]
              k = i:i+dk-1
              y = view(C,k,l)
              if kk > ll
                 ia = j:i-1
                 # y += A[ia,k]'*C[ia,l]
                 mul!(y,transpose(view(A,ia,k)),view(C,ia,l),ONE,ONE)
               end
               if i == j
                  lyapc2!(adj,y,dk,view(A,k,k),Xw,Yw)
               else
                  lyapcsylv2!(adj,y,dk,dl,view(A,k,k),view(A,l,l),Xw,Yw)
                  transpose!(view(C,l,k),y)
               end
               i += dk
          end
          j += dl
          if j <= n
             for jr = j:n
                for ir = jr:n
                    for lll = l
                        C[ir,jr] += C[ir,lll]*A[lll,jr] + A[lll,ir]*C[jr,lll]
                    end
                    C[jr,ir] = C[ir,jr]
                end
             end
          end
       end
   else
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-right corner column by column by

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

      # where
      #           N                     N
      # R(K,L) = SUM [A(K,I)*X(I,L)] + SUM [X(K,J)*A(L,J)'].
      #         I=K+1                 J=L+1
      # """
      j = n
      for ll = p:-1:1
          dl = ba[ll]
          l = j-dl+1:j
          i = j
          for kk = ll:-1:1
              dk = ba[kk]
              i1 = i-dk+1
              k = i1:i
              y = view(C,k,l)
              if kk < ll
                 ia = i+1:j
                 #  y += A[k,ia]*C[ia,l]
                 mul!(y,view(A,k,ia),view(C,ia,l),ONE,ONE)
              end
              if i == j
                 lyapc2!(adj,y,dk,view(A,k,k),Xw,Yw)
              else
                 lyapcsylv2!(adj,y,dk,dl,view(A,k,k),view(A,l,l),Xw,Yw)
                 transpose!(view(C,l,k),y)
              end
              i -= dk
           end
           j -= dl
           if j >= 0
              for jr = 1:j
                 for ir = 1:jr
                     for lll = l
                         C[ir,jr] += C[ir,lll]*A[jr,lll] + A[ir,lll]*C[lll,jr]
                      end
                      C[jr,ir] = C[ir,jr]
                 end
              end
           end
       end
   end
end
function lyapc2!(adj,C::AbstractMatrix{T},na::Int,A::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
# speed and reduced allocation oriented implementation of a solver for 1x1 or 2x2 continuous Lyapunov equations
#      A'*X + X*A = -C if adj = true  -> R = kron(I,A')+kron(A',I) = (kron(I,A)+kron(A,I))'
#      A*X + X*A' = -C if adj = false -> R = kron(I,A)+kron(A,I)
   if na == 1
      rmul!(C,inv(-2*A[1,1]))
      any(!isfinite, C) && throw("ME:SingularException: A has eigenvalue(s) ≈ 0")
      return C
   end
   ZERO = zero(T)
   i1 = 1:3
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   @inbounds  Y[1] = -C[1,1]/2
   @inbounds  Y[2] = -C[2,1]
   @inbounds  Y[3] = -C[2,2]/2
   if adj
      # @inbounds R = [ A[1,1]    A[2,1]         0;
      #                 A[1,2]    A[1,1]+A[2,2]  A[2,1];
      #                 0         A[1,2]         A[2,2] ]
      @inbounds  R[1,1] = A[1,1]
      @inbounds  R[1,2] = A[2,1]
      @inbounds  R[1,3] = ZERO
      @inbounds  R[2,1] = A[1,2]
      @inbounds  R[2,2] = A[1,1]+A[2,2]
      @inbounds  R[2,3] = A[2,1]
      @inbounds  R[3,1] = ZERO
      @inbounds  R[3,2] = A[1,2]
      @inbounds  R[3,3] = A[2,2]
   else
      # @inbounds R = [ A[1,1]    A[1,2]         0;
      #                 A[2,1]    A[1,1]+A[2,2]  A[1,2];
      #                 0         A[2,1]         A[2,2] ]
      @inbounds  R[1,1] = A[1,1]
      @inbounds  R[1,2] = A[1,2]
      @inbounds  R[1,3] = ZERO
      @inbounds  R[2,1] = A[2,1]
      @inbounds  R[2,2] = A[1,1]+A[2,2]
      @inbounds  R[2,3] = A[1,2]
      @inbounds  R[3,1] = ZERO
      @inbounds  R[3,2] = A[2,1]
      @inbounds  R[3,3] = A[2,2]
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
   @inbounds C[1,1] = Y[1]
   @inbounds C[1,2] = Y[2]
   @inbounds C[2,1] = Y[2]
   @inbounds C[2,2] = Y[3]
   return C
end
function lyapcsylv2!(adj,C::AbstractMatrix{T},na::Int,nb::Int,A::AbstractMatrix{T},B::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
# speed and reduced allocation oriented implementation of a solver for 1x1 and 2x2 Sylvester equations
# encountered in solving continuous Lyapunov equations:
#      A'*X + X*B = -C if adj = true  -> R = kron(I,A')+kron(B',I) = (kron(I,A)+kron(B,I))'
#      A*X + X*B' = -C if adj = false -> R = kron(I,A)+kron(B,I)
   if na == 1 && nb == 1
      temp = A[1,1] + B[1,1]
      rmul!(C,inv(-temp))
      any(!isfinite, C) && throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
      return C
   end
   ZERO = zero(T)
   i1 = 1:na*nb
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   if adj
      if na == 1
         # @inbounds R = [ A[1,1]+B[1,1]    B[2,1];
         #                      B[1,2]      A[1,1]+B[2,2]]
         @inbounds  R[1,1] = A[1,1]+B[1,1]
         @inbounds  R[1,2] = B[2,1]
         @inbounds  R[2,1] = B[1,2]
         @inbounds  R[2,2] = A[1,1]+B[2,2]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # @inbounds R = [ A[1,1]+B[1,1]       A[2,1];
            #                  A[1,2]          A[2,2]+B[1,1] ]
            @inbounds  R[1,1] = A[1,1]+B[1,1]
            @inbounds  R[1,2] = A[2,1]
            @inbounds  R[2,1] = A[1,2]
            @inbounds  R[2,2] = A[2,2]+B[1,1]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # @inbounds R = [ A[1,1]+B[1,1]       A[2,1]       B[2,1]           0;
            #                 A[1,2]        A[2,2]+B[1,1]      0             B[2,1];
            #                 B[1,2]              0       A[1,1]+B[2,2]     A[2,1];
            #                   0               B[1,2]       A[1,2]       A[2,2]+B[2,2]]
            @inbounds  R[1,1] = A[1,1]+B[1,1]
            @inbounds  R[1,2] = A[2,1]
            @inbounds  R[1,3] = B[2,1]
            @inbounds  R[1,4] = ZERO
            @inbounds  R[2,1] = A[1,2]
            @inbounds  R[2,2] = A[2,2]+B[1,1]
            @inbounds  R[2,3] = ZERO
            @inbounds  R[2,4] = B[2,1]
            @inbounds  R[3,1] = B[1,2]
            @inbounds  R[3,2] = ZERO
            @inbounds  R[3,3] = A[1,1]+B[2,2]
            @inbounds  R[3,4] = A[2,1]
            @inbounds  R[4,1] = ZERO
            @inbounds  R[4,2] = B[1,2]
            @inbounds  R[4,3] = A[1,2]
            @inbounds  R[4,4] = A[2,2]+B[2,2]
            @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
   else
      if na == 1
         # @inbounds R = [ A[1,1]+B[1,1]       B[1,2];
         #                    B[2,1]       A[1,1] + B[2,2]]
         @inbounds  R[1,1] = A[1,1]+B[1,1]
         @inbounds  R[1,2] = B[1,2]
         @inbounds  R[2,1] = B[2,1]
         @inbounds  R[2,2] = A[1,1]+B[2,2]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # @inbounds R = [ A[1,1]+B[1,1]       A[1,2];
            #                 A[2,1]        A[2,2]+B[1,1]]
            @inbounds  R[1,1] = A[1,1]+B[1,1]
            @inbounds  R[1,2] = A[1,2]
            @inbounds  R[2,1] = A[2,1]
            @inbounds  R[2,2] = A[2,2]+B[1,1]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # @inbounds R = [ A[1,1]+B[1,1]       A[1,2]       B[1,2]         0;
            #                 A[2,1]      A[2,2]+B[1,1]       0         B[1,2];
            #                  B[2,1]            0      A[1,1]+B[2,2]     A[1,2];
            #                  0             B[2,1]        A[2,1]    A[2,2]+B[2,2]]
            @inbounds  R[1,1] = A[1,1]+B[1,1]
            @inbounds  R[1,2] = A[1,2]
            @inbounds  R[1,3] = B[1,2]
            @inbounds  R[1,4] = ZERO
            @inbounds  R[2,1] = A[2,1]
            @inbounds  R[2,2] = A[2,2]+B[1,1]
            @inbounds  R[2,3] = ZERO
            @inbounds  R[2,4] = B[1,2]
            @inbounds  R[3,1] = B[2,1]
            @inbounds  R[3,2] = ZERO
            @inbounds  R[3,3] = A[1,1]+B[2,2]
            @inbounds  R[3,4] = A[1,2]
            @inbounds  R[4,1] = ZERO
            @inbounds  R[4,2] = B[2,1]
            @inbounds  R[4,3] = A[2,1]
            @inbounds  R[4,4] = A[2,2]+B[2,2]
            @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
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
   C[:,:] = Y
   return C
end
function lyapcs!(A::AbstractMatrix{T1},C::AbstractMatrix{T1}; adj = false) where {T1<:BlasComplex}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n hermitian matrix"))

   if adj
      # """
      # The (K,L)th element of X is determined starting from
      # upper-left corner column by column by

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

      # where
      #          K-1                    L-1
      # R(K,L) = SUM [A(I,K)'*X(I,L)] + SUM [X(K,J)*A(J,L)].
      #          I=1                    J=1
      # """
      for l = 1:n
          for k = l:n
              y = C[k,l]
              for ia = l:k-1
                  y += A[ia,k]'*C[ia,l]
              end
              C[k,l] = -y/(A[k,k]'+A[l,l])
              isfinite(C[k,l]) || throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              k == l ? (C[k,k] = real(C[k,k])) : (C[l,k] = C[k,l]')
           end
           for jr = l+1:n
               for ir = jr:n
                   C[ir,jr] += C[ir,l]*A[l,jr] + (A[l,ir]*C[jr,l])'
                   if jr != ir
                      C[jr,ir] = C[ir,jr]'
                   end
               end
           end
      end
   else
      # """
      # The (K,L)th element of X is determined starting from
      # bottom-right corner column by column by

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

      # where
      #           N                     N
      # R(K,L) = SUM [A(K,I)*X(I,L)] + SUM [X(K,J)*A(L,J)'].
      #         I=K+1                 J=L+1
      # """
      for l = n:-1:1
          for k = l:-1:1
              y = C[k,l]
              for ia = k+1:l
                 y += A[k,ia]*C[ia,l]
              end
              C[k,l] = -y/(A[k,k]+A[l,l]')
              isfinite(C[k,l]) || throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              k == l ? (C[k,k] = real(C[k,k])) : (C[l,k] = C[k,l]')
           end
           for jr = 1:l-1
               for ir = 1:jr
                  C[ir,jr] += C[ir,l]*A[jr,l]' + A[ir,l]*C[l,jr]
                  if ir != jr
                     C[jr,ir] = C[ir,jr]'
                  end
               end
           end
       end
   end
end
"""
    lyapcs!(A, E, C; adj = false)

Solve the generalized continuous Lyapunov matrix equation

                op(A)Xop(E)' + op(E)Xop(A)' + C = 0,

where `op(A) = A` and `op(E) = E` if `adj = false` and `op(A) = A'` and
`op(E) = E'` if `adj = true`. The pair `(A,E)` is in a generalized real or
complex Schur form and `C` is a symmetric or hermitian matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapcs!(A::AbstractMatrix{T1},E::Union{AbstractMatrix{T1},UniformScaling{Bool}}, C::AbstractMatrix{T1}; adj = false) where {T1<:BlasReal}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && issymmetric(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapcs!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   ONE = one(T1)
   ZERO = zero(T1)

   # determine the structure of the real Schur form
   ba, p = sfstruct(A)

   G = Matrix{T1}(undef,2,2)
   W = Array{T1,2}(undef,n,2)
   Xw = Matrix{T1}(undef,4,4)
   Yw = Vector{T1}(undef,4)
   if adj
      # """
      # 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)*E(L,L) + E(K,K)'*X(K,L)*A(L,L) = -C(K,L) - R(K,L),

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

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

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

      #          K-1
      #         {SUM [E(I,K)'*X(I,L)]}*A(L,L).
      #          I=1
      # """
      i = 1
      for kk = 1:p
          dk = ba[kk]
          dkk = 1:dk
          k = i:i+dk-1
          j = 1
          ir = 1:i-1
          for ll = 1:kk
             dl = ba[ll]
             j1 = j+dl-1
             l = j:j1
             y = view(G,1:dk,1:dl)
             Ckl = view(C,k,l)
             copyto!(y,Ckl)
             if kk > 1
                 #   C[l,k] = C[l,ir]*A[ir,k]
                 #   W[l,dkk] = C[l,ir]*E[ir,k]
                 mul!(view(C,l,k),view(C,l,ir),view(A,ir,k))
                 mul!(view(W,l,dkk),view(C,l,ir),view(E,ir,k))
                 ic = 1:j1
                 #   y += C[ic,k]'*E[ic,l] + W[ic,dkk]'*A[ic,l]
                 mul!(y,transpose(view(C,ic,k)),view(E,ic,l),ONE,ONE)
                 mul!(y,transpose(view(W,ic,dkk)),view(A,ic,l),ONE,ONE)
             end
             if i == j
                lyapc2!(adj,y,dk,view(A,k,k),view(E,k,k),Xw,Yw)
                copyto!(Ckl,y)
               else
                lyapcsylv2!(adj,y,dk,dl,view(A,k,k),view(E,k,k),view(A,l,l),view(E,l,l),Xw,Yw)
                copyto!(Ckl,y)
             end
             j += dl
             if j <= i
                #  C[l,k] += C[k,l]'*A[k,k]
                #  W[l,dkk] += C[k,l]'*E[k,k]
                mul!(view(C,l,k),transpose(Ckl),view(A,k,k),ONE,ONE)
                mul!(view(W,l,dkk),transpose(Ckl),view(E,k,k),ONE,ONE)
             end
          end
          if kk > 1
             ir = 1:i-1
             #C[ir,k] = C[k,ir]'
             transpose!(view(C,ir,k),view(C,k,ir))
          end
          i += dk
      end
   else
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-right corner column by column by

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

      # where

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

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

      #             N
      #          { SUM [A(K,J)*X(J,L)]}*E(L,L)' +
      #           J=K+1

      #             N
      #          { SUM [E(K,J)*X(J,L)]}*A(L,L)'
      #           J=K+1
      # """
      j = n
      for ll = p:-1:1
        dl = ba[ll]
        l = j-dl+1:j
        dll = 1:dl
        i = n
        ir = j+1:n
        for kk = p:-1:ll
            dk = ba[kk]
            i1 = i-dk+1
            k = i1:i
            Clk = view(C,l,k)
            y = view(G,1:dl,1:dk)
            copyto!(y,Clk)
            if ll < p
               # C[k,l] = C[k,ir]*A[l,ir]'
               # W[k,dll] = C[k,ir]*E[l,ir]'
               mul!(view(C,k,l),view(C,k,ir),transpose(view(A,l,ir)))
               mul!(view(W,k,dll),view(C,k,ir),transpose(view(E,l,ir)))
               ic = i1:n
               # y += (E[k,ic]*C[ic,l] + A[k,ic]*W[ic,dll])'
               mul!(y,transpose(view(C,ic,l)),transpose(view(E,k,ic)),ONE,ONE)
               mul!(y,transpose(view(W,ic,dll)),transpose(view(A,k,ic)),ONE,ONE)
            end
            if i == j
               lyapc2!(adj,y,dk,view(A,k,k),view(E,k,k),Xw,Yw)
               copyto!(Clk,y)
            else
               lyapcsylv2!(adj,y,dl,dk,view(A,l,l),view(E,l,l),view(A,k,k),view(E,k,k),Xw,Yw)
               copyto!(Clk,y)
            end
            i -= dk
            if i >= j
               # C[k,l] += (A[l,l]*C[l,k])'
               # W[k,dll] += (E[l,l]*C[l,k])'
               mul!(view(C,k,l),transpose(Clk),transpose(view(A,l,l)),ONE,ONE)
               mul!(view(W,k,dll),transpose(Clk),transpose(view(E,l,l)),ONE,ONE)
            else
               break
            end
        end
        if ll < p
           ir = j+1:n
           # C[ir,l] = C[l,ir]'
           transpose!(view(C,ir,l),view(C,l,ir))
         end
        j -= dl
      end
   end
end
@inline function lyapc2!(adj,C::AbstractMatrix{T},na::Int,A::AbstractMatrix{T},E::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
# speed and reduced allocation oriented implementation of a solver for 1x1 or 2x2 generalized continuous Lyapunov equations
# LAPACK generated diagonal structure of E is exploited when possible
#      A'*X*E + E'*X*A = -C if adj = true  -> R = kron(E',A')+kron(A',E') = (kron(E,A)+kron(A,E))'
#      A*X*E' + E*X*A' = -C if adj = false -> R = kron(E,A)+kron(A,E)
   if na == 1
      temp = E[1,1]*A[1,1]
      rmul!(C,inv(-2*temp))
      any(!isfinite, C) &&  throw("ME:SingularException: A-λE has zero or infinite eigenvalue(s)")
      return C
   end
   ZERO = zero(T)
   i1 = 1:3
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   @inbounds  Y[1] = -C[1,1]/2
   @inbounds  Y[2] = -C[2,1]
   @inbounds  Y[3] = -C[2,2]/2
   if adj
      # Rt =
      #  [        a11*e11,                   a21*e11,         0
      #  a11*e12 + a12*e11, a11*e22 + a21*e12 + a22*e11,   a21*e22
      #          a12*e12,       a12*e22 + a22*e12, a22*e22]
      # iszero(E[1,2]) ?
      #    (@inbounds R = [  A[1,1]*E[1,1]             A[2,1]*E[1,1]               0 ;
      #                      A[1,2]*E[1,1]    A[1,1]*E[2,2]+A[2,2]*E[1,1]    A[2,1]*E[2,2];
      #                           0                   A[1,2]*E[2,2]          A[2,2]*E[2,2]]) :
      #    (@inbounds R = [  A[1,1]*E[1,1]                    A[2,1]*E[1,1]                              0 ;
      #          A[1,1]*E[1,2]+A[1,2]*E[1,1]    A[1,1]*E[2,2]+A[2,1]*E[1,2]+A[2,2]*E[1,1]    A[2,1]*E[2,2];
      #          A[1,2]*E[1,2]                    A[1,2]*E[2,2]+A[2,2]*E[1,2]                A[2,2]*E[2,2]] )
      @inbounds  R[1,1] = A[1,1]*E[1,1]
      @inbounds  R[1,2] = A[2,1]*E[1,1]
      @inbounds  R[1,3] = ZERO
      @inbounds  R[2,1] = A[1,1]*E[1,2]+A[1,2]*E[1,1]
      @inbounds  R[2,2] = A[1,1]*E[2,2]+A[2,1]*E[1,2]+A[2,2]*E[1,1]
      @inbounds  R[2,3] = A[2,1]*E[2,2]
      @inbounds  R[3,1] = A[1,2]*E[1,2]
      @inbounds  R[3,2] = A[1,2]*E[2,2]+A[2,2]*E[1,2]
      @inbounds  R[3,3] = A[2,2]*E[2,2]
   else
      # R =
      # [ a11*e11,       a11*e12 + a12*e11,         a12*e12
      # a21*e11, a11*e22 + a21*e12 + a22*e11, a12*e22 + a22*e12
      #       0,                   a21*e22,         a22*e22]
      # iszero(E[1,2]) ?
      #    (@inbounds  R = [  A[1,1]*E[1,1]            A[1,2]*E[1,1]                    0;
      #                       A[2,1]*E[1,1]    A[1,1]*E[2,2]+A[2,2]*E[1,1]     A[1,2]*E[2,2];
      #                             0                  A[2,1]*E[2,2]           A[2,2]*E[2,2]] ) :
      #    (@inbounds  R = [  A[1,1]*E[1,1]    A[1,1]*E[1,2]+A[1,2]*E[1,1]             A[1,2]*E[1,2];
      #           A[2,1]*E[1,1]    A[1,1]*E[2,2]+A[2,1]*E[1,2]+A[2,2]*E[1,1]     A[1,2]*E[2,2]+A[2,2]*E[1,2];
      #                 0                   A[2,1]*E[2,2]                              A[2,2]*E[2,2]] )
      @inbounds  R[1,1] = A[1,1]*E[1,1]
      @inbounds  R[1,2] = A[1,1]*E[1,2]+A[1,2]*E[1,1]
      @inbounds  R[1,3] = A[1,2]*E[1,2]
      @inbounds  R[2,1] = A[2,1]*E[1,1]
      @inbounds  R[2,2] = A[1,1]*E[2,2]+A[2,1]*E[1,2]+A[2,2]*E[1,1]
      @inbounds  R[2,3] = A[1,2]*E[2,2]+A[2,2]*E[1,2]
      @inbounds  R[3,1] = ZERO
      @inbounds  R[3,2] = A[2,1]*E[2,2]
      @inbounds  R[3,3] = A[2,2]*E[2,2]
   end
   luslv!(R,Y) && throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
   @inbounds C[1,1] = Y[1]
   @inbounds C[1,2] = Y[2]
   @inbounds C[2,1] = Y[2]
   @inbounds C[2,2] = Y[3]
   return C
end
@inline function lyapcsylv2!(adj,C::AbstractMatrix{T},na::Int,nb::Int,A::AbstractMatrix{T},E::AbstractMatrix{T},B::AbstractMatrix{T},F::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
# speed and reduced allocation oriented implementation of a solver for 1x1 and 2x2 Sylvester equations
# encountered in solving generalized continuous Lyapunov equations:
#      A'*X*F + E'*X*B = -C if adj = true  -> R = kron(F',A')+kron(B',E') = (kron(F,A)+kron(B,E))'
#      A*X*F' + E*X*B' = -C if adj = false -> R = kron(F,A)+kron(B,E)
   if na == 1 && nb == 1
      temp = A[1,1]*F[1,1] + E[1,1]*B[1,1]
      rmul!(C,inv(-temp))
      any(!isfinite, C) && throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
      return C
   end
   ZERO = zero(T)
   i1 = 1:na*nb
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   if adj
      if na == 1
         # R12t =
         # [ a11*f11+b11*e11           b21*e11]
         # [ a11*f12+b12*e11 a11*f22+b22*e11]
         # @inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1]           B[2,1]*E[1,1];
         #                 A[1,1]*F[1,2]+B[1,2]*E[1,1] A[1,1]*F[2,2]+B[2,2]*E[1,1]]
         @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
         @inbounds  R[1,2] = B[2,1]*E[1,1]
         @inbounds  R[2,1] = A[1,1]*F[1,2]+B[1,2]*E[1,1]
         @inbounds  R[2,2] = A[1,1]*F[2,2]+B[2,2]*E[1,1]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # R21t =
            # [ a11*f11+b11*e11           a21*f11]
            # [ a12*f11+b11*e12 a22*f11+b11*e22]
            # @inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1]           A[2,1]*F[1,1];
            #                 A[1,2]*F[1,1]+B[1,1]*E[1,2] A[2,2]*F[1,1]+B[1,1]*E[2,2] ]
            @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
            @inbounds  R[1,2] = A[2,1]*F[1,1]
            @inbounds  R[2,1] = A[1,2]*F[1,1]+B[1,1]*E[1,2]
            @inbounds  R[2,2] = A[2,2]*F[1,1]+B[1,1]*E[2,2]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # Rt =
            # [ a11*f11+b11*e11           a21*f11           b21*e11                 0]
            # [ a12*f11+b11*e12 a22*f11+b11*e22           b21*e12           b21*e22]
            # [ a11*f12+b12*e11           a21*f12 a11*f22+b22*e11           a21*f22]
            # [ a12*f12+b12*e12 a22*f12+b12*e22 a12*f22+b22*e12 a22*f22+b22*e22]
            # (iszero(E[1,2]) && iszero(F[1,2])) ?
            #  (@inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1]           A[2,1]*F[1,1]           B[2,1]*E[1,1]                 0;
            #                A[1,2]*F[1,1]  A[2,2]*F[1,1]+B[1,1]*E[2,2]           0           B[2,1]*E[2,2];
            #                B[1,2]*E[1,1]        0  A[1,1]*F[2,2]+B[2,2]*E[1,1]           A[2,1]*F[2,2];
            #                0  B[1,2]*E[2,2] A[1,2]*F[2,2]+B[2,2]*E[1,2] A[2,2]*F[2,2]+B[2,2]*E[2,2]]) :
            #  (@inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1]           A[2,1]*F[1,1]           B[2,1]*E[1,1]                 0;
            #                A[1,2]*F[1,1]+B[1,1]*E[1,2] A[2,2]*F[1,1]+B[1,1]*E[2,2]           B[2,1]*E[1,2]           B[2,1]*E[2,2];
            #                A[1,1]*F[1,2]+B[1,2]*E[1,1]           A[2,1]*F[1,2] A[1,1]*F[2,2]+B[2,2]*E[1,1]           A[2,1]*F[2,2];
            #                A[1,2]*F[1,2]+B[1,2]*E[1,2] A[2,2]*F[1,2]+B[1,2]*E[2,2] A[1,2]*F[2,2]+B[2,2]*E[1,2] A[2,2]*F[2,2]+B[2,2]*E[2,2]])
            @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
            @inbounds  R[1,2] = A[2,1]*F[1,1]
            @inbounds  R[1,3] = B[2,1]*E[1,1]
            @inbounds  R[1,4] = ZERO
            @inbounds  R[2,1] = A[1,2]*F[1,1]+B[1,1]*E[1,2]
            @inbounds  R[2,2] = A[2,2]*F[1,1]+B[1,1]*E[2,2]
            @inbounds  R[2,3] = B[2,1]*E[1,2]
            @inbounds  R[2,4] = B[2,1]*E[2,2]
            @inbounds  R[3,1] = A[1,1]*F[1,2]+B[1,2]*E[1,1]
            @inbounds  R[3,2] = A[2,1]*F[1,2]
            @inbounds  R[3,3] = A[1,1]*F[2,2]+B[2,2]*E[1,1]
            @inbounds  R[3,4] = A[2,1]*F[2,2]
            @inbounds  R[4,1] = A[1,2]*F[1,2]+B[1,2]*E[1,2]
            @inbounds  R[4,2] = A[2,2]*F[1,2]+B[1,2]*E[2,2]
            @inbounds  R[4,3] = A[1,2]*F[2,2]+B[2,2]*E[1,2]
            @inbounds  R[4,4] = A[2,2]*F[2,2]+B[2,2]*E[2,2]
            @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
   else
      if na == 1
         # R12 =
         # [ a11*f11+b11*e11 a11*f12+b12*e11]
         # [           b21*e11 a11*f22+b22*e11]
         # @inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1] A[1,1]*F[1,2]+B[1,2]*E[1,1];
         #                    B[2,1]*E[1,1] A[1,1]*F[2,2]+B[2,2]*E[1,1]]
         @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
         @inbounds  R[1,2] = A[1,1]*F[1,2]+B[1,2]*E[1,1]
         @inbounds  R[2,1] = B[2,1]*E[1,1]
         @inbounds  R[2,2] = A[1,1]*F[2,2]+B[2,2]*E[1,1]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*f11+b11*e11 a12*f11+b11*e12]
            # [           a21*f11 a22*f11+b11*e22]
            # @inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1] A[1,2]*F[1,1]+B[1,1]*E[1,2];
            #                  A[2,1]*F[1,1] A[2,2]*F[1,1]+B[1,1]*E[2,2]]
            @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
            @inbounds  R[1,2] = A[1,2]*F[1,1]+B[1,1]*E[1,2]
            @inbounds  R[2,1] = A[2,1]*F[1,1]
            @inbounds  R[2,2] = A[2,2]*F[1,1]+B[1,1]*E[2,2]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # R =
            # [ a11*f11+b11*e11 a12*f11+b11*e12 a11*f12+b12*e11 a12*f12+b12*e12]
            # [           a21*f11 a22*f11+b11*e22           a21*f12 a22*f12+b12*e22]
            # [           b21*e11           b21*e12 a11*f22+b22*e11 a12*f22+b22*e12]
            # [                 0           b21*e22           a21*f22 a22*f22+b22*e22]
            # (iszero(E[1,2]) && iszero(F[1,2])) ?
            # (@inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1] A[1,2]*F[1,1] B[1,2]*E[1,1] 0;
            #                A[2,1]*F[1,1] A[2,2]*F[1,1]+B[1,1]*E[2,2]           0   B[1,2]*E[2,2];
            #                B[2,1]*E[1,1]           0   A[1,1]*F[2,2]+B[2,2]*E[1,1] A[1,2]*F[2,2];
            #                0           B[2,1]*E[2,2]           A[2,1]*F[2,2] A[2,2]*F[2,2]+B[2,2]*E[2,2]]) :
            # (@inbounds R = [ A[1,1]*F[1,1]+B[1,1]*E[1,1] A[1,2]*F[1,1]+B[1,1]*E[1,2] A[1,1]*F[1,2]+B[1,2]*E[1,1] A[1,2]*F[1,2]+B[1,2]*E[1,2];
            #                A[2,1]*F[1,1] A[2,2]*F[1,1]+B[1,1]*E[2,2]           A[2,1]*F[1,2] A[2,2]*F[1,2]+B[1,2]*E[2,2];
            #                B[2,1]*E[1,1]           B[2,1]*E[1,2] A[1,1]*F[2,2]+B[2,2]*E[1,1] A[1,2]*F[2,2]+B[2,2]*E[1,2];
            #                0           B[2,1]*E[2,2]           A[2,1]*F[2,2] A[2,2]*F[2,2]+B[2,2]*E[2,2]])
            @inbounds  R[1,1] = A[1,1]*F[1,1]+B[1,1]*E[1,1]
            @inbounds  R[1,2] = A[1,2]*F[1,1]+B[1,1]*E[1,2]
            @inbounds  R[1,3] = A[1,1]*F[1,2]+B[1,2]*E[1,1]
            @inbounds  R[1,4] = A[1,2]*F[1,2]+B[1,2]*E[1,2]
            @inbounds  R[2,1] = A[2,1]*F[1,1]
            @inbounds  R[2,2] = A[2,2]*F[1,1]+B[1,1]*E[2,2]
            @inbounds  R[2,3] = A[2,1]*F[1,2]
            @inbounds  R[2,4] = A[2,2]*F[1,2]+B[1,2]*E[2,2]
            @inbounds  R[3,1] = B[2,1]*E[1,1]
            @inbounds  R[3,2] = B[2,1]*E[1,2]
            @inbounds  R[3,3] = A[1,1]*F[2,2]+B[2,2]*E[1,1]
            @inbounds  R[3,4] = A[1,2]*F[2,2]+B[2,2]*E[1,2]
            @inbounds  R[4,1] = ZERO
            @inbounds  R[4,2] = B[2,1]*E[2,2]
            @inbounds  R[4,3] = A[2,1]*F[2,2]
            @inbounds  R[4,4] = A[2,2]*F[2,2]+B[2,2]*E[2,2]
            @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
   end
   luslv!(R,Y) && throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
   C[:,:] = Y
   return C
end
function lyapcs!(A::AbstractMatrix{T1},E::Union{AbstractMatrix{T1},UniformScaling{Bool}}, C::AbstractMatrix{T1}; adj = false) where {T1<:BlasComplex}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapcs!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   W = Array{T1,1}(undef,n)
   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               W[l] = C[l,1]*E[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
                  W[l] += C[l,ir]*E[ir,k]
               end
               for ic = 1:l
                   y += C[ic,k]'*E[ic,l] + W[ic]'*A[ic,l]
                end
            end
            C[k,l] = -y/(A[k,k]'*E[l,l]+E[k,k]'*A[l,l])
            isfinite(C[k,l]) || throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
            k == l && (C[k,k] = real(C[k,k]))
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
               W[l] += C[k,l]'*E[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               W[k] = C[k,l+1]*E[l,l+1]'
               for ir = l+2:n
                   C[k,l] += C[k,ir]*A[l,ir]'
                   W[k] += C[k,ir]*E[l,ir]'
               end
               for ic = k:n
                   y += (E[k,ic]*C[ic,l] + A[k,ic]*W[ic])'
               end
            end
            C[l,k] = -y/(A[k,k]'*E[l,l]+E[k,k]'*A[l,l])
            isfinite(C[l,k]) ||  throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
            k == l && (C[k,k] = real(C[k,k]))
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
               W[k] += (E[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end

"""
    lyapds!(A, C; adj = false)

Solve the discrete Lyapunov matrix equation

                op(A)Xop(A)' - X + C = 0,

where `op(A) = A` if `adj = false` and `op(A) = A'` if `adj = true`.
`A` is in a real or complex Schur form and `C` is a symmetric or hermitian matrix.
`A` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapds!(A::AbstractMatrix{T1},C::AbstractMatrix{T1}; adj = false) where {T1<:BlasReal}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && issymmetric(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric matrix"))

   ONE = one(T1)

   # determine the dimensions of the diagonal blocks of real Schur form

   G = Matrix{T1}(undef,2,2)
   Xw = Matrix{T1}(undef,4,4)
   Yw = Vector{T1}(undef,4)
   ba, p = sfstruct(A)
   if adj
      # """
      # 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)*A(L,L) - X(K,L) = -C(K,L) - R(K,L),

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

      #           K-1
      #          {SUM [A(I,K)'*X(I,L)]}*A(L,L).
      #           I=1
      # """
      i = 1
      @inbounds  for kk = 1:p
          dk = ba[kk]
          k = i:i+dk-1
          j = 1
          ir = 1:i-1
          for ll = 1:kk
              dl = ba[ll]
              j1 = j+dl-1
              l = j:j1
              Ckl = view(C,k,l)
              y = view(G,1:dk,1:dl)
              copyto!(y,Ckl)
              if kk > 1
                # C[l,k] = C[l,ir]*A[ir,k]
                 mul!(view(C,l,k),view(C,l,ir),view(A,ir,k))
                 ic = 1:j1
                 #y += C[ic,k]'*A[ic,l]
                 mul!(y,transpose(view(C,ic,k)),view(A,ic,l),ONE,ONE)
              end
              if kk == ll
                 lyapd2!(adj,y,dk,view(A,k,k),Xw,Yw)
                 copyto!(Ckl,y)
              else
                 lyapdsylv2!(adj,y,dk,dl,view(A,k,k),view(A,l,l),Xw,Yw)
                 copyto!(Ckl,y)
              end
              j += dl
              if ll < kk
                 # C[l,k] += C[k,l]'*A[k,k]
                 mul!(view(C,l,k),transpose(Ckl),view(A,k,k),ONE,ONE)
               end
          end
          if kk > 1
             # C[ir,k] = C[k,ir]'
             transpose!(view(C,ir,k),view(C,k,ir))
          end
          i += dk
      end
   else
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-right corner column by column by

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

      # where

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

      #             N
      #          { SUM [A(K,J)*X(J,L)]}*A(L,L)'
      #           J=K+1
      # """
      j = n
      for ll = p:-1:1
          dl = ba[ll]
          l = j-dl+1:j
          i = n
          ir = j+1:n
          for kk = p:-1:ll
              dk = ba[kk]
              i1 = i-dk+1
              k = i1:i
              Clk = view(C,l,k)
              y = view(G,1:dl,1:dk)
              copyto!(y,Clk)
              if ll < p
                 # C[k,l] = C[k,ir]*A[l,ir]'
                 mul!(view(C,k,l),view(C,k,ir),transpose(view(A,l,ir)))
                 ic = i1:n
                 # y += (A[k,ic]*C[ic,l])'
                 mul!(y,transpose(view(C,ic,l)),transpose(view(A,k,ic)),ONE,ONE)
              end
              if i == j
                 lyapd2!(adj,y,dk,view(A,k,k),Xw,Yw)
                 copyto!(Clk,y)
              else
                 lyapdsylv2!(adj,y,dl,dk,view(A,l,l),view(A,k,k),Xw,Yw)
                 copyto!(Clk,y)
              end
              i -= dk
              if i >= j
                 #C[k,l] += (A[l,l]*C[l,k])'
                 mul!(view(C,k,l),transpose(Clk),transpose(view(A,l,l)),ONE,ONE)
              else
                 break
              end
          end
          if ll < p
             ir = i+2:n
             # C[ir,l] = C[l,ir]'
             transpose!(view(C,ir,l),view(C,l,ir))
          end
          j -= dl
      end
   end
end

function lyapd2!(adj,C::AbstractMatrix{T},na::Int,A::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
   # speed and reduced allocation oriented implementation of a solver for 1x1 or 2x2 continuous Lyapunov equations
   #      A'*X*A - X = -C if adj = true  -> R = kron(A',A')-I = (kron(A,A)-I)'
   #      A*X*A' - X = -C if adj = false -> R = kron(A,A)-I
   MONE = -one(T)
   if na == 1
      temp = A[1,1]^2+MONE
      rmul!(C,inv(-temp))
      any(!isfinite, C) && throw("ME:SingularException: A has eigenvalue(s) with moduli ≈ 1")
      return C
   end
   TWO = 2*one(T)
   i1 = 1:3
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   @inbounds  Y[1] = -C[1,1]
   @inbounds  Y[2] = -C[2,1]
   @inbounds  Y[3] = -C[2,2]
   if adj
      # Rt =
      # [ a11^2-1              2*a11*a21      a21^2]
      # [   a11*a12  a11*a22+a12*a21-1    a21*a22]
      # [     a12^2              2*a12*a22  a22^2-1]
      # @inbounds R = [ A[1,1]^2+MONE    TWO*A[1,1]*A[2,1]         A[2,1]^2;
      #                 A[1,1]*A[1,2]    A[1,1]*A[2,2]+A[1,2]*A[2,1]+MONE  A[2,1]*A[2,2];
      #                 A[1,2]^2         TWO*A[1,2]*A[2,2]         A[2,2]^2+MONE ]
      @inbounds  R[1,1] = A[1,1]^2+MONE
      @inbounds  R[1,2] = TWO*A[1,1]*A[2,1]
      @inbounds  R[1,3] = A[2,1]^2
      @inbounds  R[2,1] = A[1,1]*A[1,2]
      @inbounds  R[2,2] = A[1,1]*A[2,2]+A[1,2]*A[2,1]+MONE
      @inbounds  R[2,3] = A[2,1]*A[2,2]
      @inbounds  R[3,1] = A[1,2]^2
      @inbounds  R[3,2] = TWO*A[1,2]*A[2,2]
      @inbounds  R[3,3] = A[2,2]^2+MONE
   else
      # R =
      # [ a11^2-1              2*a11*a12      a12^2]
      # [   a11*a21  a11*a22+a12*a21-1    a12*a22]
      # [     a21^2              2*a21*a22  a22^2-1]

      # @inbounds R = [ A[1,1]^2+MONE    TWO*A[1,1]*A[1,2]         A[1,2]^2;
      #                 A[1,1]*A[2,1]    A[1,1]*A[2,2]+A[1,2]*A[2,1]+MONE  A[1,2]*A[2,2];
      #                 A[2,1]^2         TWO*A[2,1]*A[2,2]         A[2,2]^2+MONE ]
      @inbounds  R[1,1] = A[1,1]^2+MONE
      @inbounds  R[1,2] = TWO*A[1,1]*A[1,2]
      @inbounds  R[1,3] = A[1,2]^2
      @inbounds  R[2,1] = A[1,1]*A[2,1]
      @inbounds  R[2,2] = A[1,1]*A[2,2]+A[1,2]*A[2,1]+MONE
      @inbounds  R[2,3] = A[1,2]*A[2,2]
      @inbounds  R[3,1] = A[2,1]^2
      @inbounds  R[3,2] = TWO*A[2,1]*A[2,2]
      @inbounds  R[3,3] = A[2,2]^2+MONE
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
   @inbounds C[1,1] = Y[1]
   @inbounds C[1,2] = Y[2]
   @inbounds C[2,1] = Y[2]
   @inbounds C[2,2] = Y[3]
   return C
end
function lyapdsylv2!(adj::Bool,C::AbstractMatrix{T},na::Int,nb::Int,A::AbstractMatrix{T},B::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::AbstractVector{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 adj = true  -> R = kron(B',A') - I = (kron(B,A)-I)'
#    #      A*X*B' - X = -C if adj = false -> R = kron(B,A)-I

   MONE = -one(T)
   if na == 1 && nb == 1
      temp = A[1,1]*B[1,1] + MONE
      rmul!(C,inv(-temp))
      any(!isfinite, C) && throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
      return C
   end
   i1 = 1:na*nb
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   if adj
      # @inbounds  R = kron(transpose(B),transpose(A)) - I
      if na == 1
         # R12t =
         # [ a11*b11-1      a11*b21]
         # [     a11*b12  a11*b22-1]
         # @inbounds R = [ A[1,1]*B[1,1]+MONE      A[1,1]*B[2,1];
         #                 A[1,1]*B[1,2]  A[1,1]*B[2,2]+MONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
         @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]+MONE
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
   else
         if nb == 1
            # R21t =
            # [ a11*b11-1      a21*b11]
            # [     a12*b11  a22*b11-1]
            # @inbounds R = [ A[1,1]*B[1,1]+MONE      A[2,1]*B[1,1];
            #                 A[1,2]*B[1,1]  A[2,2]*B[1,1]+MONE ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
            @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]+MONE
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # Rt =
            # [ 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]+MONE      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]+MONE      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]+MONE      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]+MONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
            @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]+MONE
            @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]+MONE
            @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]+MONE
            @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
   else
      # @inbounds  R = kron(B,A) - I
      if na == 1
         # R12 =
         # [ a11*b11-1      a11*b12]
         # [     a11*b21  a11*b22-1]
         # @inbounds R = [ A[1,1]*B[1,1]+MONE      A[1,1]*B[1,2];
         #                 A[1,1]*B[2,1]  A[1,1]*B[2,2]+MONE]
         @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
         @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]+MONE
         @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]+MONE      A[1,2]*B[1,1];
            #                 A[2,1]*B[1,1]  A[2,2]*B[1,1]+MONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
            @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]+MONE
            @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]+MONE      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]+MONE      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]+MONE      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]+MONE]
            @inbounds  R[1,1] = A[1,1]*B[1,1]+MONE
            @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]+MONE
            @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]+MONE
            @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]+MONE
            @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
   end
   luslv!(R,Y) && throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
   C[:,:] = Y
   return C
end
function lyapds!(A::AbstractMatrix{T1},C::AbstractMatrix{T1}; adj = false) where {T1<:BlasComplex}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n hermitian matrix"))

   ONE = one(T1)

   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
               end
               for ic = 1:l
                  y += C[ic,k]'*A[ic,l]
               end
            end
            temp = ONE-A[k,k]'*A[l,l]
            C[k,l] = y/temp
            isfinite(C[k,l]) || throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
            k == l && (C[k,k] = real(C[k,k]))
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               for ir = l+2:n
                  C[k,l] += C[k,ir]*A[l,ir]'
               end
               for ic = k:n
                   y += (A[k,ic]*C[ic,l])'
               end
            end
            temp = ONE-A[k,k]'*A[l,l]
            C[l,k] = y/temp
            isfinite(C[l,k]) || throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
            k == l && (C[k,k] = real(C[k,k]))
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end
"""
    lyapds!(A, E, C; adj = false)

Solve the generalized discrete Lyapunov matrix equation

                op(A)Xop(A)' - op(E)Xop(E)' + C = 0,

where `op(A) = A` and `op(E) = E` if `adj = false` and `op(A) = A'` and
`op(E) = E'` if `adj = true`. The pair `(A,E)` is in a generalized real or
complex Schur form and `C` is a symmetric or hermitian matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapds!(A::AbstractMatrix{T1},E::Union{AbstractMatrix{T1},UniformScaling{Bool}}, C::AbstractMatrix{T1}; adj = false) where {T1<:BlasReal}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && issymmetric(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapds!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   ONE = one(T1)
   MONE = -ONE

   # determine the structure of the real Schur form
   ba, p = sfstruct(A)

   G = Matrix{T1}(undef,2,2)
   W = Array{T1,2}(undef,n,2)
   Xw = Matrix{T1}(undef,4,4)
   Yw = Vector{T1}(undef,4)
   if adj
      # """
      # 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)*A(L,L) - E(K,K)'*X(K,L)*E(L,L) = -C(K,L) - R(K,L),

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

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

      #           K-1
      #          {SUM [A(I,K)'*X(I,L)]}*A(L,L) -
      #           I=1

      #           K-1
      #          {SUM [E(I,K)'*X(I,L)]}*E(L,L).
      #           I=1
      # """
      i = 1
      for kk = 1:p
          dk = ba[kk]
          dkk = 1:dk
          k = i:i+dk-1
          j = 1
          ir = 1:i-1
          for ll = 1:kk
             dl = ba[ll]
             j1 = j+dl-1
             l = j:j1
             #y = C[k,l]
             y = view(G,1:dk,1:dl)
             Ckl = view(C,k,l)
             copyto!(y,Ckl)
             if kk > 1
                 # C[l,k] = C[l,ir]*A[ir,k]
                 # W[l,dkk] = C[l,ir]*E[ir,k]
                 mul!(view(C,l,k),view(C,l,ir),view(A,ir,k))
                 mul!(view(W,l,dkk),view(C,l,ir),view(E,ir,k))
                 ic = 1:j1
                 # y += C[ic,k]'*A[ic,l] -W[ic,dkk]'*E[ic,l]
                 mul!(y,transpose(view(C,ic,k)),view(A,ic,l),ONE,ONE)
                 mul!(y,transpose(view(W,ic,dkk)),view(E,ic,l),MONE,ONE)
             end
             if i == j
                lyapd2!(adj,y,dk,view(A,k,k),view(E,k,k),Xw,Yw)
                copyto!(Ckl,y)
             else
                lyapdsylv2!(adj,y,dk,dl,view(A,k,k),view(E,k,k),view(A,l,l),view(E,l,l),Xw,Yw)
                copyto!(Ckl,y)
             end
             j += dl
             if j <= i
                 # C[l,k] += C[k,l]'*A[k,k]
                 # W[l,dkk] += C[k,l]'*E[k,k]
                 mul!(view(C,l,k),transpose(view(C,k,l)),view(A,k,k),ONE,ONE)
                 mul!(view(W,l,dkk),transpose(view(C,k,l)),view(E,k,k),ONE,ONE)
              end
          end
          if kk > 1
             # C[ir,k] = C[k,ir]'
             transpose!(view(C,ir,k),view(C,k,ir))
          end
          i += dk
      end
   else
      # """
      # The (K,L)th block of X is determined starting from
      # bottom-right corner column by column by

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

      # where

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

      #           N            N
      #          SUM {E(K,I)* SUM [X(I,J)*E(L,J)']} +
      #          I=K         J=L+1

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

      #             N
      #          { SUM [E(K,J)*X(J,L)]}*E(L,L)'
      #           J=K+1
      # """
      j = n
      for ll = p:-1:1
        dl = ba[ll]
        l = j-dl+1:j
        dll = 1:dl
        i = n
        ir = j+1:n
        for kk = p:-1:ll
            dk = ba[kk]
            i1 = i-dk+1
            k = i1:i
            Clk = view(C,l,k)
            y = view(G,1:dl,1:dk)
            copyto!(y,Clk)
            if ll < p
               # C[k,l] = C[k,ir]*A[l,ir]'
               # W[k,dll] = C[k,ir]*E[l,ir]'
               mul!(view(C,k,l),view(C,k,ir),transpose(view(A,l,ir)))
               mul!(view(W,k,dll),view(C,k,ir),transpose(view(E,l,ir)))
               ic = i1:n
               # y += (A[k,ic]*C[ic,l] - E[k,ic]*W[ic,dll])'
               mul!(y,transpose(view(C,ic,l)),transpose(view(A,k,ic)),ONE,ONE)
               mul!(y,transpose(view(W,ic,dll)),transpose(view(E,k,ic)),MONE,ONE)
            end
            if i == j
               lyapd2!(adj,y,dk,view(A,k,k),view(E,k,k),Xw,Yw)
               copyto!(Clk,y)
            else
               lyapdsylv2!(adj,y,dl,dk,view(A,l,l),view(E,l,l),view(A,k,k),view(E,k,k),Xw,Yw)
               copyto!(Clk,y)
            end
            i -= dk
            if i >= j
               # C[k,l] += (A[l,l]*C[l,k])'
               # W[k,dll] += (E[l,l]*C[l,k])'
               mul!(view(C,k,l),transpose(view(C,l,k)),transpose(view(A,l,l)),ONE,ONE)
               mul!(view(W,k,dll),transpose(view(C,l,k)),transpose(view(E,l,l)),ONE,ONE)
            else
               break
            end
        end
        if ll < p
           ir = i+2:n
           # C[ir,l] = C[l,ir]'
           transpose!(view(C,ir,l),view(C,l,ir))
         end
        j -= dl
      end
   end
end
@inline function lyapd2!(adj,C::AbstractMatrix{T},na::Int,A::AbstractMatrix{T},E::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{T}) where T <:BlasReal
   # speed and reduced allocation oriented implementation of a solver for 1x1 or 2x2 continuous Lyapunov equations
   #      A'*X*A - E'*X*E = -C if adj = true  -> R = kron(A',A')-kron(E',E') = (kron(A,A)-kron(E,E))'
   #      A*X*A' - E*X*E' = -C if adj = false -> R = kron(A,A)-kron(E,E)
   MONE = -one(T)
   if na == 1
      temp = E[1,1]^2-A[1,1]^2
      rmul!(C,inv(temp))
      any(!isfinite, C) &&  throw("ME:SingularException: A-λE has eigenvalue(s) with moduli equal to one")
      return C
   end
   TWO = 2*one(T)
   i1 = 1:3
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   @inbounds  Y[1] = -C[1,1]
   @inbounds  Y[2] = -C[2,1]
   @inbounds  Y[3] = -C[2,2]
   if adj
      # Rt =
      # [     a11^2 - e11^2,                   2*a11*a21,         a21^2]
      # [ a11*a12 - e11*e12, a11*a22 + a12*a21 - e11*e22,       a21*a22]
      # [     a12^2 - e12^2,       2*a12*a22 - 2*e12*e22, a22^2 - e22^2]
      # iszero(E[1,2]) ?
      #    (@inbounds R = [ A[1,1]^2-E[1,1]^2    TWO*A[1,1]*A[2,1]         A[2,1]^2;
      #                 A[1,1]*A[1,2]    A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]  A[2,1]*A[2,2];
      #                 A[1,2]^2         TWO*A[1,2]*A[2,2]         A[2,2]^2-E[2,2]^2 ] ) :
         # (@inbounds R = [ A[1,1]^2-E[1,1]^2    TWO*A[1,1]*A[2,1]         A[2,1]^2;
         #              A[1,1]*A[1,2]-E[1,1]*E[1,2]    A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]  A[2,1]*A[2,2];
         #              A[1,2]^2-E[1,2]^2         TWO*(A[1,2]*A[2,2]-E[1,2]*E[2,2])         A[2,2]^2-E[2,2]^2 ])
      @inbounds  R[1,1] = A[1,1]^2-E[1,1]^2
      @inbounds  R[1,2] = TWO*A[1,1]*A[2,1]
      @inbounds  R[1,3] = A[2,1]^2
      @inbounds  R[2,1] = A[1,1]*A[1,2]-E[1,1]*E[1,2]
      @inbounds  R[2,2] = A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]
      @inbounds  R[2,3] = A[2,1]*A[2,2]
      @inbounds  R[3,1] = A[1,2]^2-E[1,2]^2
      @inbounds  R[3,2] = TWO*(A[1,2]*A[2,2]-E[1,2]*E[2,2])
      @inbounds  R[3,3] = A[2,2]^2-E[2,2]^2
   else
      # R =
      # [ a11^2 - e11^2,       2*a11*a12 - 2*e11*e12,     a12^2 - e12^2]
      # [       a11*a21, a11*a22 + a12*a21 - e11*e22, a12*a22 - e12*e22]
      # [         a21^2,                   2*a21*a22,     a22^2 - e22^2]
      # iszero(E[1,2]) ?
      #    (@inbounds R = [ A[1,1]^2-E[1,1]^2    TWO*A[1,1]*A[1,2]         A[1,2]^2;
      #                 A[1,1]*A[2,1]    A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]  A[1,2]*A[2,2];
      #                 A[2,1]^2         TWO*A[2,1]*A[2,2]         A[2,2]^2-E[2,2]^2 ] ) :

      #    (@inbounds R = [ A[1,1]^2-E[1,1]^2    TWO*(A[1,1]*A[1,2]-E[1,1]*E[1,2])         A[1,2]^2-E[1,2]^2;
      #                 A[1,1]*A[2,1]    A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]  A[1,2]*A[2,2]-E[1,2]*E[2,2];
      #                 A[2,1]^2         TWO*A[2,1]*A[2,2]         A[2,2]^2-E[2,2]^2 ] )
      @inbounds  R[1,1] = A[1,1]^2-E[1,1]^2
      @inbounds  R[1,2] = TWO*(A[1,1]*A[1,2]-E[1,1]*E[1,2])
      @inbounds  R[1,3] = A[1,2]^2-E[1,2]^2
      @inbounds  R[2,1] = A[1,1]*A[2,1]
      @inbounds  R[2,2] = A[1,1]*A[2,2]+A[1,2]*A[2,1]-E[1,1]*E[2,2]
      @inbounds  R[2,3] = A[1,2]*A[2,2]-E[1,2]*E[2,2]
      @inbounds  R[3,1] = A[2,1]^2
      @inbounds  R[3,2] = TWO*A[2,1]*A[2,2]
      @inbounds  R[3,3] = A[2,2]^2-E[2,2]^2
   end
   luslv!(R,Y) && throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
   @inbounds C[1,1] = Y[1]
   @inbounds C[1,2] = Y[2]
   @inbounds C[2,1] = Y[2]
   @inbounds C[2,2] = Y[3]
   return C
end
@inline function lyapdsylv2!(adj,C::AbstractMatrix{T},na::Int,nb::Int,A::AbstractMatrix{T},E::AbstractMatrix{T},B::AbstractMatrix{T},F::AbstractMatrix{T},Xw::AbstractMatrix{T},Yw::StridedVector{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 - E'*X*F = -C if adj = true  -> R = kron(B',A') - kron(F',E') = (kron(B,A)-kron(F,E))'
   #      A*X*B' - E*X*F' = -C if adj = false -> R = kron(B,A)-kron(F,E)
   if na == 1 && nb == 1
      temp = E[1,1]*F[1,1] - A[1,1]*B[1,1]
      rmul!(C,inv(temp))
      any(!isfinite, C) && throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
      return C
   end
   i1 = 1:na*nb
   R = view(Xw,i1,i1)
   Y = view(Yw,i1)
   if adj
      if na == 1
         # R12t =
         # [ a11*b11 - e11*f11,           a11*b21]
         # [ a11*b12 - e11*f12, a11*b22 - e11*f22]
         # @inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[1,1]      A[1,1]*B[2,1];
         #                 A[1,1]*B[1,2]-E[1,1]*F[1,2]  A[1,1]*B[2,2]-E[1,1]*F[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
         @inbounds  R[1,2] = A[1,1]*B[2,1]
         @inbounds  R[2,1] = A[1,1]*B[1,2]-E[1,1]*F[1,2]
         @inbounds  R[2,2] = A[1,1]*B[2,2]-E[1,1]*F[2,2]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # R21t =
            # [ a11*b11 - e11*f11,           a21*b11]
            # [ a12*b11 - e12*f11, a22*b11 - e22*f11]
            # @inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[1,1]      A[2,1]*B[1,1];
            #                 A[1,2]*B[1,1]-E[1,2]*F[1,1]  A[2,2]*B[1,1]-E[2,2]*F[1,1] ]
            @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
            @inbounds  R[1,2] = A[2,1]*B[1,1]
            @inbounds  R[2,1] = A[1,2]*B[1,1]-E[1,2]*F[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]-E[2,2]*F[1,1]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # Rt =
            # [ a11*b11 - e11*f11,           a21*b11,           a11*b21,           a21*b21]
            # [ a12*b11 - e12*f11, a22*b11 - e22*f11,           a12*b21,           a22*b21]
            # [ a11*b12 - e11*f12,           a21*b12, a11*b22 - e11*f22,           a21*b22]
            # [ a12*b12 - e12*f12, a22*b12 - e22*f12, a12*b22 - e12*f22, a22*b22 - e22*f22]
            # (iszero(E[1,2]) && iszero(F[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[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]-E[2,2]*F[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]-E[1,1]*F[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]-E[2,2]*F[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[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]-E[1,1]*F[1,2]  A[2,2]*B[1,1]-E[2,2]*F[1,1]      A[1,2]*B[2,1]      A[2,2]*B[2,1];
            # A[1,1]*B[1,2]-E[1,1]*F[1,2]      A[2,1]*B[1,2]  A[1,1]*B[2,2]-E[1,1]*F[2,2]      A[2,1]*B[2,2];
            # A[1,2]*B[1,2]-E[1,2]*F[1,2]      A[2,2]*B[1,2]-E[2,2]*F[1,2]      A[1,2]*B[2,2]-E[1,2]*F[2,2]  A[2,2]*B[2,2]-E[2,2]*F[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
            @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]-E[1,1]*F[1,2]
            @inbounds  R[2,2] = A[2,2]*B[1,1]-E[2,2]*F[1,1]
            @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]-E[1,1]*F[1,2]
            @inbounds  R[3,2] = A[2,1]*B[1,2]
            @inbounds  R[3,3] = A[1,1]*B[2,2]-E[1,1]*F[2,2]
            @inbounds  R[3,4] = A[2,1]*B[2,2]
            @inbounds  R[4,1] = A[1,2]*B[1,2]-E[1,2]*F[1,2]
            @inbounds  R[4,2] = A[2,2]*B[1,2]-E[2,2]*F[1,2]
            @inbounds  R[4,3] = A[1,2]*B[2,2]-E[1,2]*F[2,2]
            @inbounds  R[4,4] = A[2,2]*B[2,2]-E[2,2]*F[2,2]
            @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
   else
      if na == 1
         # R12 =
         # [ a11*b11 - e11*f11, a11*b12 - e11*f12]
         # [           a11*b21, a11*b22 - e11*f22]
         # @inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[1,1]      A[1,1]*B[1,2]-E[1,1]*F[1,2];
         #                 A[1,1]*B[2,1]  A[1,1]*B[2,2]-E[1,1]*F[2,2]]
         @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
         @inbounds  R[1,2] = A[1,1]*B[1,2]-E[1,1]*F[1,2]
         @inbounds  R[2,1] = A[1,1]*B[2,1]
         @inbounds  R[2,2] = A[1,1]*B[2,2]-E[1,1]*F[2,2]
         @inbounds  Y[1] = -C[1,1]
         @inbounds  Y[2] = -C[1,2]
      else
         if nb == 1
            # R21 =
            # [ a11*b11 - e11*f11, a12*b11 - e12*f11]
            # [           a21*b11, a22*b11 - e22*f11]
            # @inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[1,1]      A[1,2]*B[1,1]-E[1,2]*F[1,1];
            #                 A[2,1]*B[1,1]  A[2,2]*B[1,1]-E[2,2]*F[1,1]]
            @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]-E[1,2]*F[1,1]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]-E[2,2]*F[1,1]
            @inbounds  Y[1] = -C[1,1]
            @inbounds  Y[2] = -C[2,1]
         else
            # R =
            # [ a11*b11 - e11*f11, a12*b11 - e12*f11, a11*b12 - e11*f12, a12*b12 - e12*f12]
            # [           a21*b11, a22*b11 - e22*f11,           a21*b12, a22*b12 - e22*f12]
            # [           a11*b21,           a12*b21, a11*b22 - e11*f22, a12*b22 - e12*f22]
            # [           a21*b21,           a22*b21,           a21*b22, a22*b22 - e22*f22]
            # (iszero(E[1,2]) && iszero(F[1,2])) ?
            # (@inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[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]-E[2,2]*F[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]-E[1,1]*F[2,2]      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]-E[2,2]*F[2,2]]) :
            # (@inbounds R = [ A[1,1]*B[1,1]-E[1,1]*F[1,1]      A[1,2]*B[1,1]-E[1,2]*F[1,1]      A[1,1]*B[1,2]-E[1,1]*F[1,2]      A[1,2]*B[1,2]-E[1,2]*F[1,2];
            # A[2,1]*B[1,1]  A[2,2]*B[1,1]-E[2,2]*F[1,1]      A[2,1]*B[1,2]      A[2,2]*B[1,2]-E[2,2]*F[1,2];
            # A[1,1]*B[2,1]      A[1,2]*B[2,1]  A[1,1]*B[2,2]-E[1,1]*F[2,2]      A[1,2]*B[2,2]-E[1,2]*F[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]-E[2,2]*F[2,2]])
            @inbounds  R[1,1] = A[1,1]*B[1,1]-E[1,1]*F[1,1]
            @inbounds  R[1,2] = A[1,2]*B[1,1]-E[1,2]*F[1,1]
            @inbounds  R[1,3] = A[1,1]*B[1,2]-E[1,1]*F[1,2]
            @inbounds  R[1,4] = A[1,2]*B[1,2]-E[1,2]*F[1,2]
            @inbounds  R[2,1] = A[2,1]*B[1,1]
            @inbounds  R[2,2] = A[2,2]*B[1,1]-E[2,2]*F[1,1]
            @inbounds  R[2,3] = A[2,1]*B[1,2]
            @inbounds  R[2,4] = A[2,2]*B[1,2]-E[2,2]*F[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]-E[1,1]*F[2,2]
            @inbounds  R[3,4] = A[1,2]*B[2,2]-E[1,2]*F[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]-E[2,2]*F[2,2]
            @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
   end
   luslv!(R,Y) && throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
   C[:,:] = Y
   return C
end
function lyapds!(A::AbstractMatrix{T1},E::Union{AbstractMatrix{T1},UniformScaling{Bool}}, C::AbstractMatrix{T1}; adj = false) where {T1<:BlasComplex}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapds!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   W = Array{T1,1}(undef,n)
   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               W[l] = C[l,1]*E[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
                  W[l] += C[l,ir]*E[ir,k]
               end
               for ic = 1:l
                   y += C[ic,k]'*A[ic,l] - W[ic]'*E[ic,l]
                end
            end
            C[k,l] = y/(E[k,k]'*E[l,l]-A[k,k]'*A[l,l])
            isfinite(C[k,l]) || throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
            k == l && (C[k,k] = real(C[k,k]))
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
               W[l] += C[k,l]'*E[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               W[k] = C[k,l+1]*E[l,l+1]'
               for ir = l+2:n
                  C[k,l] += C[k,ir]*A[l,ir]'
                  W[k] += C[k,ir]*E[l,ir]'
               end
               for ic = k:n
                   y += (A[k,ic]*C[ic,l] - E[k,ic]*W[ic])'
               end
            end
            C[l,k] = y/(E[k,k]'*E[l,l]-A[k,k]'*A[l,l])
            isfinite(C[l,k]) || throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
            k == l && (C[k,k] = real(C[k,k]))
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
               W[k] += (E[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
struct IndicatorBridge{T, B<:MOIBC.SetMapBridge{T}, A, IS_MOI_OR_CS, F, S} <: MOIBC.AbstractBridge
    con_idx::CI
end

function MOI.supports_constraint(
    ::Type{<:IndicatorBridge{T, B}},
    ::Type{F},
    ::Type{MOI.IndicatorSet{A,S}}
) where {T, B, F<:MOI.VectorAffineFunction, A, S}
    is_supported = MOI.supports_constraint(B, MOIU.scalar_type(F), S)
    !is_supported && return false
    S <: AbstractBoolSet && return true

    concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    return supports_concreteB(concrete_B)
end

function MOI.supports_constraint(
    ::Type{<:IndicatorBridge{T, B}},
    ::Type{F},
    ::Type{CS.IndicatorSet{A,IF,S}}
) where {T, B, F<:MOI.VectorAffineFunction, A, IF, S}
    is_supported = MOI.supports_constraint(B, IF, S)
    !is_supported && return false
    S <: AbstractBoolSet && return true

    concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    return supports_concreteB(concrete_B)
end

function MOIBC.concrete_bridge_type(
    IB::Type{<:IndicatorBridge{T,B}},
    G::Type{<:MOI.VectorAffineFunction},
    ::Type{IS},
) where {T,B,A,S,IS<:MOI.IndicatorSet{A,S}}
    concrete_B = get_concrete_B(IB,S)
    return IndicatorBridge{T,concrete_B,A,MOI.IndicatorSet,MOI.ScalarAffineFunction,S}
end

function MOIBC.concrete_bridge_type(
    IB::Type{<:IndicatorBridge{T,B}},
    G::Type{<:MOI.VectorAffineFunction},
    ::Type{IS},
) where {T,B,A,S,IF,IS<:CS.IndicatorSet{A,IF,S}}
    concrete_B = get_concrete_B(IB,S)
    return IndicatorBridge{T,concrete_B,A,CS.IndicatorSet,IF,S}
end

function get_concrete_B(
    ::Type{<:IndicatorBridge{T,B}},
    S
) where {T,B}
    if S <: AbstractBoolSet
        return B
    else
        return MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    end
end

function MOIB.added_constraint_types(
    ::Type{<:IndicatorBridge{T,B,A,IS_MOI_OR_CS,F,S}}
) where {T,B,A,IS_MOI_OR_CS<:MOI.IndicatorSet,F,S}
    added_constraints = added_constraint_types(B, S)
    return [(MOI.VectorAffineFunction{T}, MOI.IndicatorSet{A,added_constraints[1][2]})]
end

function MOIB.added_constraint_types(
    ::Type{<:IndicatorBridge{T,B,A,IS_MOI_OR_CS,F,S}}
) where {T,B,A,IS_MOI_OR_CS<:CS.IndicatorSet,F,S}
    added_constraints = added_constraint_types(B, S)
    return [(MOI.VectorAffineFunction{T}, CS.IndicatorSet{A,F,added_constraints[1][2]})]
end

function MOIB.added_constrained_variable_types(::Type{<:IndicatorBridge{T,B}}) where {T,B}
    return MOIB.added_constrained_variable_types(B)
end

function MOIBC.bridge_constraint(::Type{<:IndicatorBridge{T, B, A, IS_MOI_OR_CS, F, S}}, model, func, set) where {T, B, A, IS_MOI_OR_CS<:MOI.IndicatorSet, F, S}
    f = MOIU.eachscalar(func)
    new_func = MOIU.operate(vcat, T, f[1], map_function(B, f[2:end], set.set))
    new_inner_set = MOIBC.map_set(B, set.set)
    new_set = MOI.IndicatorSet{A}(new_inner_set)
    return IndicatorBridge{T,B,A,IS_MOI_OR_CS,F,S}(MOI.add_constraint(model, new_func, new_set))
end

function MOIBC.bridge_constraint(::Type{<:IndicatorBridge{T, B, A, IS_MOI_OR_CS, F, S}}, model, func, set) where {T, B, A, IS_MOI_OR_CS<:CS.IndicatorSet, F, S}
    f = MOIU.eachscalar(func)
    new_func = MOIU.operate(vcat, T, f[1], map_function(B, f[2:end], set.set))
    new_inner_set = MOIBC.map_set(B, set.set)
    new_set = CS.IndicatorSet{A,F}(new_inner_set)
    return IndicatorBridge{T,B,A,IS_MOI_OR_CS,F,S}(MOI.add_constraint(model, new_func, new_set))
end


text:
import KernelAbstractions: CPU

using MPI
using Printf

using ..BalanceLaws
using ..MPIStateArrays
using ..DGMethods: SpaceDiscretization
using ..VariableTemplates

"""
    VTKFieldWriter(
        name::String,
        FT::DataType,
        fields::Vector{<:Tuple{String, <:Function}};
        path_prefix = ".",
        number_sample_points = 0,
    )

Construct a callable type that computes the specified `fields` and
writes them to a VTK file in `path_prefix`/. `number_sample_points`
are passed to [`writevtk`](@ref).

Intended for use in a callback passed to a running simulation; files
are suffixed with a running number beginning with 0.

# Example
```julia
    function pres_fun(atmos::AtmosModel, prognostic::Vars, auxiliary::Vars)
        ts = recover_thermo_state(atmos, prognostic, auxiliary)
        air_pressure(ts)
    end
    fwriter = VTKFieldWriter(solver_config.name, [("pressure", pres_fun)])
    cbfw = GenericCallbacks.EveryXSimulationTime(60) do
        fwriter(solver_config.dg, solver_config.Q)
    end
```
"""
mutable struct VTKFieldWriter
    path_prefix::String
    name::String
    number_sample_points::Int
    nfields::Int
    field_names::Vector{String}
    field_funs::Vector{<:Function}
    vars_type::DataType
    num::Int

    function VTKFieldWriter(
        name::String,
        FT::DataType,
        fields::Vector{<:Tuple{String, <:Function}};
        path_prefix = ".",
        number_sample_points = 0,
    )
        nfields = length(fields)
        field_names = [name for (name, _) in fields]
        field_funs = [fun for (_, fun) in fields]
        vars_type = NamedTuple{
            tuple(Symbol.(field_names)...),
            Tuple{[FT for _ in 1:length(fields)]...},
        }
        new(
            path_prefix,
            name * "_fields",
            number_sample_points,
            nfields,
            field_names,
            field_funs,
            vars_type,
            0,
        )
    end
end
function (vfw::VTKFieldWriter)(dg::SpaceDiscretization, Q::MPIStateArray)
    bl = dg.balance_law
    fQ = similar(Q, Array; vars = vfw.vars_type, nstate = vfw.nfields)

    if array_device(Q) isa CPU
        prognostic_array = Q.realdata
        auxiliary_array = dg.state_auxiliary.realdata
    else
        prognostic_array = Array(Q.realdata)
        auxiliary_array = Array(dg.state_auxiliary.realdata)
    end
    FT = eltype(prognostic_array)

    for e in 1:size(prognostic_array, 3)
        for n in 1:size(prognostic_array, 1)
            prognostic = Vars{vars_state(bl, Prognostic(), FT)}(view(
                prognostic_array,
                n,
                :,
                e,
            ),)
            auxiliary = Vars{vars_state(bl, Auxiliary(), FT)}(view(
                auxiliary_array,
                n,
                :,
                e,
            ),)
            for i in 1:(vfw.nfields)
                fQ[n, i, e] = vfw.field_funs[i](bl, prognostic, auxiliary)
            end
        end
    end

    mpirank = MPI.Comm_rank(fQ.mpicomm)
    vprefix = @sprintf("%s_mpirank%04d_num%04d", vfw.name, mpirank, vfw.num,)
    outprefix = joinpath(vfw.path_prefix, vprefix)
    writevtk(outprefix, fQ, dg, number_sample_points = vfw.number_sample_points)

    # Generate the pvtu file for these vtk files
    if mpirank == 0
        pprefix = @sprintf("%s_num%04d", vfw.name, vfw.num)
        pvtuprefix = joinpath(vfw.path_prefix, pprefix)

        prefixes = ntuple(MPI.Comm_size(fQ.mpicomm)) do i
            @sprintf("%s_mpirank%04d_num%04d", vfw.name, i - 1, vfw.num,)
        end
        writepvtu(pvtuprefix, prefixes, tuple(vfw.field_names...), eltype(fQ))
    end

    vfw.num += 1
    nothing
end