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

Solve the continuous Sylvester matrix equation

                AX + XB = C

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

The following particular cases are also adressed:

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

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

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

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

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

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

    x = sylvc(α,β,γ)

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

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

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

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

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

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

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

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

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

   realcase = eltype(A) <: AbstractFloat
   if adjA
      RA, QA = schur(A.parent)
   else
      RA, QA = schur(A)
   end
   if adjB
      RB, QB = schur(B.parent)
   else
      RB, QB = schur(B)
   end
   D = adjoint(QA) * (C*QB)
   Y = sylvcs!(RA, RB, D, adjA = adjA, adjB = adjB)
   QA*(Y * adjoint(QB))
end
# solve X(B+α) = C or (α+β)X = C
sylvc(A::Union{Real,Complex,UniformScaling},B::Union{AbstractMatrix,UniformScaling},C::AbstractMatrix) = C/(A*I+B)
# solve (A+β)X = C
sylvc(A::AbstractMatrix,B::Union{Real,Complex,UniformScaling},C::AbstractMatrix) = (A+B*I)\C
# solve (α+β)X = C
sylvc(A::Union{Real,Complex},B::Union{Real,Complex},C::AbstractMatrix) = A+B == 0 ? throw(SingularException(1)) : C/(A+B)
# solve (α+β)x = γ
sylvc(A::Union{Real,Complex}, B::Union{Real,Complex}, C::Union{Real,Complex}) = A+B == 0 ? throw(SingularException(1)) : C/(A+B)
"""
    X = sylvd(A,B,C)

Solve the discrete Sylvester matrix equation

                AXB + X = C

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

The following particular cases are also adressed:

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

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

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

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

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

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

    x = sylvd(α,β,γ)

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

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

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

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

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

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

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

   adjA = isa(A,Adjoint)
   adjB = isa(B,Adjoint)
   if adjA
      RA, QA = schur(A.parent)
   else
      RA, QA = schur(A)
   end
   if adjB
      RB, QB = schur(B.parent)
   else
      RB, QB = schur(B)
   end
   D = adjoint(QA) * (C*QB)
   Y = sylvds!(RA, RB, D, adjA = adjA, adjB = adjB)
   QA*(Y * adjoint(QB))
end
# solve X(αB+I) = C or X(αβ+1) = C
sylvd(A::Union{Real,Complex,UniformScaling},B::Union{AbstractMatrix,UniformScaling},C::AbstractMatrix) = C/(A*B+I)
# solve (Aβ+I)X = C
sylvd(A::AbstractMatrix,B::Union{Real,Complex,UniformScaling},C::AbstractMatrix) = (A*B+I)\C
# solve (αβ+1)X = C
sylvd(A::Union{Real,Complex},B::Union{Real,Complex},C::AbstractMatrix) = A*B+1 == 0 ? throw(SingularException(1)) : C/(A*B+1)
# solve (αβ+1)x = γ
sylvd(A::Union{Real,Complex}, B::Union{Real,Complex}, C::Union{Real,Complex}) = A*B+1 == 0 ? throw(SingularException(1)) : C/(A*B+one(C))
"""
    X = gsylv(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

              AXB + CXD = E

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

The following particular cases are also adressed:

    X = gsylv(A,B,E)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if adjAC
       AS, CS, Z1, Q1 = schur(A.parent,C.parent)
    else
       if adjA
          A = copy(A)
       end
       if adjC
          C = copy(C)
       end
       AS, CS, Q1, Z1 = schur(A,C)
    end
    if adjBD
       BS, DS, Z2, Q2 = schur(B.parent,D.parent)
    else
      if adjB
          B = copy(B)
      end
      if adjD
          D = copy(D)
      end
       BS, DS, Q2, Z2 = schur(B,D)
    end
    Y = adjoint(Q1) * (E*Z2)
    gsylvs!(AS, BS, CS, DS, Y, adjAC = adjAC, adjBD = adjBD)
    Z1*(Y * adjoint(Q2))
end
# solve AXB = C
gsylv(A::Union{AbstractMatrix,UniformScaling,Real,Complex},B::Union{AbstractMatrix,UniformScaling,Real,Complex},E::AbstractMatrix) = (A\E)/B
# solve AXB+γX = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},E::AbstractMatrix) =
size(A,1) == size(A,2) && size(B,1) == size(B,2) ? gsylv(A,B,Matrix{eltype(C)}(C*I,size(A)),Matrix{eltype(C)}(I,size(B)),E) :
throw(DimensionMismatch("A and B must be square matrices"))
# solve AXB+γXδ = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},D::Union{UniformScaling,Real,Complex},E::AbstractMatrix) = gsylv(A,B,C*D,E)
# solve AXB+γXD = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::Union{UniformScaling,Real,Complex},D::AbstractMatrix,E::AbstractMatrix) =
size(A,1) == size(A,2) ? gsylv(A,B,Matrix{eltype(C)}(C*I,size(A)),D,E) :
throw(DimensionMismatch("A must be a square matrix"))
# solve AXB+CXδ = E
gsylv(A::AbstractMatrix,B::AbstractMatrix,C::AbstractMatrix,D::Union{UniformScaling,Real,Complex},E::AbstractMatrix) =
size(B,1) == size(B,2) ? gsylv(A,B,C,Matrix{eltype(D)}(D*I,size(B)),E) :
throw(DimensionMismatch("B must be a square matrix"))

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

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F,

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

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

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

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

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

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

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

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

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

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

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

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

    if isa(A,Adjoint)
      A = copy(A)
    end
    if isa(B,Adjoint)
      B = copy(B)
    end
    if isa(D,Adjoint)
      D = copy(D)
    end
    if isa(E,Adjoint)
      E = copy(E)
    end
    AS, DS, Q1, Z1 = schur(A,D)
    BS, ES, Q2, Z2 = schur(B,E)

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

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

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

Solve the dual Sylvester system of matrix equations

       AX + DY = C
       XB + YE = F ,

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Solve the continuous Sylvester matrix equation

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

where `op(A) = A` or `op(A) = A'` if `adjA = false` or `adjA = true`, respectively,
and `op(B) = B` or `op(B) = B'` if `adjB = false` or `adjB = true`, respectively.
`A` and `B` are square matrices in Schur forms, and `A` and `-B` must not have
common eigenvalues. `C` contains on output the solution `X`.

_Note:_ This is an enhanced interface to `LAPACK.trsyl!` also covering the case when
`A` and `B` are real matrices and `C` is a complex matrix.
"""
function sylvcs!(A::T1, B::T1, C::Union{T1,T2}; adjA = false, adjB = false) where
     {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   """
   This is a wrapper to the LAPACK.trsylv! function for real matrices A and B,
   based on the Bartels-Stewart Schur form based approach.
   For a complex C, two calls of trsyl! are performed (for the real and imaginary parts).
   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   T = eltype(A)
   TC = eltype(C)
   TR = real(TC)
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")
    cmplx = TC <: Complex
   adjA ? TA = 'T' : TA = 'N'
   adjB ? TB = 'T' : TB = 'N'
   if cmplx
      YR = real(C)
      YR, scale = LAPACK.trsyl!(TA, TB, A, B, YR)
      rmul!(YR, inv(scale))
      YI = imag(C)
      YI, scale = LAPACK.trsyl!(TA, TB, A, B, YI)
      rmul!(YI, inv(scale))
      C[:,:] = YR + im*YI
      return C
   else
      C, scale = LAPACK.trsyl!(TA, TB, A, B, C)
      rmul!(C, inv(scale))
      return C
   end
end
function sylvcs!(A::T, B::T, C::T; adjA = false, adjB = false) where  {T<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}}
   """
   This is a wrapper to the LAPACK.trsylv! function, based on the Bartels-Stewart Schur form based approach.
   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   adjA ? TA = 'C' : TA = 'N'
   adjB ? TB = 'C' : TB = 'N'
   C, scale = LAPACK.trsyl!(TA, TB, A, B, C)
   rmul!(C, inv(scale))
end
"""
    sylvds!(A,B,C; adjA = false, adjB = false)

Solve the discrete Sylvester matrix equation

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

where `op(A) = A` or `op(A) = A'` if `adjA = false` or `adjA = true`, respectively,
and `op(B) = B` or `op(B) = B'` if `adjB = false` or `adjB = true`, respectively.
`A` and `B` are square matrices in Schur forms, and `A` and `-B` must not have
common reciprocal eigenvalues. `C` contains on output the solution `X`.
"""
function sylvds!(A::T1, B::T1, C::Union{T1,T2}; adjA = false, adjB = false) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
#function sylvds!(A::T, B::T, C::T; adjA = false, adjB = false) where  {T<:Union{Array{Float64,2},Array{Float32,2}}}
   """
   An extension of the Bartels-Stewart Schur form based approach is employed.

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

   T = eltype(A)
   TC = eltype(C)
   TR = real(TC)
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")

   # determine the structure of the real Schur form of A
   ba = fill(1,m,1)
   pa = 1
   if m > 1
      d = [diag(A,-1);zeros(1)]
      i = 1
      pa = 0
      while i <= m
         pa += 1
         if d[i] != 0
            ba[pa] = 2
            i += 1
         end
         i += 1
      end
   end
   # determine the structure of the real Schur form of B
   bb = fill(1,n,1)
   pb = 1
   if n > 1
      d = [diag(B,-1);zeros(1)]
      i = 1
      pb = 0
      while i <= n
         pb += 1
         if d[i] != 0
            bb[pb] = 2
            i += 1
         end
         i += 1
      end
   end

   W = fill(zero(eltype(C)),m,2)
   if !adjA && !adjB
      """
      The (K,L)th block of X is determined starting from
      bottom-left corner column by column by

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

      where
                             M
                 R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L) +
                           J=K+1
                             M             L-1
                            SUM { A(K,J) * SUM [X(J,I)*B(I,L)] }.
                            J=K            I=1
      """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = m
          for kk = pa:-1:1
              dk = ba[kk]
              dkk = 1:dk
              i1 = i-dk+1
              k = i1:i
              y = C[k,l]
              if kk < pa
                 ir = i+1:m
                 W1 = A[k,ir]*C[ir,l]
                 y -= W1*B[l,l]
              end
              if ll > 1
                 ic = i1:m
                 W[k,dll] = C[k,il1]*B[il1,l]
                 y -= A[k,ic]*W[ic,dll]
              end
              Z = (kron(transpose(B[l,l]),A[k,k])+I)\(y[:])
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("SingularException: A and -B have common or close reciprocal eigenvalues")
              i -= dk
          end
          j += dl
      end
   elseif !adjA && adjB
         """
         The (K,L)th block of X is determined starting from
         bottom-right corner column by column by

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

         where
                                M
                    R(K,L) = { SUM [A(K,J)*X(J,L)] } * B(L,L)' +
                              J=K+1
                                M              N
                               SUM { A(K,J) * SUM [X(J,I)*B(L,I)'] }.
                               J=K           I=L+1
         """
         j = n
         for ll = pb:-1:1
             dl = bb[ll]
             dll = 1:dl
             il1 = j+1:n
             l = j-dl+1:j
             i = m
             for kk = pa:-1:1
                 dk = ba[kk]
                 dkk = 1:dk
                 i1 = i-dk+1
                 k = i1:i
                 y = C[k,l]
                 if kk < pa
                    ir = i+1:m
                    W1 = A[k,ir]*C[ir,l]
                    y -= W1*B[l,l]'
                 end
                 if ll < pb
                    ic = i1:m
                    W[k,dll] = C[k,il1]*B[l,il1]'
                    y -= A[k,ic]*W[ic,dll]
                 end
                 Z = (kron(B[l,l],A[k,k])+I)\(y[:])
                 isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("SingularException: A and -B have common or close reciprocal eigenvalues")
                 i -= dk
             end
             j -= dl
         end
   elseif adjA && !adjB
      """
      The (K,L)th block of X is determined starting from the
      upper-left corner column by column by

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

      where
                            K-1
                 R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L) +
                            J=1
                             K              L-1
                            SUM A(J,K)' * { SUM [X(J,I)*B(I,L)] }.
                            J=1             I=1
      """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = C[k,l]
              if kk > 1
                 ir = 1:i-1
                 W1 = A[ir,k]'*C[ir,l]
                 y -= W1*B[l,l]
              end
              if ll > 1
                 ic = 1:i1
                 W[k,dll] = C[k,il1]*B[il1,l]
                 y -= A[ic,k]'*W[ic,dll]
              end
              Z = (kron(transpose(B[l,l]),transpose(A[k,k]))+I)\(y[:])
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("SingularException: A and -B have common or close reciprocal eigenvalues")
              i += dk
          end
          j += dl
      end
   elseif adjA && adjB
      """
      The (K,L)th block of X is determined starting from the
      lower-left corner column by column by

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

      where
                            K-1
                 R(K,L) = { SUM [A(J,K)'*X(J,L)] } * B(L,L)' +
                            J=1
                             K               N
                            SUM A(J,K)' * { SUM [X(J,I)*B(L,I)'] }.
                            J=1            I=L+1
      """
      j = n
      for ll = pb:-1:1
          dl = bb[ll]
          dll = 1:dl
          il1 = j+1:n
          l = j-dl+1:j
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = C[k,l]
              if kk > 1
                 ir = 1:i-1
                 W1 = A[ir,k]'*C[ir,l]
                 y -= W1*B[l,l]'
              end
              if ll < pb
                 ic = 1:i1
                 W[k,dll] = C[k,il1]*B[l,il1]'
                 y -= A[ic,k]'*W[ic,dll]
              end
              Z = (kron(B[l,l],transpose(A[k,k]))+I)\(y[:])
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("SingularException: A and -B have common or close reciprocal eigenvalues")
              i += dk
          end
          j -= dl
      end
   end
   return C
end
function sylvds!(A::T, B::T, C::T; adjA = false, adjB = false) where {T<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}}
   """
   An extension of the Bartels-Stewart Schur form based approach is employed.

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

   W = fill(zero(eltype(C)),m,1)
   if !adjA && !adjB
      """
      The (K,L)th element of X is determined starting from
      bottom-left corner column by column by

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

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

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

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

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

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

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

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

Solve the generalized Sylvester matrix equation

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

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

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

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

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

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

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

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

   T = eltype(A)
   TE = eltype(E)
   TR = real(TE)
   T == TR || error("TypeError: for real part of E expected Type{$T}, got Type{$TR}")

   # determine the structure of the generalized real Schur form of (A,C)
   ba = fill(1,m,1)
   pa = 1
   if m > 1
      CASchur ? d = [diag(C,-1);zeros(1)] : d = [diag(A,-1);zeros(1)]
      i = 1
      pa = 0
      while i <= m
         pa += 1
         if d[i] != 0
            ba[pa] = 2
            i += 1
         end
         i += 1
      end
   end
   # determine the structure of the generalized real Schur form of (B,D)
   bb = fill(1,n,1)
   pb = 1
   if n > 1
      DBSchur ? d = [diag(D,-1);zeros(1)] : d = [diag(B,-1);zeros(1)]
      i = 1
      pb = 0
      while i <= n
         pb += 1
         if d[i] != 0
            bb[pb] = 2
            i += 1
         end
         i += 1
      end
   end

   WB = fill(zero(eltype(E)),m,2)
   WD = fill(zero(eltype(E)),m,2)
   if !adjAC && !adjBD
      """
      The (K,L)th block of X is determined starting from
      bottom-left corner column by column by

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

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

                             M
                          { SUM [C(K,J)*X(J,L)] } * D(L,L) +
                           J=K+1
                             M             L-1
                            SUM { C(K,J) * SUM [X(J,I)*D(I,L)] }.
                            J=K            I=1
      """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = m
          for kk = pa:-1:1
              dk = ba[kk]
              dkk = 1:dk
              i1 = i-dk+1
              k = i1:i
              y = E[k,l]
              if kk < pa
                 ir = i+1:m
                 W1 = A[k,ir]*E[ir,l]
                 y -= W1*B[l,l]
                 W1 = C[k,ir]*E[ir,l]
                 y -= W1*D[l,l]
              end
              if ll > 1
                 ic = i1:m
                 WB[k,dll] = E[k,il1]*B[il1,l]
                 WD[k,dll] = E[k,il1]*D[il1,l]
                 y -= (A[k,ic]*WB[ic,dll] + C[k,ic]*WD[ic,dll])
              end
              Z = (kron(transpose(B[l,l]),A[k,k])+kron(transpose(D[l,l]),C[k,k]))\(y[:])
              isfinite(maximum(abs.(Z))) ? E[k,l] = Z : throw("SingularException: A-λC and D+λB have common or close eigenvalues")
              i -= dk
          end
          j += dl
      end
   elseif !adjAC && adjBD
         """
          The (K,L)th block of X is determined starting from
          bottom-right corner column by column by

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

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

                               M
                            { SUM [C(K,J)*X(J,L)] } * D(L,L)' +
                             J=K+1
                               M              N
                              SUM { C(K,J) * SUM [X(J,I)*D(L,I)'] }.
                              J=K           I=L+1
         """
         j = n
         for ll = pb:-1:1
             dl = bb[ll]
             dll = 1:dl
             il1 = j+1:n
             l = j-dl+1:j
             i = m
             for kk = pa:-1:1
                 dk = ba[kk]
                 dkk = 1:dk
                 i1 = i-dk+1
                 k = i1:i
                 y = E[k,l]
                 if kk < pa
                    ir = i+1:m
                    W1 = A[k,ir]*E[ir,l]
                    y -= W1*B[l,l]'
                    W2 = C[k,ir]*E[ir,l]
                    y -= W2*D[l,l]'
                 end
                 if ll < pb
                    ic = i1:m
                    WB[k,dll] = E[k,il1]*B[l,il1]'
                    WD[k,dll] = E[k,il1]*D[l,il1]'
                    y -= (A[k,ic]*WB[ic,dll]+C[k,ic]*WD[ic,dll])
                 end
                 Z = (kron(B[l,l],A[k,k])+kron(D[l,l],C[k,k]))\(y[:])
                 isfinite(maximum(abs.(Z))) ? E[k,l] = Z : throw("SingularException: A-λC and D+λB have common or close eigenvalues")
                 i -= dk
             end
             j -= dl
         end
   elseif adjAC && !adjBD
      """
      The (K,L)th block of X is determined starting from the
      upper-left corner column by column by

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

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

                            K-1
                          { SUM [C(J,K)'*X(J,L)] } * D(L,L) +
                            J=1
                             K              L-1
                            SUM C(J,K)' * { SUM [X(J,I)*D(I,L)] }.
                            J=1             I=1
      """
      j = 1
      for ll = 1:pb
          dl = bb[ll]
          dll = 1:dl
          il1 = 1:j-1
          j1 = j+dl-1
          l = j:j1
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = E[k,l]
              if kk > 1
                 ir = 1:i-1
                 W1 = A[ir,k]'*E[ir,l]
                 y -= W1*B[l,l]
                 W2 = C[ir,k]'*E[ir,l]
                 y -= W2*D[l,l]
              end
              if ll > 1
                 ic = 1:i1
                 WB[k,dll] = E[k,il1]*B[il1,l]
                 y -= A[ic,k]'*WB[ic,dll]
                 WD[k,dll] = E[k,il1]*D[il1,l]
                 y -= C[ic,k]'*WD[ic,dll]
              end
              Z = (kron(transpose(B[l,l]),transpose(A[k,k]))+kron(transpose(D[l,l]),transpose(C[k,k])))\(y[:])
              isfinite(maximum(abs.(Z))) ? E[k,l] = Z : throw("SingularException: A-λC and D+λB have common or close eigenvalues")
              i += dk
          end
          j += dl
      end
   elseif adjAC && adjBD
      """
      The (K,L)th block of X is determined starting from
      upper-rght corner column by column by

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

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

                            K-1
                          { SUM [C(J,K)'*X(J,L)] } * D(L,L)' +
                            J=1
                             K               N
                            SUM C(J,K)' * { SUM [X(J,I)*D(L,I)'] }.
                            J=1            I=L+1
      """
      j = n
      for ll = pb:-1:1
          dl = bb[ll]
          dll = 1:dl
          il1 = j+1:n
          l = j-dl+1:j
          i = 1
          for kk = 1:pa
              dk = ba[kk]
              dkk = 1:dk
              i1 = i+dk-1
              k = i:i1
              y = E[k,l]
              if kk > 1
                 ir = 1:i-1
                 W1 = A[ir,k]'*E[ir,l]
                 y -= W1*B[l,l]'
                 W2 = C[ir,k]'*E[ir,l]
                 y -= W2*D[l,l]'
              end
              if ll < pb
                 ic = 1:i1
                 WB[k,dll] = E[k,il1]*B[l,il1]'
                 WD[k,dll] = E[k,il1]*D[l,il1]'
                 y -= (A[ic,k]'*WB[ic,dll] + C[ic,k]'*WD[ic,dll])
              end
              Z = (kron(B[l,l],transpose(A[k,k]))+kron(D[l,l],transpose(C[k,k])))\(y[:])
              isfinite(maximum(abs.(Z))) ? E[k,l] = Z : throw("SingularException: A-λC and D+λB have common or close eigenvalues")
              i += dk
          end
          j -= dl
      end
   end
   return E
end
function gsylvs!(A::T, B::T, C::T, D::T, E::T; adjAC = false, adjBD = false, CASchur = false, DBSchur = false) where {T<:Union{Array{Complex{Float64},2},Array{Complex{Float32},2}}}
   """
   An extension proposed in [1] of the Bartels-Stewart Schur form based approach [2] is employed.

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

   WB = fill(zero(eltype(E)),m,1)
   WD = fill(zero(eltype(E)),m,1)
   if !adjAC && !adjBD
      """
      The (K,L)th element of X is determined starting from
      bottom-left corner column by column by

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

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

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

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

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

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

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

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

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

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

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

                            K-1
                          { SUM [C(J,K)'*X(J,L)] } * D(L,L)' +
                            J=1
                             K               N
                            SUM C(J,K)' * { SUM [X(J,I)*D(L,I)'] }.
                            J=1            I=L+1
      """
      for l = n:-1:1
          ll = l:l
          il1 = l+1:n
          for k = 1:m
              kk = k:k
              y = E[k,l]
              if k > 1
                 ir = 1:k-1
                 W1 = A[ir,kk]'*E[ir,ll]
                 W2 = C[ir,kk]'*E[ir,ll]
                 y -= (W1[1]*B[l,l]'+W2[1]*D[l,l]')
              end
              if l < n
                 ic = 1:m
                 ZB = E[kk,il1]*B[ll,il1]'
                 ZD = E[kk,il1]*D[ll,il1]'
                 WB[k,1] = ZB[1]
                 WD[k,1] = ZD[1]
                 TB = A[ic,kk]'*WB[ic,1]
                 TD = C[ic,kk]'*WD[ic,1]
                 y -= (TB[1]+TD[1])
              end
              Z = y/(B[l,l]'*A[k,k]'+D[l,l]'*C[k,k]')
              isfinite(Z) ? E[k,l] = Z : throw("SingularException: A-λC and D+λB have common or close eigenvalues")
          end
      end
   end
   return E
end
"""
    (X,Y) = sylvsyss!(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F,

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

_Note:_ This is an enhanced interface to the `LAPACK.tgsyl!` function to also cover the case when
`A`, `B`, `D` and `E` are real matrices and `C` and `F` are complex matrices.
"""
function sylvsyss!(A::T, B::T, C::T, D::T, E::T, F::T) where {T<:Union{Array{Complex{Float64},2},Array{Complex{Float32},2}}}
   """
   This is a wrapper to the complex LAPACK.tgsyl! function with `trans = 'N'`.
   """
   C, F, scale =  tgsyl!('N',A,B,C,D,E,F)
   return rmul!(C,inv(scale)), rmul!(F,inv(-scale))
end
function sylvsyss!(A::T1, B::T1, C::Union{T1,T2}, D::T1, E::T1, F::Union{T1,T2}) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   """
   This is an enhanced interface to the real LAPACK.tgsyl! function with `trans = 'N'`.
   For a complex pair (C,F), two calls of tgsyl! are performed (for the real and imaginary parts).
   """
   T = eltype(A)
   TC = promote_type(eltype(C),eltype(F))
   TR = real(TC)
   T == TR || error("TypeError: for real parts of (C,F) expected Type{$T}, got Type{$TR}")
   cmplx = TC <: Complex
   if cmplx
      XR = real(C)
      XI = imag(C)
      YR = real(F)
      YI = imag(F)
      XR, YR, scale =  tgsyl!('N',A,B,XR,D,E,YR)
      rmul!(XR,inv(scale))
      rmul!(YR,inv(-scale))
      XI, YI, scale =  tgsyl!('N',A,B,XI,D,E,YI)
      rmul!(XI,inv(scale))
      rmul!(YI,inv(-scale))
      C = XR+im*XI
      F = YR+im*YI
      return C[:,:], F[:,:]
   else
      C, F, scale =  tgsyl!('N',A,B,C,D,E,F)
      return rmul!(C,inv(scale)), rmul!(F,inv(-scale))
   end
end
"""
    (X,Y) = dsylvsyss!(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

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

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

_Note:_ This is an enhanced interface to the `LAPACK.tgsyl!` function to also cover the case when
`A`, `B`, `D` and `E` are real matrices and `C` and `F` are complex matrices.
"""
function dsylvsyss!(A::T, B::T, C::T, D::T, E::T, F::T) where {T<:Union{Array{Complex{Float64},2},Array{Complex{Float32},2}}}
   """
   This is an interface to the complex LAPACK.tgsyl! function with `trans = 'C'`.
   """
   MF = -F
   E, F, scale =  tgsyl!('C',A,B,C,D,E,MF)
   F = MF
   return rmul!(C[:,:],inv(scale)), rmul!(F[:,:],inv(scale))
end
function dsylvsyss!(A::T1, B::T1, C::Union{T1,T2}, D::T1, E::T1, F::Union{T1,T2}) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   """
   This is an interface to the real LAPACK.tgsyl! function with `trans = 'T'`.
   For a complex pair (C,F), two calls of tgsyl! are performed (for the real and imaginary parts).
   """
   T = eltype(A)
   TC = promote_type(eltype(C),eltype(F))
   TR = real(TC)
   T == TR || error("TypeError: for real parts of (C,F) expected Type{$T}, got Type{$TR}")
   cmplx = TC <: Complex
   if cmplx
      XR = real(C)
      XI = imag(C)
      YR = -real(F)
      YI = -imag(F)
      XR, YR, scale =  tgsyl!('T',A,B,XR,D,E,YR)
      rmul!(XR,inv(scale))
      rmul!(YR,inv(scale))
      XI, YI, scale =  tgsyl!('T',A,B,XI,D,E,YI)
      rmul!(XI,inv(scale))
      rmul!(YI,inv(scale))
      C = XR+im*XI
      F = YR+im*YI
      return C[:,:], F[:,:]
   else
      MF = -F
      C, MF, scale =  tgsyl!('T',A,B,C,D,E,MF)
      F = MF
      return rmul!(C[:,:],inv(scale)), rmul!(F[:,:],inv(scale))
   end
end


text:
function find_arrest(S₀, S₁; rtol = eps(), artol = sqrt(eps()))
    grid = S₁.grid
    Λᵤ = Interpolator(S₀, S₁)
    ws = ApproximationGrids.weights(grid)
    avars, kvars, D₀ = initialize_asymptotics(S₁)

    K = nodes(grid)
    S = kvars.svars.S
    υ = kvars.svars.υ
    w = kvars.svars.w
    d = dimensionality(S₁.f.liquid)
    Z = similar(w)
    T = eltype(K)

    i  = argmax(Λᵤ.A)
    u₀ = zero(T)
    u₁ = log(floatmax(T)) / abs(Λᵤ.A[i])
    uₘ = u₁ / 2

    g = LiquidsDynamics.FixedPoint(Z, kvars, D₀)

    while true
        interpolate!(S, Λᵤ, uₘ)
        LiquidsDynamics.asymptotics_weights!(w, ws, K, S, d, grid)
        ζ = υ * sum(w)

        asymptotics!(g, avars, ζ; rtol = artol)

        if iszero(avars.ζ∞[])
            u₀ = uₘ
        else
            u₁ = uₘ
        end

        if isapprox(u₀, u₁, rtol = rtol, nans = true)
            if iszero(avars.ζ∞[])
                interpolate!(S, Λᵤ, u₁)
                LiquidsDynamics.asymptotics_weights!(w, ws, K, S, d, grid)
                asymptotics!(g, avars, υ * sum(w); rtol = artol)
            end
            break
        end

        uₘ = (u₀ + u₁) / 2
    end

    return u₀, u₁, avars.ζ∞[]
end

ext: jl


text:
# Autogenerated wrapper script for JpegTurbo_jll for x86_64-w64-mingw32
export cjpeg, djpeg, jpegtran, libjpeg, libturbojpeg

JLLWrappers.@generate_wrapper_header("JpegTurbo")
JLLWrappers.@declare_executable_product(cjpeg)
JLLWrappers.@declare_executable_product(djpeg)
JLLWrappers.@declare_executable_product(jpegtran)
JLLWrappers.@declare_library_product(libjpeg, "libjpeg-62.dll")
JLLWrappers.@declare_library_product(libturbojpeg, "libturbojpeg.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cjpeg,
        "bin\\cjpeg.exe",
    )

    JLLWrappers.@init_executable_product(
        djpeg,
        "bin\\djpeg.exe",
    )

    JLLWrappers.@init_executable_product(
        jpegtran,
        "bin\\jpegtran.exe",
    )

    JLLWrappers.@init_library_product(
        libjpeg,
        "bin\\libjpeg-62.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libturbojpeg,
        "bin\\libturbojpeg.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#MDDatasets: Support broadcast of DataRS types
#-------------------------------------------------------------------------------

#==result_type: Figure out result type for a function call with given arguments
===============================================================================#
#=NOTES:
 - result_type never returns concrete DataRS{T}.  Only abstract "DataRS".

 - The default return types specified here are not necessarily correct.
   Methods should be overwritten for special functions.
   Sadly: can't currently dispatch on function.  Need to think about this.
=#

result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:DataRS, T2<:DataRS} = DataRS
result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:DataRS, T2} = DataRS
result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1, T2<:DataRS} = DataRS

#Reducting functions
#-------------------------------------------------------------------------------
result_type(CT::CastTypeRed, fn::Function, ::Type{DataRS{T1}}) where T1<:Number = T1
result_type(CT::CastTypeRed, fn::Function, ::Type{DataRS{T}}) where T = DataRS

#==eltype/valtype: Figure out result type for a function call with given arguments
===============================================================================#

#Element type of an operation with a DataRS (result: DataRS):
Base.eltype(fn::Function, d1::DataRS) = result_type(fn, eltype(d1))

Base.eltype(fn::Function, d1::DataRS, d2::DataRS) =
	result_type(fn, eltype(d1), eltype(d2))
Base.eltype(fn::Function, d1::DataRS, d2) =
	result_type(fn, eltype(d1), typeof(d2))
Base.eltype(fn::Function, d1, d2::DataRS) =
	result_type(fn, typeof(d1), eltype(d2))

#Reducting functions
#-------------------------------------------------------------------------------
#Element type afer reducing fn on Array{Number}:
Base.eltype(CT::CastTypeRed, fn::Function, a::Vector{T}) where T<:Number = result_type(CT, fn, typeof(a))

#Element type afer reducing fn on Array{DataF1}:
function Base.eltype(CT::CastTypeRed, fn::Function, a::Array{DataF1})
	result = Set{DataType}()
	for elem in a
		RT = result_type(CT, fn, typeof(elem))
		push!(result, RT)
	end
	return promote_type(result...)
end

#Element type afer reducing fn on Array{DataRS}:
function Base.eltype(CT::CastTypeRed, fn::Function, a::Vector{DataRS})
	result = Set{DataType}()
	RT = result_type(CT, fn, typeof(a[1]))
	if DataRS == RT; return DataRS; end
	for elem in a
		RT = result_type(CT, fn, typeof(elem))
		push!(result, RT)
	end
	return promote_type(result...)
end

#==Broadcast tools for fn(DataMD) - Can broadcast on 1 argument
===============================================================================#

#Broadcast functions capable of operating directly on 1 base type (Number):
#-------------------------------------------------------------------------------
#fn(DataRS) - core: fn(Number):
function broadcastMD(CT::CastType1{Number,1}, fn::Function, d::DataRS, args...; kwargs...)
	result = DataRS{eltype(fn, d)}(d.sweep)
	for i in 1:length(d.sweep)
		result.elem[i] = broadcastMD(CT, fn, d.elem[i], args...; kwargs...)
	end
	return result
end

#Data reducing fn(DataRS{Number}) - core: fn(Number):
broadcastMD(CT::CastTypeRed1{Number,1}, fn::Function, d::DataRS{T}, args...; kwargs...) where T<:Number =
	fn(d.elem, args...; kwargs...)

#Data reducing fn(DataRS) - core: fn(Number):
function broadcastMD(CT::CastTypeRed1{Number,1}, fn::Function, d::DataRS, args...; kwargs...)
	RT = eltype(CT, fn, d.elem)
	result = DataRS{RT}(d.sweep)
	for i in 1:length(d.sweep)
		result.elem[i] = broadcastMD(CT, fn, d.elem[i], args...; kwargs...)
	end
	return result
end

#Broadcast functions capable of operating only on a dataF1 value:
#-------------------------------------------------------------------------------
#fn(DataRS) - core: fn(DataF1):
function broadcastMD(CT::CastType1{DataF1,1}, fn::Function, d::DataRS{T}, args...; kwargs...) where T
	result = DataRS{T}(d.sweep)
	for i in 1:length(d.sweep)
		result.elem[i] = broadcastMD(CT, fn, d.elem[i], args...; kwargs...)
	end
	return result
end
#fn(???, DataRS) - core: fn(???, DataF1):
function broadcastMD(CT::CastType1{DataF1,2}, fn::Function, dany1, d::DataRS{T}, args...; kwargs...) where T
	result = DataRS{T}(d.sweep)
	for i in 1:length(d.sweep)
		result.elem[i] = broadcastMD(CT, fn, dany1, d.elem[i], args...; kwargs...)
	end
	return result
end

#Data reducing fn(DataRS) - core: fn(Number):
broadcastMD(CT::CastTypeRed1{DataF1,1}, fn::Function, d::DataRS{T}, args...; kwargs...) where T<:Number =
	fn(DataF1(d.sweep.v, d.elem), args...; kwargs...)

#Data reducing fn(DataRS) - core: fn(Number):
function broadcastMD(CT::CastTypeRed1{DataF1,1}, fn::Function, d::DataRS, args...; kwargs...)
	RT = eltype(CT, fn, d.elem)
	result = DataRS{RT}(d.sweep)
	for i in 1:length(d.sweep)
		result.elem[i] = broadcastMD(CT, fn, d.elem[i], args...; kwargs...)
	end
	return result
end

#==Broadcast tools for fn(DataMD, DataMD) - Can broadcast on 2 arguments
===============================================================================#

#Broadcast functions capable of operating directly on base types (Number, Number):
#-------------------------------------------------------------------------------
#fn(DataRS, DataRS) - core: fn(Number, Number):
function broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function,
	d1::DataRS, d2::DataRS, args...; kwargs...)
	if d1.sweep != d2.sweep
		msg = "Sweeps do not match (not yet supported):"
		throw(ArgumentError(string(msg, "\n", d1.sweep, "\n", d2.sweep)))
	end
	result = DataRS{eltype(fn, d1, d2)}(d1.sweep)
	for i in 1:length(d1.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1.elem[i], d2.elem[i], args...; kwargs...)
	end
	return result
end

#fn(DataRS, DataF1/Number) - core: fn(Number, Number):
function broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function,
	d1::DataRS, d2::DF1_Num, args...; kwargs...)
	result = DataRS{eltype(fn, d1, d2)}(d1.sweep)
	for i in 1:length(d1.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1.elem[i], d2, args...; kwargs...)
	end
	return result
end

#fn(DataF1/Number, DataRS) - core: fn(Number, Number):
function broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function,
	d1::DF1_Num, d2::DataRS, args...; kwargs...)
	result = DataRS{eltype(fn, d1, d2)}(d2.sweep)
	for i in 1:length(d2.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1, d2.elem[i], args...; kwargs...)
	end
	return result
end

#Broadcast functions capable of operating on DataF1 values:
#-------------------------------------------------------------------------------
const DF1_DRS = Union{DataF1,DataRS}
#fn(DataRS, DataRS) - core: fn(DataF1, DataF1):
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataRS{T1}, d2::DataRS{T2}, args...; kwargs...) where {T1<:DF1_DRS,T2<:DF1_DRS}
	if d1.sweep != d2.sweep
		msg = "Sweeps do not match (not yet supported):"
		throw(ArgumentError(string(msg, "\n", d1.sweep, "\n", d2.sweep)))
	end
	result = DataRS{eltype(fn, d1, d2)}(d1.sweep)
	for i in 1:length(d1.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1.elem[i], d2.elem[i], args...; kwargs...)
	end
	return result
end
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataRS{T1}, d2::DataRS{T2}, args...; kwargs...) where {T1<:Number,T2<:Number}

	if d1.sweep.id != d2.sweep.id
		msg = "Sweep ids do not match:"
		throw(ArgumentError(string(msg, "\n", d1.sweep, "\n", d2.sweep)))
	end
	return fn(DataF1(d1.sweep.v, d1.elem), DataF1(d2.sweep.v, d2.elem), args...; kwargs...)
end
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataRS{T1}, d2::DataMD, args...; kwargs...) where T1<:Number
	return broadcastMD(CT, fn, DataF1(d1.sweep.v, d1.elem), d2, args...; kwargs...)
end
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataMD, d2::DataRS{T2}, args...; kwargs...) where T2<:Number
	return broadcastMD(CT, fn, d1, DataF1(d2.sweep.v, d2.elem), args...; kwargs...)
end

#fn(DataRS, DataF1) - core: fn(DataF1, DataF1):
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataRS{TRS}, d2::DataF1, args...; kwargs...) where TRS<:DF1_DRS
	result = DataRS{eltype(fn, d1, d2)}(d1.sweep)
	for i in 1:length(d1.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1.elem[i], d2, args...; kwargs...)
	end
	return result
end
#fn(DataF1, DataRS) - core: fn(DataF1, DataF1):
function broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function,
	d1::DataF1, d2::DataRS{TRS}, args...; kwargs...) where TRS<:DF1_DRS
	result = DataRS{eltype(fn, d1, d2)}(d2.sweep)
	for i in 1:length(d2.sweep)
		result.elem[i] = broadcastMD(CT, fn, d1, d2.elem[i], args...; kwargs...)
	end
	return result
end
#Last line


text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setpath.jl");

using Preprocessing, Potential

# Read xyz file to obtain configurations
filename = "/Users/ngoccuongnguyen/Dropbox (MIT)/MDP/data/LJCluster/curated_lj_cluster.xyz";
species = ["Ar"];
config = Preprocessing.readEXTXYZ(filename, species);

# shift atom positions to put them in the box [0,10]x[0,10]x[0,10]
config.x = config.x .+ 5.0

# cumalative sum of numbers of atoms
natom = [0; cumsum(config.natom[:])];

ci = 1; # configuration index in the data
x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
t = config.t[:,(natom[ci]+1):natom[ci+1]]                        # atom types of configuration ci
q = [];                                   # atom charges of configuration ci

# 3 vectors define the simulation box
a = [10.0; 0.0; 0.0]
b = [0.0; 10.0; 0.0]
c = [0.0; 0.0; 10.0]

# periodic boundary conditions
pbc = [0; 0; 0]

# LJ parameters
epsilon = 1.0;
sigma = 1.0;
A = 2*epsilon*sigma^12;
B = 2*epsilon*sigma^6;
rcut = 4.0;

# define a list of potentials
potential = Array{Any}(undef, 1)
potential[1] = Potential.addpotential("pair", "nonbonded", [0], "ljpotential", [A; B], rcut);
for i = 1:length(potential)
    include(potential[i].potentialfunction * ".jl");  # include the potential file
end

# compute total energy, peratom energies, and forces
rcutmax = rcut; eta = [0]; kappa = [0];
etot, e, f, v, vtot = Potential.empiricalpotential(config.x[:,1:13], q, t, a, b, c, pbc, rcutmax, eta, kappa, potential);

# using LAMMPS.jl to compute energy and forces
unitstyle = "lj"
pair_style = "lj/cut $rcut"
pair_coeff = ["* * $epsilon $sigma"]
pe, forces, id, pos, type = Potential.lammpspotential(x, t, a, b, c, pbc, unitstyle, pair_style, pair_coeff)

# compare with LAMMPS.jl
print("Error in potential energy: "); display(abs(etot - pe));
print("Error in force: "); display(maximum(abs.(f[:,id] .- forces)));
print("Error in stress: "); display(maximum(abs.(vtot .- config.stress[:,ci])));

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Dates
using AxisArrays
using ClimatePlots
using ClimateTools
using PyPlot
using Test

@testset "Mapping" begin

# Shapefile test
filenc = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")

C = load(filenc, "tas")
# msk = inpolygrid(C.longrid, C.latgrid, P)
# status, figh = ClimatePlots.mapclimgrid(C); @test status == true; PyPlot.close()

regions = ["World", "Mollweide", "eck1", "eck2", "eck3", "eck4", "eck5", "eck6","Canada", "Quebec", "USA","NorthAmerica", "SouthAmerica", "Arctic", "Antarctic", "Africa", "Europe", "Asia","Mauricie", "auto"]

for iregion in regions
    status, figh = ClimatePlots.contourf(C, region=iregion);@test status == true; PyPlot.close()
    status, figh = ClimatePlots.contour(C, region=iregion);@test status == true; PyPlot.close()
    status, figh = ClimatePlots.pcolormesh(C, region=iregion);@test status == true; PyPlot.close()
end

# precip
C = load(filenc, "pr", data_units="mm") + 2.0

for iregion in regions
    status, figh = ClimatePlots.contourf(C, region=iregion);@test status == true; PyPlot.close()
    status, figh = ClimatePlots.contour(C, region=iregion);@test status == true; PyPlot.close()
    status, figh = ClimatePlots.pcolormesh(C, region=iregion);@test status == true; PyPlot.close()
end

# ua wind
C = load(filenc, "ua")
status, figh = ClimatePlots.contourf(C, level = 3);@test status == true; PyPlot.close() # given level
status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close() # feeding a 4D field

# DUMMY MAPS
lon = collect(-180.0:2.0:180.0)
lat = collect(-90.0:2.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = DateTime(2003,1,1):Day(1):DateTime(2003,01,31)
data = randn(181, 91, 31)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()
status, figh = ClimatePlots.contour(C);@test status == true; PyPlot.close()
status, figh = ClimatePlots.pcolormesh(C);@test status == true; PyPlot.close()

status, figh = ClimatePlots.contourf(C, center_cs=true, caxis=[200, 250]);@test status == true; PyPlot.close()
status, figh = ClimatePlots.contour(C, center_cs=true, caxis=[200, 250]);@test status == true; PyPlot.close()
status, figh = ClimatePlots.pcolormesh(C, center_cs=true, caxis=[200, 250]);@test status == true; PyPlot.close()

status, figh = ClimatePlots.contourf(C, center_cs=true);@test status == true; PyPlot.close()
status, figh = ClimatePlots.contour(C, center_cs=true);@test status == true; PyPlot.close()
status, figh = ClimatePlots.pcolormesh(C, center_cs=true);@test status == true; PyPlot.close()

lon = collect(-180.0:2.0:180.0)
lat = collect(-90.0:2.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = Dates.year(DateTime(2003, 01, 01)):Dates.year(1):Dates.year(DateTime(2004, 01, 01))
data = randn(181, 91, 2)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()

C = ClimateTools.ClimGrid(axisdata, variable = "clt", typeofvar="clt", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()

C = ClimateTools.ClimGrid(axisdata, variable = "rlut", typeofvar="rlut", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()

C = ClimateTools.ClimGrid(axisdata, variable = "hurs", typeofvar="hurs", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()

C = ClimateTools.ClimGrid(axisdata, variable = "other", typeofvar="other", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.contourf(C);@test status == true; PyPlot.close()

lon = collect(-180.0:2.0:180.0)
lat = collect(-90.0:2.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = DateTime(2003, 01, 01):Day(1):DateTime(2004, 01, 01)
data = randn(181, 91, 366)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = plot(C); @test status == true; PyPlot.close()
status, figh = plot(C, xlimit=[2 5]); @test status==true;PyPlot.close()
status, figh = plot(C, label = "dummy", titlestr="dummytest", gridfig=true); @test status == true; PyPlot.close()

status, figh = hist(C); @test status == true; PyPlot.close()
status, figh = hist(C, ylimit=[2, 5]); @test status == true; PyPlot.close()
status, figh = hist(C, label = "dummy", titlestr="dummytest", gridfig=false); @test status == true; PyPlot.close()

end


text:
#=
Queue
=#

"""
Space complexity: O(n)
"""
struct SimpleQueue{T}
    l::SimpleLinkedList{T}

    SimpleQueue{T}() where T = new{T}(SimpleLinkedList{T}())
end

#=
Core methods
=#

"""
Time complexity: O(1)
"""
function Base.push!(sq::SimpleQueue{T}, item::T) where T
    push!(sq.l, item)
    return sq
end

"""
Time complexity: O(1)
"""
Base.pop!(sq::SimpleQueue) = popfirst!(sq.l)


text:
export DoubleCategory, FreeDoubleCategory, HomH, HomV, Hom2, idH, idV,
  id2, id2V, id2H, composeH, composeV, ⋆, MonoidalDoubleCategory,
  SymmetricMonoidalDoubleCategory, FreeSymmetricMonoidalDoubleCategory,
  braidH, braidV, σH, σV

# Double category
#################

""" Theory of (strict) *double categories*
"""
@theory DoubleCategory{Ob,HomV,HomH,Hom2} begin
  # """ Object in a category """
  Ob::TYPE
  """ Vertical Morphism in a double category """
  HomV(dom::Ob, codom::Ob)::TYPE
  """ Horizontal Morphism in a double category """
  HomH(dom::Ob, codom::Ob)::TYPE
  """ 2-cell in a double category """
  Hom2(top::HomH(A,B),
       bottom::HomH(C,D),
       left::HomV(A,C),
       right::HomV(B,D))::TYPE ⊣ (A::Ob, B::Ob, C::Ob, D::Ob)
  @op begin
    (→) := HomH
    (↓) := HomV
    (⇒) := Hom2
    (⋆) := composeH
    (⋅) := composeV
  end

  idH(A::Ob)::(A → A) ⊣ (A::Ob)
  idV(A::Ob)::(A ↓ A) ⊣ (A::Ob)
  composeH(f::(A → B), g::(B → C))::(A → C) ⊣ (A::Ob, B::Ob, C::Ob)
  composeV(f::(A ↓ B), g::(B ↓ C))::(A ↓ C) ⊣ (A::Ob, B::Ob, C::Ob)

  # Category axioms for Horizontal morphisms
  ((f ⋆ g) ⋆ h == f ⋆ (g ⋆ h)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B), g::(B → C), h::(C → D)))
  f ⋆ idH(B) == f ⊣ (A::Ob, B::Ob, f::(A → B))
  idH(A) ⋆ f == f ⊣ (A::Ob, B::Ob, f::(A → B))

  # Category axioms for Vertical morphisms
  ((f ⋅ g) ⋅ h == f ⋅ (g ⋅ h)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, f::(A ↓ B), g::(B ↓ C), h::(C ↓ D)))
  f ⋅ idV(B) == f ⊣ (A::Ob, B::Ob, f::(A ↓ B))
  idV(A) ⋅ f == f ⊣ (A::Ob, B::Ob, f::(A ↓ B))

  # identity two cell on 1 object
  id2(X::Ob)::Hom2(idH(X), idH(X), idV(X), idV(X)) ⊣ (X::Ob)
  # identity two cell for vertical composition
  id2V(f::(X→Y))::Hom2(f, f, idV(X), idV(Y)) ⊣ (X::Ob, Y::Ob)
  # identity two cell for horizontal composition
  id2H(f::(X↓Y))::Hom2(idH(X), idH(Y), f, f) ⊣ (X::Ob, Y::Ob)

  # Vertical composition of 2-cells
  composeV(
    α::Hom2(t,b,l,r),
    β::Hom2(b,b2,l2,r2)
  )::Hom2(t, b2, l⋅l2, r⋅r2)  ⊣ (A::Ob, B::Ob, X::Ob, Y::Ob, C::Ob, D::Ob,
                                 t::(A→B), b::(X→Y), l::(A↓X), r::(B↓Y),
                                 b2::(C→D), l2::(X↓C), r2::(Y↓D))

  # Horizontal composition of 2-cells
  composeH(
    α::Hom2(t,b,l,r),
    β::Hom2(t2,b2,r,r2)
  )::Hom2(t⋆t2, b⋆b2, l, r2)  ⊣ (A::Ob, B::Ob, X::Ob, Y::Ob, C::Ob, D::Ob,
                                 t::(A→X), b::(B→Y), l::(A↓B), r::(X↓Y),
                                 t2::(X→C), b2::(Y→D), r2::(C↓D))
end

# Convenience constructors
composeH(αs::Vector) = foldl(composeH, αs)
composeV(αs::Vector) = foldl(composeV, αs)
composeH(α, β, γ, αs...) = composeH([α, β, γ, αs...])
composeV(α, β, γ, αs...) = composeV([α, β, γ, αs...])

""" Syntax for a double category.

Checks domains of morphisms but not 2-morphisms.
"""
@syntax FreeDoubleCategory{ObExpr,HomVExpr, HomHExpr,Hom2Expr} DoubleCategory begin
  compose(f::HomV, g::HomV) = associate(new(f,g; strict=true))
  compose(f::HomH, g::HomH) = associate(new(f,g; strict=true))
  composeH(α::Hom2, β::Hom2) = associate(new(α,β))
  composeV(α::Hom2, β::Hom2) = associate(new(α,β))
end

function show_unicode(io::IO, expr::Hom2Expr{:composeV}; kw...)
  Syntax.show_unicode_infix(io, expr, "⋅"; kw...)
end
function show_unicode(io::IO, expr::Hom2Expr{:composeH}; kw...)
  Syntax.show_unicode_infix(io, expr, "*"; kw...)
end

function show_latex(io::IO, expr::Hom2Expr{:composeV}; kw...)
  Syntax.show_latex_infix(io, expr, "\\cdot"; kw...)
end
function show_latex(io::IO, expr::Hom2Expr{:composeH}; kw...)
  Syntax.show_latex_infix(io, expr, "\\star"; kw...)
end

# Monoidal category
###################

""" Theory of *monoidal double categories*

To avoid associators and unitors, we assume the monoidal double category is
*strict*. By the coherence theorem there is no loss of generality, but we may
add a theory for weak monoidal categories later.
"""
@theory MonoidalDoubleCategory{Ob,HomV,HomH,Hom2} <: DoubleCategory{Ob,HomV,HomH,Hom2} begin
  otimes(A::Ob, B::Ob)::Ob
  otimes(f::(A → B), g::(C → D))::((A ⊗ C) → (B ⊗ D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  otimes(f::(A ↓ B), g::(C ↓ D))::((A ⊗ C) ↓ (B ⊗ D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  otimes(f::Hom2(t,b,l,r),g::Hom2(t′,b′,l′,r′))::Hom2(t⊗t′,b⊗b′,l⊗l′,r⊗r′) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob, G::Ob, H::Ob,
     t::(A → B), b::(C → D), l::(A ↓ C), r::(B ↓ D),
     t′::(E → F), b′::(G → H), l′::(E ↓ G), r′::(F ↓ H))

  @op (⊗) := otimes
  munit()::Ob

  # Monoid axioms.
  (A ⊗ B) ⊗ C == A ⊗ (B ⊗ C) ⊣ (A::Ob, B::Ob, C::Ob)
  A ⊗ munit() == A ⊣ (A::Ob)
  munit() ⊗ A == A ⊣ (A::Ob)
  (f ⊗ g) ⊗ h == f ⊗ (g ⊗ h) ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
                                 f::(A → X), g::(B → Y), h::(C → Z))
  (f ⊗ g) ⊗ h == f ⊗ (g ⊗ h) ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
                                 f::(A ↓ X), g::(B ↓ Y), h::(C ↓ Z))
  (α ⊗ β) ⊗ γ == α ⊗ (β ⊗ γ) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                E::Ob, F::Ob, G::Ob, H::Ob,
                                I::Ob, J::Ob, K::Ob, L::Ob,
                                t1::(A → B), b1::(C → D), l1::(A ↓ C), r1::(B ↓ D),
                                t2::(E → F), b2::(G → H), l2::(E ↓ G), r2::(F ↓ H),
                                t3::(I → J), b3::(K → L), l3::(I ↓ K), r3::(J ↓ L),
                                α::Hom2(t1, b1, l1, r1),
                                β::Hom2(t2, b2, l2, r2),
                                γ::Hom2(t3, b3, l3, r3))

  # Functorality axioms.
  ((f ⊗ g) ⋅ (h ⊗ k) == (f ⋅ h) ⊗ (g ⋅ k)
    ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
       f::(A → B), h::(B → C), g::(X → Y), k::(Y → Z)))
  ((f ⊗ g) ⋅ (h ⊗ k) == (f ⋅ h) ⊗ (g ⋅ k)
    ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
       f::(A ↓ B), h::(B ↓ C), g::(X ↓ Y), k::(Y ↓ Z)))
  ((α ⊗ β) ⋅ (γ ⊗ δ) == (α ⋅ γ) ⊗ (β ⋅ δ)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
       E::Ob, F::Ob, G::Ob, H::Ob,
       I::Ob, J::Ob, K::Ob, L::Ob,
       M::Ob, N::Ob, O::Ob, P::Ob,
       t1::(A → B), b1::(C → D), l1::(A ↓ C), r1::(B ↓ D),
       t2::(E → F), b2::(G → H), l2::(E ↓ G), r2::(F ↓ H),
       t3::(I → J), b3::(K → L), l3::(I ↓ K), r3::(J ↓ L),
       t4::(M → N), b4::(O → P), l4::(M ↓ O), r4::(N ↓ P),
       α::Hom2(t1, b1, l1, r1),
       β::Hom2(t2, b2, l2, r2),
       γ::Hom2(t3, b3, l3, r3),
       δ::Hom2(t4, b4, l4, r4)))
  idH(A ⊗ B) == idH(A) ⊗ idH(B) ⊣ (A::Ob, B::Ob)
  idV(A ⊗ B) == idV(A) ⊗ idV(B) ⊣ (A::Ob, B::Ob)
  id2(A ⊗ B) == id2(A) ⊗ id2(B) ⊣ (A::Ob, B::Ob)
  id2H(f ⊗ g) == id2H(f) ⊗ id2H(g) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                      f::(A↓C), g::(B↓D))
  id2V(f ⊗ g) == id2V(f) ⊗ id2V(g) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                      f::(A→C), g::(B→D))
end

""" Theory of (strict) *symmetric monoidal double categories*
"""
@theory SymmetricMonoidalDoubleCategory{Ob,HomV,HomH,Hom2} <: MonoidalDoubleCategory{Ob,HomV,HomH,Hom2} begin
  braidH(A::Ob, B::Ob)::((A ⊗ B) → (B ⊗ A))
  braidV(A::Ob, B::Ob)::((A ⊗ B) ↓ (B ⊗ A))
  @op (σH) := braidH
  @op (σV) := braidV

  # Involutivity axiom.
  σH(A,B) ⋅ σH(B,A) == idH(A ⊗ B) ⊣ (A::Ob, B::Ob)
  σV(A,B) ⋅ σV(B,A) == idV(A ⊗ B) ⊣ (A::Ob, B::Ob)

  # Coherence axioms.
  σH(A,B⊗C) == (σH(A,B) ⊗ idH(C)) ⋅ (idH(B) ⊗ σH(A,C)) ⊣ (A::Ob, B::Ob, C::Ob)
  σV(A⊗B,C) == (idV(A) ⊗ σV(B,C)) ⋅ (σV(A,C) ⊗ idV(B)) ⊣ (A::Ob, B::Ob, C::Ob)

  # Naturality axiom.
  (f ⊗ g) ⋅ σH(B,D) == σH(A,C) ⋅ (g ⊗ f) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                           f::(A → B), g::(C → D))
  (f ⊗ g) ⋅ σV(B,D) == σV(A,C) ⋅ (g ⊗ f) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                           f::(A ↓ B), g::(C ↓ D))
end

@syntax FreeSymmetricMonoidalDoubleCategory{ObExpr,HomVExpr,HomHExpr,Hom2Expr} SymmetricMonoidalDoubleCategory begin
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::HomV, g::HomV) = associate(new(f,g))
  otimes(f::HomH, g::HomH) = associate(new(f,g))
  otimes(f::Hom2, g::Hom2) = associate(new(f,g))
  compose(f::HomV, g::HomV) = associate(new(f,g; strict=true))
  compose(f::HomH, g::HomH) = associate(new(f,g; strict=true))
  composeH(α::Hom2, β::Hom2) = associate(new(α,β))
  composeV(α::Hom2, β::Hom2) = associate(new(α,β))
end

function show_unicode(io::IO, expr::HomHExpr{:braidH}; kw...)
  Syntax.show_unicode_infix(io, expr, "σH"; kw...)
end
function show_unicode(io::IO, expr::HomVExpr{:braidV}; kw...)
  Syntax.show_unicode_infix(io, expr, "σV"; kw...)
end

function show_latex(io::IO, expr::HomHExpr{:braidH}; kw...)
  Syntax.show_latex_script(io, expr, "\\sigmaH")
end

function show_latex(io::IO, expr::HomVExpr{:braidV}; kw...)
  Syntax.show_latex_script(io, expr, "\\sigmaV")
end