text:
"""
    preduceBF(M, N; fast = true, atol = 0, rtol, withQ = true, withZ = true) -> (F, G, Q, Z, n, m, p)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'*(M - λN)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the pencil `M - λN` is transformed
into the following standard form

                   | B  | A-λE |
          F - λG = |----|------| ,
                   | D  |  C   |

where `E` is an `nxn` non-singular matrix, and `A`, `B`, `C`, `D` are `nxn`-, `nxm`-, `pxn`- and `pxm`-dimensional matrices,
respectively. The order `n` of `E` is equal to the numerical rank of `N` determined using the absolute tolerance `atol` and
relative tolerance `rtol`. `M` and `N` are overwritten by `F` and `G`, respectively.

If `fast = true`, `E` is determined upper triangular using a rank revealing QR-decomposition with column pivoting of `N`
and `n` is evaluated as the number of nonzero diagonal elements of the `R` factor, whose magnitudes are greater than
`tol = max(atol,abs(R[1,1])*rtol)`.
If `fast = false`,  `E` is determined diagonal using a rank revealing SVD-decomposition of `N` and
`n` is evaluated as the number of singular values greater than `tol = max(atol,smax*rtol)`, where `smax`
is the largest singular value.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function preduceBF(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true,
                   atol::Real = zero(real(eltype(M))),
                   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(atol),
                   withQ::Bool = true, withZ::Bool = true)
   # In interest of performance, no dimensional checks are performed
   mM, nM = size(M)
   (mM,nM) == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))
   T = promote_type(eltype(M), eltype(N))
   T <: BlasFloat || (T = promote_type(Float64,T))
   M1 = copy_oftype(M,T)
   N1 = copy_oftype(N,T)

   withQ ? (Q = Matrix{T}(I,mM,mM)) : (Q = nothing)
   withZ ? (Z = Matrix{T}(I,nM,nM)) : (Z = nothing)

   # Step 0: Reduce to the standard form
   n, m, p = _preduceBF!(M1, N1, Q, Z; atol = atol, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ)

   return M1, N1, Q, Z, n, m, p
end
"""
    klf_rlsplit(M, N; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (F, G, Q, Z, ν, μ, n, m, p)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'*(M - λN)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in one of the
following Kronecker-like forms:

(1) if `finite_infinite = false`, then

                   | Mri-λNri |     *      |
          F - λG = |----------|------------|,
                   |    O     | Mfl-λNfl   |

where the subpencil `Mri-λNri` contains the right Kronecker structure and infinite elementary
divisors and the subpencil `Mfl-λNfl` contains the finite and left Kronecker structure of the pencil `M-λN`.

The full row rank subpencil `Mri-λNri` is in a staircase form.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mri-λNri` such that `i`-th block has dimensions `ν[i] x μ[i]` and
has full row rank.
The difference `μ[i]-ν[i]` for `i = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size `(i-1) x i`.
The difference `ν[i]-μ[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `μ[nb+1] = 0`).

The full column rank subpencil `Mfl-λNfl` is in the form

                     | A-λE |
          Mfl-λNfl = |------| ,
                     |  C   |

where `E` is an nxn non-singular upper triangular matrix, and `A` and `C` are `nxn`- and `pxn`-dimensional matrices,
respectively, and `m = 0`.

(2) if `finite_infinite = true`, then

                   | Mrf-λNrf |     *      |
          F - λG = |----------|------------|,
                   |    O     | Mil-λNil   |

where the subpencil `Mrf-λNrf` contains the right Kronecker and finite Kronecker structures and
the subpencil `Mil-λNil` contains the left Kronecker structures and infinite elementary
divisors of the pencil `M-λN`.

The full row rank subpencil `Mrf-λNrf` is in the form

          Mrf-λNrf = | B  | A-λE | ,

where `E` is an `nxn` non-singular upper triangular matrix, and `A` and `B` are `nxn`- and `nxm`-dimensional matrices,
respectively, and `p = 0`.

The full column rank sub pencil `Mil-λNil` is in a staircase form.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mil-λNil` such that the `i`-th block has dimensions `ν[i] x μ[i]` and has full column rank.
The difference `ν[nb-j+1]-μ[nb-j+1]` for `j = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size
`j x (j-1)`. The difference `μ[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function klf_rlsplit(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, finite_infinite::Bool = false,
                     atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                     rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)),
                     withQ::Bool = true, withZ::Bool = true)

   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast)

   mM, nM = size(M)
   maxmn = max(mM,nM)
   μ = Vector{Int}(undef,maxmn)
   ν = Vector{Int}(undef,maxmn)

   # fast returns for null dimensions
   if mM == 0 && nM == 0
      return M1, N1, Q, Z, ν, μ, n, m, p
   end

   if finite_infinite
      if mM == 0
         return M1, N1, Q, Z, ν[1:0], μ[1:0], n, m, p
      elseif nM == 0
         ν[1] = mM
         μ[1] = 0
         p = 0
         return M1, N1, Q, Z, ν[1:1], μ[1:1], n, m, p
      end

      # Reduce M-λN to a KLF which exhibits the splitting of the right-finite and infinite-left structures
      #
      #                  | Mrf - λ Nrf |     *        |
      #      M1 - λ N1 = |-------------|--------------|
      #                  |    0        | Mil -  λ Nil |
      #
      # where Mil -  λ Nil is in a staircase form.

      mrinf = 0
      nrinf = 0
      rtrail = 0
      ctrail = 0
      i = 0
      tol1 = max(atol1, rtol*opnorm(M,1))
      while p > 0
         # Step 1 & 2: Dual algorithm PREDUCE
         τ, ρ  = _preduce2!(n, m, p, M1, N1, Q, Z, tol1; fast = fast,
                            roff = mrinf, coff = nrinf, rtrail = rtrail, ctrail = ctrail,
                            withQ = withQ, withZ = withZ)
         i += 1
         ν[i] = p
         μ[i] = ρ+τ
         ctrail += ρ+τ
         rtrail += p
         n -= ρ
         p = ρ
         m -= τ
      end
      return M1, N1, Q, Z, reverse(ν[1:i]), reverse(μ[1:i]), n, m, p
   else
      if mM == 0
         ν[1] = 0
         μ[1] = nM
         m = 0
         return M1, N1, Q, Z, ν[1:1], μ[1:1], n, m, p
      elseif nM == 0
         return M1, N1, Q, Z, ν[1:0], μ[1:0], n, m, p
      end

      # Reduce M-λN to a KLF which exhibits the splitting the right-infinite and finite-left structures
      #
      #                  | Mri - λ Nri |     *        |
      #      M1 - λ N1 = |-------------|--------------|
      #                  |    0        | Mfl -  λ Nfl |
      # where Mri - λ Nri is in a staircase form.

      mrinf = 0
      nrinf = 0
      i = 0
      tol1 = max(atol1, rtol*opnorm(M,1))

      while m > 0
         # Steps 1 & 2: Standard algorithm PREDUCE
          τ, ρ = _preduce1!( n, m, p, M1, N1, Q, Z, tol1; fast = fast,
                            roff = mrinf, coff = nrinf, withQ = withQ, withZ = withZ)
         i += 1
         ν[i] = ρ+τ
         μ[i] = m
         mrinf += ρ+τ
         nrinf += m
         n -= ρ
         m = ρ
         p -= τ
      end
      return M1, N1, Q, Z, ν[1:i], μ[1:i], n, m, p
   end
end
"""
    klf(M, N; fast = true, finite_infinite = false, ut = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (F, G, Q, Z, νr, μr, νi, nf, νl, μl)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'(M - λN)Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in the
following Kronecker-like form exhibiting the complete Kronecker structure:

                   |  Mr-λNr  |     *      |    *    |
                   |----------|------------|---------|
         F - λG =  |    O     | Mreg-λNreg |    *    |
                   |----------|------------|---------|
                   |    O     |     0      |  Ml-λNl |

The full row rank pencil `Mr-λNr` is in a staircase form, contains the right Kronecker indices
of the pencil `M-λN`and has the form

         Mr-λNr  = | Br | Ar-λEr |,

where `Er` is upper triangular and nonsingular.
The `nr`-dimensional vectors `νr` and `μr` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mr-λNr` such that the `i`-th block has dimensions `νr[i] x μr[i]` and
has full row rank. The difference `μr[i]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker
blocks of size `(i-1) x i`.
If `ut = true`, the full row rank diagonal blocks of `Mr` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nr` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.

The pencil `Mreg-λNreg` is regular, contains the infinite and finite elementary divisors of `M-λN` and has the form

                       | Mi-λNi |    *    |
         Mreg-λNreg  = |--------|---------|, if `finite_infinite = false`, or the form
                       |   0    |  Mf-λNf |

                       | Mf-λNf |    *    |
         Mreg-λNreg  = |--------|---------|, if `finite_infinite = true`,
                       |   0    |  Mi-λNi |

where: (1) `Mi-λNi`, in staircase form, contains the infinite elementary divisors of `M-λN`,
`Mi` upper triangular if `ut = true` and nonsingular, and `Ni` is upper triangular and nilpotent;
(2) `Mf-λNf` contains the infinite elementary divisors of `M-λN` and
`Nf` is upper triangular and nonsingular.
The `ni`-dimensional vector `νi` contains the dimensions of the square blocks of the staircase form  `Mi-λNi`
such that the `i`-th block has dimensions `νi[i] x νi[i]`.
If `finite_infinite = true`, the difference `νi[i]-νi[i+1]` for `i = 1, 2, ..., ni`
is the number of infinite elementary divisors of degree `i` (with `νi[ni] = 0`).
If `finite_infinite = false`, the difference `νi[ni-i+1]-νi[ni-i]` for `i = 1, 2, ..., ni`
is the number of infinite elementary divisors of degree `i` (with `νi[0] = 0`).

The full column rank pencil `Ml-λNl`, in a staircase form, contains the left Kronecker indices of `M-λN` and has the form

                   | Al-λEl |
         Ml-λNl  = |--------|,
                   |   Cl   |

where `El` is upper triangular and nonsingular.
The `nl`-dimensional vectors `νl` and `μl` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Ml-λNl` such that the `j`-th block has dimensions `νl[j] x μl[j]` and has full column rank.
The difference `νl[nl-j+1]-μl[nl-j+1]` for `j = 1, 2, ..., nl` is the number of elementary Kronecker blocks of size
`j x (j-1)`.
If `ut = true`, the full column rank diagonal blocks of `Ml` are reduced to the form `[X; 0]`
with `X` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nl` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function klf(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, finite_infinite::Bool = false, ut::Bool = false,
             atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
             rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(min(atol1,atol2)), withQ::Bool = true, withZ::Bool = true)

   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast)
   mM, nM = size(M)

   if finite_infinite

      # Reduce M-λN to a KLF exhibiting the right and finite structures
      #                  [ Mr - λ Nr  |   *        |     *        ]
      #      M1 - λ N1 = [    0       | Mf -  λ Nf |     *        ]
      #                  [    0       |    0       | Mil -  λ Nil ]

      νr, μr, nf, ν, μ, tol1 = klf_right!(n, m, p, M1, N1, Q, Z, atol = atol1, rtol = rtol,
                                          withQ = withQ, withZ = withZ, fast = fast)
      if mM == 0 || nM == 0
          return  M1, N1, Q, Z, νr, μr, ν[1:0], nf, ν, μ
      end

      mr = sum(νr)+nf
      nr = sum(μr)+nf
      ut && klf_right_refineut!(νr, μr, M1, N1, Q, Z, ctrail = nM-nr, withQ = withQ, withZ = withZ)

      # Reduce Mil-λNil to a KLF exhibiting the infinite and left structures and update M1 - λ N1 to
      #                  [ Mr - λ Nr  |   *        |     *      |    *       ]
      #     M2 -  λ N2 = [    0       | Mf -  λ Nf |     *      |    *       ]
      #                  [    0       |    0       | Mi -  λ Ni |    *       ]
      #                  [    0       |    0       |    0       | Ml -  λ Nl ]

      jM2 = nr+1:nr+sum(μ)
      M2 = view(M1,:,jM2)
      N2 = view(N1,:,jM2)
      withZ ? (Z2 = view(Z,:,jM2)) : (Z2 = nothing)
      νi, νl, μl = klf_left_refine!(ν, μ, M2, N2, Q, Z2, tol1, roff = mr,
                                    withQ = withQ, withZ = withZ, fast = fast, ut = ut)
   else

      # Reduce M-λN to a KLF exhibiting the left and finite structures
      #                  [ Mri - λ Nri  |   *        |     *      ]
      #      M1 - λ N1 = [     0        | Mf -  λ Nf |     *      ]
      #                  [     0        |    0       | Ml -  λ Nl ]

      ν, μ, nf, νl, μl, tol1 = klf_left!(n, m, p, M1, N1, Q, Z, atol = atol1, rtol = rtol,
                                         withQ = withQ, withZ = withZ, fast = fast)
      if mM == 0 || nM == 0
         return  M1, N1, Q, Z, ν, μ, ν[1:0], nf, νl, μl
      end

      ut && klf_left_refineut!(νl, μl, M1, N1, Q, Z, roff = mM-sum(νl), coff = nM-sum(μl), withQ = withQ, withZ = withZ)

      # Reduce Mri-λNri to a KLF exhibiting the right and infinite structures and update M1 - λ N1 to
      #                  [ Mr - λ Nr  |   *        |     *      |    *       ]
      #     M2 -  λ N2 = [    0       | Mi -  λ Ni |     *      |    *       ]
      #                  [    0       |    0       | Mf -  λ Nf |    *       ]
      #                  [    0       |    0       |    0       | Ml -  λ Nl ]

      iM11 = 1:sum(ν)
      M11 = view(M1,iM11,:)
      N11 = view(N1,iM11,:)
      νr, μr, νi = klf_right_refine!(ν, μ, M11, N11, Q, Z, tol1, ctrail = nM-sum(μ),
                                     withQ = withQ, withZ = withZ, fast = fast, ut = ut)
   end
   return  M1, N1, Q, Z, νr, μr, νi, nf, νl, μl
end
"""
    klf_right(M, N; fast = true, ut = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (F, G, Q, Z, νr, μr, nf, ν, μ)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'(M - λN)Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in the
following Kronecker-like form exhibiting the right and finite Kronecker structures:

                   |  Mr-λNr    |     *      |    *     |
                   |------------|------------|----------|
         F - λG =  |    O       |   Mf-λNf   |    *     |
                   |------------|------------|----------|
                   |    O       |     0      | Mil-λNil |

The full row rank pencil `Mr-λNr`, in a staircase form, contains the right Kronecker indices of `M-λN` and has the form

         Mr-λNr  = | Br | Ar-λEr |,

where `Er` is upper triangular and nonsingular.
The `nr`-dimensional vectors `νr` and `μr` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mr-λNr` such that `i`-th block has dimensions `νr[i] x μr[i]` and
has full row rank. The difference `μr[i]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks
of size `(i-1) x i`.
If `ut = true`, the full row rank diagonal blocks of `Mr` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nr` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.

The `nf x nf` pencil `Mf-λNf` is regular and contains the finite elementary divisors of `M-λN`.
`Nf` is upper triangular and nonsingular.

The full column rank pencil `Mil-λNil` is in a staircase form, and contains the  left Kronecker indices
and infinite elementary divisors of the pencil `M-λN`.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mil-λNil` such that the `i`-th block has dimensions `ν[i] x μ[i]` and has full column rank.
The difference `ν[nb-j+1]-μ[nb-j+1]` for `j = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size
`j x (j-1)`. The difference `μ[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).
If `ut = true`, the full column rank diagonal blocks of `Mil` are reduced to the form `[X; 0]`
with `X` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nil` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.

`Note:` If the pencil `M - λN` has full row rank, then the regular pencil `Mil-λNil` is in a staircase form with
square upper triangular diagonal blocks (i.e.,`μ[j] = ν[j]`), and the difference `ν[nb-j+1]-ν[nb-j]` for
`j = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `j` (with `ν[0] = 0`).
"""
function klf_right(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, ut::Bool = false,
                   atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(max(atol1,atol2)),
                   withQ::Bool = true, withZ::Bool = true)

   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast)

   νr, μr, nf, ν, μ, tol1 = klf_right!(n, m, p, M1, N1, Q, Z, atol = atol1, rtol = rtol,
                                       withQ = withQ, withZ = withZ, fast = fast)

   if ut
      mM, nM = size(M)
      mr = sum(νr)+nf
      nr = sum(μr)+nf
      klf_right_refineut!(νr, μr, M1, N1, Q, Z, ctrail = nM-nr, withQ = withQ, withZ = withZ)
      klf_left_refineut!(ν, μ, M1, N1, Q, Z, roff = mM-sum(ν), coff = nM-sum(μ), withQ = withQ, withZ = withZ)
   end
   return M1, N1, Q, Z, νr, μr, nf, ν, μ

end
"""
    klf_left(M, N; fast = true, ut = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (F, G, Q, Z, ν, μ, nf, νl, μl)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'(M - λN)Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in the
following Kronecker-like form exhibiting the finite and left Kronecker structures:

                   |  Mri-λNri  |     *      |    *    |
                   |------------|------------|---------|
         F - λG =  |    O       |   Mf-λNf   |    *    |
                   |------------|------------|---------|
                   |    O       |     0      |  Ml-λNl |

The full row rank pencil `Mri-λNri` is in a staircase form, and contains the right Kronecker indices
and infinite elementary divisors of the pencil `M-λN`.

The `nf x nf` pencil `Mf-λNf` is regular and contains the finite elementary divisors of `M-λN`.
`Nf` is upper triangular and nonsingular.

The full column rank pencil `Ml-λNl`, in a staircase form, contains the left Kronecker indices of `M-λN` and has the form

                   | Al-λEl |
         Ml-λNl  = |--------|,
                   |   Cl   |

where `El` is upper triangular and nonsingular.

The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mri-λNri` such that `i`-th block has dimensions `ν[i] x μ[i]` and
has full row rank.
The difference `μ[i]-ν[i]` for `i = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size `(i-1) x i`.
The difference `ν[i]-μ[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `μ[nb+1] = 0`).
If `ut = true`, the full row rank diagonal blocks of `Mri` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nri` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.

The `nl`-dimensional vectors `νl` and `μl` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Ml-λNl` such that the `j`-th block has dimensions `νl[j] x μl[j]` and has full column rank.
The difference `νl[nl-j+1]-μl[nl-j+1]` for `j = 1, 2, ..., nl` is the number of elementary Kronecker blocks of size
`j x (j-1)`.
If `ut = true`, the full column rank diagonal blocks of `Ml` are reduced to the form `[X; 0]`
with `X` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nl` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.

`Note:` If the pencil `M - λN` has full column rank, then the regular pencil `Mri-λNri` is in a staircase form with
square upper triangular diagonal blocks (i.e.,`μ[i] = ν[i]`), and the difference `ν[i+1]-ν[i]` for `i = 1, 2, ..., nb`
is the number of infinite elementary divisors of degree `i` (with `ν[nb+1] = 0`).
"""
function klf_left(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, ut::Bool = false,
                  atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                  rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(max(atol1,atol2)),
                  withQ::Bool = true, withZ::Bool = true)
   # Step 0: Reduce to the standard form
   M1, N1, Q, Z, n, m, p = preduceBF(M, N; atol = atol2, rtol = rtol, fast = fast)

   ν, μ, nf, νl, μl, tol1 = klf_left!(n, m, p, M1, N1, Q, Z, atol = atol1, rtol = rtol,
                                      withQ = withQ, withZ = withZ, fast = fast)
   if ut
      mM, nM = size(M)
      klf_left_refineut!(νl, μl, M1, N1, Q, Z, roff = mM-sum(νl), coff = nM-sum(μl), withQ = withQ, withZ = withZ)
      klf_right_refineut!(ν, μ, M1, N1, Q, Z, ctrail = nM-sum(μ), withQ = withQ, withZ = withZ)
   end
   return M1, N1, Q, Z, ν, μ, nf, νl, μl
end
"""
    klf_rightinf(M, N; fast = true, ut = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true)
                -> (F, G, Q, Z, νr, μr, νi, n, p)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'(M - λN)Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in the
following Kronecker-like form exhibiting the infinite and left Kronecker structures:

                   |  Mr-λNr  |     *      |    *     |
                   |----------|------------|----------|
         F - λG =  |    O     |   Mi-λNi   |    *     |
                   |----------|------------|----------|
                   |    O     |     0      | Mfl-λNfl |

The full row rank pencil `Mr-λNr` is in a staircase form, contains the right Kronecker indices
of the pencil `M-λN`and has the form

         Mr-λNr  = | Br | Ar-λEr |,

where `Er` is upper triangular and nonsingular.
The `nr`-dimensional vectors `νr` and `μr` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mr-λNr` such that the `i`-th block has dimensions `νr[i] x μr[i]` and
has full row rank. The difference `μr[i]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker
blocks of size `(i-1) x i`.
If `ut = true`, the full row rank diagonal blocks of `Mr` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nr` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.

The regular pencil `Mi-λNi` is in a staircase form, contains the infinite elementary divisors of `M-λN`
with `Mi` upper triangular if `ut = true` and nonsingular, and `Ni` is upper triangular and nilpotent.
The `ni`-dimensional vector `νi` contains the dimensions of the
square diagonal blocks of the staircase form  `Mi-λNi` such that the `i`-th block has dimensions `νi[i] x νi[i]`.
The difference `νi[ni-i+1]-νi[ni-i]` for `i = 1, 2, ..., ni` is the number of infinite elementary divisors of degree `i` (with `νi[0] = 0`).

The full column rank pencil `Mfl-λNfl` contains the left Kronecker
and finite Kronecker structures of the pencil `M-λN` and is in the form

                     | A-λE |
          Mfl-λNfl = |------| ,
                     |  C   |

where `E` is an nxn non-singular upper triangular matrix, and `A` and `C` are `nxn`- and `pxn`-dimensional matrices,
respectively.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function klf_rightinf(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, ut::Bool = false,
                  atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                  rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(max(atol1,atol2)),
                  withQ::Bool = true, withZ::Bool = true)
   M1, N1, Q, Z, ν, μ, n, m, p = klf_rlsplit(M, N; fast = fast, finite_infinite = false,
                                             atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
   tol1 = max(atol1, rtol*opnorm(M,1))
   νr, μr, νi = klf_right_refine!(ν, μ, M1, N1, Q, Z, tol1, ut = ut, withQ = withQ, withZ = withZ, fast = fast)
   return M1, N1, Q, Z, νr, μr, νi, n, p
end
"""
    klf_leftinf(M, N; fast = true, ut = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true)
                -> (F, G, Q, Z, n, m, νi, νl, μl)

Reduce the matrix pencil `M - λN` to an equivalent form `F - λG = Q'(M - λN)Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `F` and `G` are in the
following Kronecker-like form exhibiting the infinite and left Kronecker structures:

                   |  Mrf-λNrf  |     *      |    *    |
                   |------------|------------|---------|
         F - λG =  |    O       |   Mi-λNi   |    *    |
                   |------------|------------|---------|
                   |    O       |     0      |  Ml-λNl |

The full row rank pencil `Mrf-λNrf` contains the right Kronecker
and finite Kronecker structures of the pencil `M-λN` and is in the standard form

          Mrf-λNrf = | B  | A-λE | ,

where `E` is an `nxn` non-singular upper triangular matrix, and `A` and `B` are `nxn`- and `nxm`-dimensional matrices,
respectively.

The regular pencil `Mi-λNi` is in a staircase form, contains the infinite elementary divisors of `M-λN`
with `Mi` upper triangular if `ut = true` and nonsingular, and `Ni` is upper triangular and nilpotent.
The `ni`-dimensional vector `νi` contains the dimensions of the
square diagonal blocks of the staircase form  `Mi-λNi` such that the `i`-th block has dimensions `νi[i] x νi[i]`.
The difference `νi[i]-νi[i-1] = ν[ni-i+1]-μ[ni-i+1]` for `i = 1, 2, ..., ni` is the number of infinite elementary
divisors of degree `i` (with `νi[0] = 0` and `μ[nb+1] = 0`).

The full column rank pencil `Ml-λNl`, in a staircase form, contains the left Kronecker indices of `M-λN` and has the form

                   | Al-λEl |
         Ml-λNl  = |--------|,
                   |   Cl   |

where `El` is upper triangular and nonsingular.

The `nl`-dimensional vectors `νl` and `μl` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Ml-λNl` such that the `j`-th block has dimensions `νl[j] x μl[j]` and has full column rank.
The difference `νl[nl-j+1]-μl[nl-j+1]` for `j = 1, 2, ..., nl` is the number of elementary Kronecker blocks of size
`j x (j-1)`.
If `ut = true`, the full column rank diagonal blocks of `Ml` are reduced to the form `[X; 0]`
with `X` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nl` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function klf_leftinf(M::AbstractMatrix, N::AbstractMatrix; fast::Bool = true, ut::Bool = false,
                  atol1::Real = zero(real(eltype(M))), atol2::Real = zero(real(eltype(M))),
                  rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(max(atol1,atol2)),
                  withQ::Bool = true, withZ::Bool = true)
   M1, N1, Q, Z, ν, μ, n, m, p = klf_rlsplit(M, N; fast = fast, finite_infinite = true,
            atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
   mr = n+p
   nr = n+m
   jM2 = nr+1:nr+sum(μ)
   M2 = view(M1,:,jM2)
   N2 = view(N1,:,jM2)
   withZ ? (Z2 = view(Z,:,jM2)) : (Z2 = nothing)
   tol1 = max(atol1, rtol*opnorm(M2,1))
   νi, νl, μl = klf_left_refine!(ν, μ, M2, N2, Q, Z2, tol1, ut = ut, roff = mr,
                                 withQ = withQ, withZ = withZ, fast = fast)
   return M1, N1, Q, Z, n, m, νi, νl, μl
end
"""
    klf_right!(M, N; fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, atol = 0, rtol, withQ = true, withZ = true) -> (Q, Z, νr, μr, nf, ν, μ, tol)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

              [   *         *        *  ] roff
     M - λN = [   0      M22-λN22    *  ] m
              [   0         0        *  ] rtrail
                coff        n     ctrail

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation matrices `Q1` and `Z1`
such that the subpencil `M22 - λN22` is transformed into the following Kronecker-like form exhibiting
its right and finite Kronecker structures:

                   |  Mr-λNr    |     *      |    *     |
                   |------------|------------|----------|
     F22 - λG22 =  |    O       |   Mf-λNf   |    *     |
                   |------------|------------|----------|
                   |    O       |     0      | Mil-λNil |

`F` and `G` are returned in `M` and `N`, respectively.

The full row rank pencil `Mr-λNr`, in a staircase form, contains the right Kronecker indices of
the subpencil `M22 - λN22` and has the form

         Mr-λNr  = | Br | Ar-λEr |,

where `Er` is upper triangular and nonsingular.
The `nr`-dimensional vectors `νr` and `μr` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mr-λNr` such that `i`-th block has dimensions `νr[i] x μr[i]` and
has full row rank. The difference `μr[i]-νr[i]` for `i = 1, 2, ..., nr` is the number of elementary Kronecker blocks
of size `(i-1) x i`.

The `nf x nf` pencil `Mf-λNf` is regular and contains the finite elementary divisors of the subpencil `M22 - λN22`.
`Nf` is upper triangular and nonsingular.

The full column rank pencil `Mil-λNil` is in a staircase form, and contains the  left Kronecker indices
and infinite elementary divisors of the subpencil `M22 - λN22`.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mil-λNil` such that the `i`-th block has dimensions `ν[nb-i+1] x μ[nb-i+1]` and has full column rank.
The difference `ν[nb-i+1]-μ[nb-i+1]` for `i = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size
`i x (i-1)`. The difference `μ[nb-i+1]-ν[nb-i]` for `i = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `i` (with `ν[0] = 0`).

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero elements of `M`, respectively. The internally employed absolute tolerance  for the
nonzero elements of `M` is returned in `tol`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` (i.e., `Q <- Q*Q1`)
if `withQ = true`. Otherwise, `Q` is not modified.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` (i.e., `Z <- Z*Z1`)
if `withZ = true`. Otherwise, `Z` is not modified.

`Note:` If the subpencil `M22 - λN22` has full row rank, then the regular pencil `Mil-λNil` is in a staircase form with
square upper triangular diagonal blocks (i.e.,`μ[i] = ν[i]`), and the difference `ν[nb-i+1]-ν[nb-i]` for
`i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i` (with `ν[0] = 0`).
"""
function klf_right!(n::Int, m::Int, p::Int, M::AbstractMatrix{T1}, N::AbstractMatrix{T1},
                    Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing};
                    fast::Bool = true, atol::Real = zero(real(eltype(M))),
                    rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(atol),
                    roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat
   mM, nM = size(M)
   (mM,nM) == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))
   (!isa(M,Adjoint) && !isa(N,Adjoint)) || error("No adjoint inputs are supported")

   # Step 0: Reduce M22-λN22 to the standard form
   #n, m, p = _preduceBF!(M, N, Q, Z; atol = atol2, rtol = rtol, fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail, withQ = withQ, withZ = withQ)

   maxmn = max(mM,nM)
   μr = Vector{Int}(undef,maxmn)
   νr = Vector{Int}(undef,maxmn)
   μ = Vector{Int}(undef,maxmn)
   ν = Vector{Int}(undef,maxmn)
   nf = 0
   tol = atol

   # fast returns for null dimensions
   if mM == 0 && nM == 0
      return νr, μr, nf, ν, μ, tol
   elseif mM == 0
      νr[1] = 0
      μr[1] = nM
      return νr[1:1], μr[1:1], nf, ν[1:0], μ[1:0], tol
   elseif nM == 0
      ν[1] = mM
      μ[1] = 0
      return νr[1:0], μr[1:0], nf, ν[1:1], μ[1:1], tol
   end

   j = 0
   tol = max(atol, rtol*opnorm(M,1))

   while p > 0
      # Steps 1 & 2: Dual algorithm PREDUCE
      τ, ρ = _preduce2!(n, m, p, M, N, Q, Z, tol;
                        fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                        withQ = withQ, withZ = withZ)
      j += 1
      ν[j] = p
      μ[j] = ρ+τ
      ctrail += τ+ρ
      rtrail += p
      n -= ρ
      p = ρ
      m -= τ
   end
   i = 0
   if m > 0
      imM11 = 1:mM-rtrail
      M11 = view(M,imM11,1:nM)
      N11 = view(N,imM11,1:nM)
   end
   while m > 0
      # Step 3: Particular case of the standard algorithm PREDUCE
      ρ = _preduce3!(n, m, M11, N11, Q, Z, tol,
                     fast = fast, coff = coff, roff = roff, ctrail = ctrail,  withQ = withQ, withZ = withZ)
      i += 1
      νr[i] = ρ
      μr[i] = m
      roff += ρ
      coff += m
      n -= ρ
      m = ρ
   end

   return νr[1:i], μr[1:i], n, reverse(ν[1:j]), reverse(μ[1:j]), tol
end
"""

     klf_right_refine!(ν, μ, M, N, tol; fast = true, ut = false, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) -> (νr, μr, νi)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [   *         *        *  ] roff
    M - λN = [   0      Mri-λNri    *  ] mri
             [   0         0        *  ] rtrail
               coff      nri     ctrail

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation matrices `Q1` and `Z1`
such that the full row rank subpencil `Mri-λNri`is transformed into the
following Kronecker-like form  exhibiting its right and infinite Kronecker structures:

                  |  Mr-λNr    |     *      |
     Fri - λGri = |------------|------------|
                  |    O       |   Mi-λNi   |

The full row rank pencil `Mri-λNri` is in a staircase form and the `nb`-dimensional vectors `ν` and `μ`
contain the row and, respectively, column dimensions of the blocks of the staircase form  `Mri-λNri` such that
the `i`-th block has dimensions `ν[i] x μ[i]` and has full row rank. The matrix Mri has full row rank and
the trailing `μ[2]+...+μ[nb]` columns of `Nri` form a full column rank submatrix.

The full row rank pencil `Mr-λNr` is in a staircase form, contains the right Kronecker indices
of the pencil `M-λN`and has the form

      Mr-λNr  = | Br | Ar-λEr |,

where `Er` is upper triangular and nonsingular.
The `nr`-dimensional vectors `νr` and `μr` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mr-λNr` such that the `i`-th block has dimensions `νr[i] x μr[i]` and
has full row rank. The difference `μr[i]-νr[i] = μ[i]-ν[i]` for `i = 1, 2, ..., nr` is the number of
elementary Kronecker blocks of size `(i-1) x i`.
If `ut = true`, the full row rank diagonal blocks of `Mr` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nr` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.

The regular pencil `Mi-λNi`is in a staircase form, contains the infinite elementary divisors of `Mri-λNri`,
`Mi` is upper triangular if `ut = true` and nonsingular and `Ni` is upper triangular and nilpotent. The `ni`-dimensional vector `νi` contains the dimensions of the
square diagonal blocks of the staircase form  `Mi-λNi` such that the `i`-th block has dimensions `νi[i] x νi[i]`.
The difference `νi[ni-i+1]-νi[ni-i] = ν[i]-μ[i+1]` for `i = 1, 2, ..., ni` is the number of infinite elementary
divisors of degree `i` (with `νi[0] = 0` and `μ[nb+1] = 0`).

`F` and `G` are returned in `M` and `N`, respectively.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` (i.e., `Q <- Q*Q1`)
if `withQ = true`. Otherwise, `Q` is not modified.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` (i.e., `Z <- Z*Z1`)
if `withZ = true`. Otherwise, `Z` is not modified.
"""
function klf_right_refine!(ν::Vector{Int}, μ::Vector{Int}, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, Q::Union{AbstractMatrix{T1},Nothing},
   Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; fast::Bool = true, ut::Bool = false, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
   withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat
   nb = length(ν)
   nb == length(μ) || throw(DimensionMismatch("ν and μ must have the same lengths"))

   nb == 0 && (return ν[1:0], ν[1:0], ν[1:0])

   mri = sum(ν)
   nri = sum(μ)
   mM = roff + mri + rtrail
   nM = coff + nri + ctrail

   rtrail0 = rtrail
   ctrail0 = ctrail

   m = μ[1]
   n = nri-m
   p = mri-n

   if n > 0
      # Step 0: Reduce Nri = [ 0 E11] to standard form, where E11 is full column rank
      #                      [ 0 0  ]
      it = roff+1:roff+mri
      jt = coff+m+1:coff+nri
      tau = similar(N,n)
      E11 = view(N,it,jt)
      LinearAlgebra.LAPACK.geqrf!(E11,tau)
      eltype(M) <: Complex ? tran = 'C' : tran = 'T'
      LinearAlgebra.LAPACK.ormqr!('L',tran,E11,tau,view(M,it,coff+1:nM))
      withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E11,tau,view(Q,:,it))
      LinearAlgebra.LAPACK.ormqr!('L',tran,E11,tau,view(N,it,coff+nri+1:nM))
      triu!(E11)
   end

   μr = Vector{Int}(undef,nb)
   νr = Vector{Int}(undef,nb)
   νi = Vector{Int}(undef,nb)

   mrinf = 0
   nrinf = 0
   j = 0
   ni = 0

   while p > 0
      # Steps 1 & 2: Dual algorithm PREDUCE to separate right-finite and infinite-left parts.

      # Reduce Mri-λNri to [ Mr1-λNr1 * ; 0 Mi-λNi], where Mr1-λNr1 contains the right Kronecker structure and the empty
      # finite part and Mi-λNi contains the infinite elementary divisors and the empty left Kronecker structure.
      τ, ρ = _preduce2!(n, m, p, M, N, Q, Z, tol;
                        fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                        withQ = withQ, withZ = withZ)
      ρ+τ == p || error("The reduced pencil must not have left structure: try to adjust the tolerances")
      j += 1
      νi[j] = p
      ni += p
      ctrail += p
      rtrail += p
      n -= ρ
      p = ρ
      m -= τ
   end

   # make Mi upper triangular
   reverse!(view(νi,1:j))
   ut && ni > 0 && klf_right_refineinf!(view(νi,1:j), M, N, Z, missing; roff = mM-rtrail0-ni, coff = nM-ctrail0-ni, withZ = withZ)

   i = 0
   if m > 0
      imM11 = 1:mM-rtrail
      M11 = view(M,imM11,1:nM)
      N11 = view(N,imM11,1:nM)
   end
   while m > 0
      # Step 3: Particular form of the standard algorithm PREDUCE to reduce Mr1-λNr1 to Mr-λNr in  staircase form.
      ρ = _preduce3!(n, m, M11, N11, Q, Z, tol, fast = fast,
                     coff = coff, roff = roff, ctrail = ctrail, withQ = withQ, withZ = withZ)
      i += 1
      νr[i] = ρ
      μr[i] = m
      roff += ρ
      coff += m
      n -= ρ
      m = ρ
   end
   ut && i > 0 && klf_right_refineut!(view(νr,1:i), view(νr,1:i), M, N, Q, Z;
                                      ctrail = nM-sum(view(μr,1:i)), withQ = withQ, withZ = withZ)

   return νr[1:i], μr[1:i], νi[1:j]
end
"""
    klf_left!(M, N; fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, atol = 0, rtol, withQ = true, withZ = true) -> (Q, Z, ν, μ, nf, νl, μl, tol)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

              [   *         *        *  ] roff
     M - λN = [   0      M22-λN22    *  ] m
              [   0         0        *  ] rtrail
                coff        n     ctrail

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation matrices `Q1` and `Z1`
such that the subpencil `M22 - λN22`  is transformed into the following Kronecker-like form exhibiting
its finite and left Kronecker structures

                   |  Mri-λNri  |     *      |    *    |
                   |------------|------------|---------|
     F22 - λG22 =  |    O       |   Mf-λNf   |    *    |
                   |------------|------------|---------|
                   |    O       |     0      |  Ml-λNl |

`F` and `G` are returned in `M` and `N`, respectively.

The full row rank pencil `Mri-λNri` is in a staircase form, and contains the right Kronecker indices
and infinite elementary divisors of the subpencil `M22 - λN22`.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mri-λNri` such that `i`-th block has dimensions `ν[i] x μ[i]` and
has full row rank.
The difference `μ[i]-ν[i]` for `i = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size `(i-1) x i`.
The difference `ν[i]-μ[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `μ[nb+1] = 0`).

The `nf x nf` pencil `Mf-λNf` is regular and contains the finite elementary divisors of `M-λN`.
`Nf` is upper triangular and nonsingular.

The full column rank pencil `Ml-λNl`, in a staircase form, contains the left Kronecker indices of `M-λN` and has the form

                   | Al-λEl |
         Ml-λNl  = |--------|,
                   |   Cl   |

where `El` is upper triangular and nonsingular.
The `nl`-dimensional vectors `νl` and `μl` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Ml-λNl` such that `j`-th block has dimensions `νl[nl-j+1] x μl[nl-j+1]` and has full column rank.
The difference `νl[nl-j+1]-μl[nl-j+1]` for `j = 1, 2, ..., nl` is the number of elementary Kronecker blocks of size
`j x (j-1)`.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero elements of `M`, respectively. The internally employed absolute tolerance  for the
nonzero elements of `M` is returned in `tol`.
The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` (i.e., `Q <- Q*Q1`)
if `withQ = true`. Otherwise, `Q` is not modified.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` (i.e., `Z <- Z*Z1`)
if `withZ = true`. Otherwise, `Z` is not modified.

`Note:` If the pencil `M22 - λN22` has full column rank, then the regular pencil `Mri-λNri` is in a staircase form with
square diagonal blocks (i.e.,`μ[i] = ν[i]`), and the difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb`
is the number of infinite elementary divisors of degree `i` (with `ν[nb+1] = 0`).
"""
function klf_left!(n::Int, m::Int, p::Int, M::AbstractMatrix{T1}, N::AbstractMatrix{T1},
                   Q::Union{AbstractMatrix{T1},Nothing}, Z::Union{AbstractMatrix{T1},Nothing};
                   fast::Bool = true, atol::Real = zero(real(eltype(M))),
                   rtol::Real = (min(size(M)...)*eps(real(float(one(eltype(M))))))*iszero(atol),
                   roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0, withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat
   mM, nM = size(M)
   (mM,nM) == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))
   (!isa(M,Adjoint) && !isa(N,Adjoint)) || error("No adjoint inputs are supported")

   maxmn = max(mM,nM)
   μ = Vector{Int}(undef,maxmn)
   ν = Vector{Int}(undef,maxmn)
   μl = Vector{Int}(undef,maxmn)
   νl = Vector{Int}(undef,maxmn)
   nf = 0
   tol = atol

   # fast returns for null dimensions
   if mM == 0 && nM == 0
      return ν, μ, nf, νl, μl, tol
   elseif mM == 0
      ν[1] = 0
      μ[1] = nM
      return ν[1:1], μ[1:1], nf, νl[1:0], μl[1:0], tol
   elseif nM == 0
      νl[1] = mM
      μl[1] = 0
      return ν[1:0], μ[1:0], nf, νl[1:1], μl[1:1], tol
   end

   mrinf = 0
   nrinf = 0
   i = 0
   tol = max(atol, rtol*opnorm(M,1))

   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      τ, ρ = _preduce1!(n, m, p, M, N, Q, Z, tol;
                        fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                        withQ = withQ, withZ = withZ)
      i += 1
      ν[i] = ρ+τ
      μ[i] = m
      roff += τ+ρ
      coff += m
      n -= ρ
      m = ρ
      p -= τ
   end
   j = 0
   while p > 0
      # Step 3: Particular form of the dual algorithm PREDUCE
      ρ = _preduce4!(n, 0, p, M, N, Q, Z, tol,
                     fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                     withQ = withQ, withZ = withZ)
      j += 1
      νl[j] = p
      μl[j] = ρ
      rtrail += p
      ctrail += ρ
      n -= ρ
      p = ρ
   end

   return ν[1:i], μ[1:i], n, reverse(νl[1:j]), reverse(μl[1:j]), tol
end
"""
    klf_left_refine!(ν, μ, M, N, tol; fast = true, ut = false, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true) -> (νi, νl, μl)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [   *         *        *  ] roff
    M - λN = [   0      Mil-λNil    *  ] mli
             [   0         0        *  ] rtrail
               coff      nli     ctrail

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation matrices `Q1` and `Z1`
such that the full column rank subpencil `Mil-λNil` is transformed into the
following Kronecker-like form  exhibiting its infinite and left Kronecker structures:

                  |  Mi-λNi    |     *      |
     Fil - λGil = |------------|------------|
                  |    O       |   Ml-λNl   |

The full column rank pencil `Mil-λNil` is in a staircase form and the `nb`-dimensional vectors `ν` and `μ`
contain the row and, respectively, column dimensions of the blocks of the staircase form  `Mil-λNil` such that
the `i`-th block has dimensions `ν[i] x μ[i]` and has full column rank. The matrix Mil has full column rank and
the leading `μ[1]+...+μ[nb-1]` columns of `Nil` form a full row rank submatrix.

The regular pencil `Mi-λNi` is in a staircase form, contains the infinite elementary divisors of `Mil-λNil`,
`Mi` is upper triangular if `ut = true` and nonsingular and `Ni` is upper triangular and nilpotent. The `nbi`-dimensional vector `νi` contains the dimensions of the
square diagonal blocks of the staircase form  `Mi-λNi` such that the `i`-th block has dimensions `νi[i] x νi[i]`.
The difference `νi[i]-νi[i-1] = ν[nbi-i+1]-μ[nbi-i+1]` for `i = 1, 2, ..., nbi` is the number of infinite elementary
divisors of degree `i` (with `νi[0] = 0` and `μ[nbi+1] = 0`).

The full column rank pencil `Mil-λNil` is in a staircase form, and contains the  left Kronecker indices
and infinite elementary divisors of the subpencil `M22 - λN22`.
The `nb`-dimensional vectors `ν` and `μ` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Mil-λNil` such that the `i`-th block has dimensions `ν[nb-i+1] x μ[nb-i+1]` and has full column rank.
The difference `ν[nb-i+1]-μ[nb-i+1]` for `i = 1, 2, ..., nb` is the number of elementary Kronecker blocks of size
`i x (i-1)`. The difference `μ[nb-i+1]-ν[nb-i]` for `i = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `i` (with `ν[0] = 0`).

The full column rank pencil `Ml-λNl`, in a staircase form, contains the left Kronecker indices of `M-λN` and has the form

               | Al-λEl |
     Ml-λNl  = |--------|,
               |   Cl   |

where `El` is upper triangular and nonsingular.
The `nl`-dimensional vectors `νl` and `μl` contain the row and, respectively, column dimensions of the blocks
of the staircase form  `Ml-λNl` such that the `j`-th block has dimensions `νl[j] x μl[j]` and has full column rank.
The difference `νl[nl-j+1]-μl[nl-j+1] = ν[nl-j+1]-μ[nl-j+1]` for `j = 1, 2, ..., nl` is the number of elementary
Kronecker blocks of size `j x (j-1)`.
If `ut = true`, the full column rank diagonal blocks of `Ml` are reduced to the form `[Y; 0]`
with `Y` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nl` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.

`F` and `G` are returned in `M` and `N`, respectively.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` (i.e., `Q <- Q*Q1`)
if `withQ = true`. Otherwise, `Q` is not modified.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` (i.e., `Z <- Z*Z1`)
if `withZ = true`. Otherwise, `Z` is not modified.
"""
function klf_left_refine!(ν::Vector{Int}, μ::Vector{Int}, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, Q::Union{AbstractMatrix{T1},Nothing},
   Z::Union{AbstractMatrix{T1},Nothing}, tol::Real; fast::Bool = true, ut::Bool = false, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
   withQ::Bool = true, withZ::Bool = true) where T1 <: BlasFloat
   nb = length(ν)
   nb == length(μ) || throw(DimensionMismatch("ν and μ must have the same lengths"))

   nb == 0 && (return ν[1:0], ν[1:0], ν[1:0])

   mli = sum(ν)
   nli = sum(μ)
   mM = roff + mli + rtrail
   nM = coff + nli + ctrail
   p = ν[nb]
   n = mli-p
   m = nli-n

   if n > 0
      # Step 0: Reduce N = [ 0 E11] to standard form, where E11 is full row rank
      #                    [ 0 0  ]
      it = roff+1:roff+n
      jt = coff+1:coff+nli
      tau = similar(N,n)
      E11 = view(N,it,jt)
      LinearAlgebra.LAPACK.gerqf!(E11,tau)
      eltype(M) <: Complex ? tran = 'C' : tran = 'T'
      LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(M,1:mM,jt))
      withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(Z,:,jt))
      LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(N,1:roff,jt))
      triu!(E11,m)
   end

   μl = Vector{Int}(undef,nb)
   νl = Vector{Int}(undef,nb)
   νi = Vector{Int}(undef,nb)

   roff0 = roff
   coff0 = coff

   i = 0
   ni = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm
      τ, ρ = _preduce1!(n, m, p, M, N, Q, Z, tol;
                        fast = fast, roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                        withQ = withQ, withZ = withZ)
      ρ+τ == m || error("The reduced pencil must not have right structure: try to adjust the tolerances")
      i += 1
      νi[i] = m
      ni += m
      roff += m
      coff += m
      n -= ρ
      m = ρ
      p -= τ
   end
   # make Mi upper triangular
   ut && ni > 0 && klf_left_refineinf!(view(νi,1:i), M, N, Q, missing; roff = roff0, coff = coff0, withQ = withQ)

   j = 0
   while p > 0
      # Step 3: Dual algorithm
      ρ = _preduce4!(n, 0, p, M, N, Q, Z, tol, fast = fast,
                     roff = roff, coff = coff, rtrail = rtrail, ctrail = ctrail,
                     withQ = withQ, withZ = withZ)
      j += 1
      νl[j] = p
      μl[j] = ρ
      rtrail += p
      ctrail += ρ
      n -= ρ
      p = ρ
   end
   reverse!(view(νl,1:j))
   reverse!(view(μl,1:j))
   # make diagonal blocks of Ml and supradiagonal blocks of Nl upper triangular
   ut && j > 0 && klf_left_refineut!(view(νl,1:j), view(μl,1:j), M, N, Q, Z; roff = roff0+ni, coff = coff0+ni, withQ = withQ, withZ = withZ)

   return νi[1:i], νl[1:j], μl[1:j]
end
"""
    klf_right_refineut!(ν, μ, M, N, Q, Z; ctrail = 0, withQ = true, withZ = true)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [ Mri-λNri  *   ]
    M - λN = [    0      *   ]
                       ctrail

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation
matrices  `Q1` and `Z1`, such that the full row rank subpencil `Mri-λNri`, in staircase form ,
is transformed as follows: the full row rank diagonal blocks of `Mri` are reduced to the form `[0 X]`
with `X` upper triangular and nonsingular and the full column rank supradiagonal blocks of `Nri` are
reduced to the form `[Y; 0]` with `Y` upper triangular and nonsingular.
It is assumed that `Mri-λNri` has `nb` diagonal blocks and the dimensions of the diagonal blocks
are specified by the `nb`-dimensional vectors `ν` and `μ` such that the `i`-th block has dimensions `ν[i] x μ[i]`.

The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`),
if `withQ = true` and `Z` (i.e., `Z <- Z*Z1`), if `withZ = true`.
"""
function klf_right_refineut!(ν::AbstractVector{Int}, μ::AbstractVector{Int}, M::AbstractMatrix{T}, N::AbstractMatrix{T},
                             Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing};
                             rtrail::Int = 0, ctrail::Int = 0, withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat

   nb = length(ν)
   nb == 0 && return
   nb == length(μ) || throw(DimensionMismatch("ν and μ must have the same lengths"))
   mri = sum(ν)
   nri = sum(μ)
   nM = nri + ctrail
   # mM, nM == size(M) || throw(DimensionMismatch("Incompatible ν, μ, rtrail and ctrail with the dimensions of M"))
   # mM, nM == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))
   ki2 = mri
   kj2 = nri
   T <: Complex ? tran = 'C' : tran = 'T'
   for i = nb:-1:1
      #i == 1 && return
      ni1 = ν[i]
      ki1 = ki2-ni1+1
      kki = ki1:ki2
      nj1 = μ[i]
      kj1 = kj2-nj1+1
      kkj = kj1:kj2
      if nj1 > 1
         Mij = view(M,kki,kkj)
         tau = similar(M,ni1)
         LinearAlgebra.LAPACK.gerqf!(Mij,tau)
         LinearAlgebra.LAPACK.ormrq!('R',tran,Mij,tau,view(M,1:ki1-1,kkj))
         LinearAlgebra.LAPACK.ormrq!('R',tran,Mij,tau,view(N,1:ki1-1,kkj))
         withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Mij,tau,view(Z,:,kkj))
         M[kki,kkj] = [ zeros(T,ni1,nj1-ni1) triu(view(Mij,:,nj1-ni1+1:nj1))]
      end
      if i > 1
         nim1 = ν[i-1]
         kim1 = ki1-nim1
         kim2 = ki1-1
         kkim = kim1:kim2
         if nim1 > 1
            Nimj = view(N,kkim,kkj)
            tau = similar(N,nj1)
            LinearAlgebra.LAPACK.geqrf!(Nimj,tau)
            njm1 = μ[i-1]
            kjm1 = kj1-njm1
            kjm2 = kj1-1
            LinearAlgebra.LAPACK.ormqr!('L',tran,Nimj,tau,view(M,kkim,kjm1:nM))
            LinearAlgebra.LAPACK.ormqr!('L',tran,Nimj,tau,view(N,kkim,kj2+1:nM))
            #withQ && LinearAlgebra.LAPACK.ormqr!('R',tran,Nimj,tau,view(Q,:,kkim))
            withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Nimj,tau,view(Q,:,kkim))
            triu!(Nimj)
         end
      end
      ki2 = ki1-1
      kj2 = kj1-1
   end
   return
end
"""
    klf_left_refineut!(ν, μ, M, N, Q, Z; roff = 0, coff = 0, withQ = true, withZ = true)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [  *     *     ] roff
    M - λN = [  0  Mil-λNil ]
              coff

to an equivalent form `F - λG = Q1'*(M - λN)*Z1` using orthogonal or unitary transformation
matrices  `Q1` and `Z1`, such that the full column rank subpencil `Mil-λNil`, in staircase form ,
is transformed as follows: the full column rank diagonal blocks of `Mil` are reduced to the form `[Y ; 0]`
with `Y` upper triangular and nonsingular and the full row rank supradiagonal blocks of `Nil` are
reduced to the form `[0 Y]` with `Y` upper triangular and nonsingular.
It is assumed that `Mil-λNil` has `nb` diagonal blocks and the dimensions of the diagonal blocks
are specified by the `nb`-dimensional vectors `ν` and `μ` such that the `i`-th block has dimensions `ν[i] x μ[i]`.

The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`),
if `withQ = true` and `Z` (i.e., `Z <- Z*Z1`), if `withZ = true`.
"""
function klf_left_refineut!(ν::AbstractVector{Int}, μ::AbstractVector{Int}, M::AbstractMatrix{T}, N::AbstractMatrix{T},
                             Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing};
                             roff::Int = 0, coff::Int = 0, withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat

   nb = length(ν)
   nb == 0 && return
   nb == length(μ) || throw(DimensionMismatch("ν and μ must have the same lengths"))
   mli = sum(ν)
   nli = sum(μ)
   mM = mli + roff
   nM = nli + coff
   # mM, nM == size(M) || throw(DimensionMismatch("Incompatible ν, μ, rtrail and ctrail with the dimensions of M"))
   # mM, nM == size(N) || throw(DimensionMismatch("M and N must have the same dimensions"))
   ki1 = roff+1
   kj1 = coff+1
   T <: Complex ? tran = 'C' : tran = 'T'
   for i = 1:nb
      ni1 = ν[i]
      ki2 = ki1+ni1-1
      kki = ki1:ki2
      nj1 = μ[i]
      kj2 = kj1+nj1-1
      kkj = kj1:kj2
      if ni1 > 1
         Mij = view(M,kki,kkj)
         tau = similar(M,nj1)
         LinearAlgebra.LAPACK.geqrf!(Mij,tau)
         LinearAlgebra.LAPACK.ormqr!('L',tran,Mij,tau,view(M,kki,kj2+1:nM))
         LinearAlgebra.LAPACK.ormqr!('L',tran,Mij,tau,view(N,kki,kj1+1:nM))
         withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Mij,tau,view(Q,:,kki))
         triu!(Mij)
      end
      if i < nb
         nip1 = ν[i+1]
         njp1 = μ[i+1]
         kjp1 = kj2+1
         kjp2 = kj2+njp1
         kkjp = kjp1:kjp2
         if nip1 > 1
            Nijp = view(N,kki,kkjp)
            tau = similar(N,ni1)
            LinearAlgebra.LAPACK.gerqf!(Nijp,tau)
            nip1 = ν[i+1]
            kip1 = ki2+1
            kip2 = ki2+nip1
            LinearAlgebra.LAPACK.ormrq!('R',tran,Nijp,tau,view(M,1:kip2,kkjp))
            LinearAlgebra.LAPACK.ormrq!('R',tran,Nijp,tau,view(N,kip1:kip2,kkjp))
            LinearAlgebra.LAPACK.ormrq!('R',tran,Nijp,tau,view(N,1:ki1-1,kkjp))
            withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Nijp,tau,view(Z,:,kkjp))
            N[kki,kkjp] = [ zeros(T,ni1,njp1-ni1) triu(view(Nijp,:,njp1-ni1+1:njp1))]
         end
      end
      ki1 = ki2+1
      kj1 = kj2+1
   end
   return
end
"""
    klf_right_refineinf!(νi, M, N, Z, R; roff = 0, coff = 0, withZ = true)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [   *         *        *  ] roff
    M - λN = [   0       Mi-λNi     *  ] ni
             [   0         0        *  ] rtrail
               coff       ni     ctrail

to an equivalent form `F - λG = (M - λN)*Z1` using an orthogonal or unitary transformation
matrix `Z1`, such that the regular subpencil `Mi-λNi`, in staircase form ,
with `Mi` nonsingular and `Ni` nillpotent and upper triangular, is transformed to obtain
`Mi` is upper-triangular and preserve `Ni` upper triangular.
It is assumed that `Mi-λNi` has `nb` diagonal blocks and the `i`-th diagonal block has dimensions
`νi[i] x νi[i]`.

The performed orthogonal or unitary transformations are accumulated in `Z` (i.e., `Z <- Z*Z1`),
if `withZ = true`.

The  matrix `R` is overwritten by `R*Z1` unless `R = missing`.
"""
function klf_right_refineinf!(νi::AbstractVector{Int}, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, Z::Union{AbstractMatrix{T1},Nothing},
   R::Union{AbstractMatrix{T1},Missing}; roff::Int = 0, coff::Int = 0, withZ::Bool = true) where T1 <: BlasFloat
   nb = length(νi)
   nb == 0 && return
   T1 <: Complex ? tran = 'C' : tran = 'T'
   ni = sum(νi)

   ki2 = roff+ni
   kj2 = coff+ni
   for k = nb:-1:1
       nk = νi[k]
       ki1 = ki2-nk+1
       kj1 = kj2-nk+1
       kki = ki1:ki2
       kkj = kj1:kj2
       if nk > 1
          Mk = view(M,kki,kkj)
          tau = similar(M,nk)
          LinearAlgebra.LAPACK.gerqf!(Mk,tau)
          LinearAlgebra.LAPACK.ormrq!('R',tran,Mk,tau,view(M,1:ki1-1,kkj))
          withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Mk,tau,view(Z,:,kkj))
          LinearAlgebra.LAPACK.ormrq!('R',tran,Mk,tau,view(N,1:ki1-1,kkj))
          ismissing(R) || LinearAlgebra.LAPACK.ormrq!('R',tran,Mk,tau,view(R,:,kkj))
          triu!(Mk)
       end
       ki2 = ki1-1
       kj2 = kj1-1
   end
   return
end
"""
    klf_left_refineinf!(νi, M, N, Q, L; roff = 0, coff = 0, withQ = true)

Reduce the partitioned matrix pencil `M - λN` (`*` stands for a not relevant subpencil)

             [   *         *        *  ] roff
    M - λN = [   0       Mi-λNi     *  ] ni
             [   0         0        *  ] rtrail
               coff       ni     ctrail

to an equivalent form `F - λG = Q1'*(M - λN)` using an orthogonal or unitary transformation
matrix `Q1`, such that the regular subpencil `Mi-λNi`, in staircase form,
with `Mi` nonsingular and `Ni` nillpotent and upper triangular, is transformed to obtain
`Mi` upper-triangular and preserve `Ni` upper triangular.
It is assumed that `Mi-λNi` has `nb` diagonal blocks and the `i`-th diagonal block has dimensions
`νi[i] x νi[i]`.

The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`),
if `withQ = true`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing`.
"""
function klf_left_refineinf!(νi::AbstractVector{Int}, M::AbstractMatrix{T1}, N::AbstractMatrix{T1}, Q::Union{AbstractMatrix{T1},Nothing},
   L::Union{AbstractMatrix{T1},Missing}; roff::Int = 0, coff::Int = 0, withQ::Bool = true) where T1 <: BlasFloat
   nb = length(νi)
   nb == 0 && return
   T1 <: Complex ? tran = 'C' : tran = 'T'
   n = size(M,2)

   ki1 = roff+1
   kj1 = coff+1
   for k = 1:nb
       nk = νi[k]
       ki2 = ki1+nk-1
       kj2 = kj1+nk-1
       kki = ki1:ki2
       kkj = kj1:kj2
       if nk > 1
          Mk = view(M,kki,kkj)
          tau = similar(M,nk)
          LinearAlgebra.LAPACK.geqrf!(Mk,tau)
          LinearAlgebra.LAPACK.ormqr!('L',tran,Mk,tau,view(M,kki,kj2+1:n))
          withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Mk,tau,view(Q,:,kki))
          ismissing(L) || LinearAlgebra.LAPACK.ormqr!('L',tran,Mk,tau,view(L,kki,:))
          LinearAlgebra.LAPACK.ormqr!('L',tran,Mk,tau,view(N,kki,kj2+1:n))
          triu!(Mk)
       end
       ki1 = ki2+1
       kj1 = kj2+1
   end
   return
end


text:
using Plots
using DifferentialEquations

@userplot threebodyplot
@recipe function f(tbp::threebodyplot)
    x, y, l, xlim, ylim = tbp.args
    xaxis --> ("x", xlim)
    yaxis --> ("y", ylim)
    seriestype --> :scatter
    markertype --> :o
    markercolor --> "black"
    label --> l
    dpi --> 300
    [x], [y]
end

function anim()

    function inertial(ddu, du, u, p, t)

        μ1 = p[1,1]
        μ2 = p[2,1]
        x1 = p[:,2]
        x2 = p[:,3]

        r1 = sqrt((x1[1] - u[1])^2 + (x1[2] - u[2])^2 )
        r2 = sqrt((x2[1] - u[1])^2 + (x2[2] - u[2])^2 )

        ddu[1] = μ1*(x1[1] - u[1])/r1^3 + μ2*(x2[1] - u[1])/r2^3
        ddu[2] = μ1*(x1[2] - u[2])/r1^3 + μ2*(x2[2] - u[2])/r2^3

    end

    function rotating(ddu, du, u, p, t)

        μ1 = p[1]
        μ2 = p[2]
        n = p[3]

        r1 = sqrt((u[1] + μ2)^2 + u[2]^2)
        r2 = sqrt((u[1] - μ1)^2 + u[2]^2)

        ddu[1] = -(μ1*(u[1] + μ2)/r1^3 + μ2*(u[1] - μ1)/r2^3) + 2n*du[2] + n^2*u[1]
        ddu[2] = -u[2]*(μ1/r1^3 + μ2/r2^3) - 2n*du[1] + n^2*u[2]

    end

    function hill(ddu, du, u, p, t)

        f = -3/(u[1]^2 + u[2]^2)^1.5

        ddu[1] = u[1]*(f + 3) + 2du[2]
        ddu[2] = u[2]*f - 2du[1]

    end

    t0 = -100.
    t1 = 150.
    tspan = (t0, t1)

    mu1 = .2
    mu2 = .8
    n = 1
    x1 = [1, 1]
    x2 = [-1, -1]

    p = zeros(2,3)
    p[1,1] = mu1
    p[2,1] = mu2
    p[:,2] = x1
    p[:,3] = x2

    du0 = zeros(2)
    u0 = zeros(2)
    u0[1] = 0.5

    probi = SecondOrderODEProblem(inertial, du0, u0, tspan, p, abstol=1e-8, reltol=1e-8)
    soli = solve(probi, Tsit5())

    p = zeros(3)
    p[1] = mu1
    p[2] = mu2
    p[3] = n

    dv0 = zeros(2)
    v0 = zeros(2)
    v0[1] = u0[1]*cos(n*t0) + u0[2]*sin(n*t0)
    v0[2] = -u0[1]*sin(n*t0) + u0[2]*cos(n*t0)

    probr = SecondOrderODEProblem(rotating, dv0, v0, tspan, p, abstol=1e-8, reltol=1e-8)
    solr = solve(probr, Tsit5())

    dw0 = zeros(2)
    w0 = zeros(2)
    w0[1] = v0[1] + 1
    w0[2] = v0[2]

    probh = SecondOrderODEProblem(hill, dw0, w0, tspan, abstol=1e-8, reltol=1e-8)
    solh = solve(probr, Tsit5())

    n = min(length(solh.t), length(solr.t), length(soli.t))

    anim = @animate for i in 1:10:n

        p1 = plot((x1[1], x1[2]), markershape=:o, markercolor="blue", label="planet 1",
                  legend=:bottomright, framestyle=:zerolines)
        p1 = plot!((x2[1], x2[2]), markershape=:o, markercolor="green", label="planet 2")
        p1 = threebodyplot!(soli[3,i], soli[4,i], "intertial", (-1.5,1.5), (-1.5,1.5))

        p2 = plot((-mu2, 0), markershape=:o, markercolor="blue",  label="planet 1",
                  legend=:bottomright, framestyle=:zerolines)
        p2 = plot!((mu1, 0), markershape=:o, markercolor="green", label="planet 2")
        p2 = threebodyplot!(solr[3,i], solr[4,i], "rotating", (-1, 1), (-1, 1))

        p3 = plot((0, 0), markershape=:o, markercolor="green", label="planet 2",
                  legend=:bottomright, framestyle=:zerolines)
        p3 = threebodyplot!(solh[3,i], solh[4,i], "hill", (-1.5,1.5), (-1.5,1.5))

        plot(p1, p2, p3, layout=(3,1), size=(800,1600))
    end

    gif(anim, "de_anim_15.gif", fps=15)

end

anim()

ext: jl

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
function get_eccentricities(obj::FrequencyPattern)
    return jcall(obj, "getEccentricities", Vector3D, ())
end

function get_phase_center_variation(obj::FrequencyPattern, arg0::Vector3D)
    return jcall(obj, "getPhaseCenterVariation", jdouble, (Vector3D,), arg0)
end


text:
# This file roughly corresponds to functions documented in the
# Files API: https://canvas.instructure.com/doc/api/files

"""
    Canvas.quota(co::Union{Course,Group,User}; kwargs...) -> Dict

Return the total and used storage quota for the course/group/user.
Return a dictionary.

**Canvas API documentation**:
[Get quota information (`GET /api/v1/courses/:course_id/files/quota`,
`GET /api/v1/groups/:group_id/files/quota`, `GET /api/v1/users/:user_id/files/quota`)]
(https://canvas.instructure.com/doc/api/files#method.files.api_quota)
"""
function quota(co::Union{Course,Group,User}; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(co))/files/quota"; kwargs...)
    return json
end

"""
    Canvas.files(co::Union{Course,Group,User,Folder}; kwargs...) -> Vector{File}, page_data

Return a paginated  list of files for the course/groupt/user/folder.
Return a vector of [`File`](@ref)s and a dictionary of page data.

**Canvas API documentation**:
[List files (`GET /api/v1/courses/:course_id/files`, `GET /api/v1/users/:user_id/files`,
`GET /api/v1/groups/:group_id/files`, `GET /api/v1/folders/:id/files`)]
(https://canvas.instructure.com/doc/api/files#method.files.api_index)
"""
function files(co::Union{Course,Group,User,Folder}; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(co))/files"; kwargs...)
    return File.(json), page_data
end

"""
    Canvas.preview_url(f::File; kwargs...) -> Dict

Determine the URL that should be used for inline preview of the file.
Return a dictionary.

**Canvas API documentation**:
[Get public inline preview url (`GET /api/v1/files/:id/public_url`)]
(https://canvas.instructure.com/doc/api/files#method.files.public_url)
"""
function preview_url(f::File; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(f))/public_url"; kwargs...)
    return json
end

"""
    Canvas.file(f::File; kwargs...) -> File

Get the details of a single file. Return a [`File`](@ref).

**Canvas API documentation**:
[Get file (`GET /api/v1/files/:id`)]
(https://canvas.instructure.com/doc/api/files#method.files.api_show)
"""
function file(f::File; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(f))"; kwargs...)
    return File(json)
end

"""
    Canvas.file(co::Union{Course,Group,User}, f::File; kwargs...) -> File

Get the details of a single file. Return a [`File`](@ref).

**Canvas API documentation**:
[Get file (`GET /api/v1/courses/:course_id/files/:id`,
`GET /api/v1/groups/:group_id/files/:id`, `GET /api/v1/users/:user_id/files/:id`)]
(https://canvas.instructure.com/doc/api/files#method.files.api_show)
"""
function file(co::Union{Course,Group,User}, f::File; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(co))$(Internals.pid(f))"; kwargs...)
    return File(json)
end

"""
    Canvas.update_file(f::File; kwargs...) -> File

Update settings for the specified file. Return a [`File`](@ref).

**Canvas API documentation**:
[Update file (`PUT /api/v1/files/:id`)](https://canvas.instructure.com/doc/api/files#method.files.api_update)
"""
function update_file(f::File; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(f))"; kwargs...)
    return File(json)
end

"""
    Canvas.delete_file(f::File; kwargs...) -> File

Delete the specified file and return the details. Return a [`File`](@ref).

**Canvas API documentation**:
[Delete file (`DELETE /api/v1/files/:id`)](https://canvas.instructure.com/doc/api/files#method.files.destroy)
"""
function delete_file(f::File; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(f))"; kwargs...)
    return File(json)
end

"""
    Canvas.folders(f::Folder; kwargs...) -> Vector{Folder}, page_data

Return a paginated list of folders in the specified folder.
Return a vector of [`Folder`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List folders (`GET /api/v1/folders/:id/folders`)]
(https://canvas.instructure.com/doc/api/files#method.folders.api_index)
"""
function folders(f::Folder; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(f))/folders"; kwargs...)
    return Folder.(json), page_data
end

"""
    Canvas.folders(co::Union{Course,Group,User}; kwargs...) -> Vector{Folder}, page_data

Return a paginated, flattened, list of all folders in the course/group/user.
Return a vector of [`Folder`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List all folders (`GET /api/v1/courses/:course_id/folders`,
`GET /api/v1/groups/:group_id/folders`, `GET /api/v1/users/:user_id/folders`)]
(https://canvas.instructure.com/doc/api/files#method.folders.list_all_folders)
"""
function folders(co::Union{Course,Group,User}; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(co))/folders"; kwargs...)
    return Folder.(json), page_data
end

"""
    Canvas.resolve_path(co::Union{Course,Group,User}, path=""; kwargs...) -> Vector{Folder}

Return the list of folders from the root of the course/group/user to the path `path`.
Return a vector of [`Folder`](@ref)s.

**Canvas API documentation**:
[Resolve path (`GET /api/v1/courses/:course_id/folders/by_path`,
`GET /api/v1/groups/:group_id/folders/by_path`, `GET /api/v1/users/:user_id/folders/by_path`)]
(https://canvas.instructure.com/doc/api/files#method.folders.resolve_path)
"""
function resolve_path(co::Union{Course,Group,User}, path=""; kwargs...)
    path = isempty(path) ? "" : "/$(lstrip(path, '/'))"
    json = request("GET", "/api/v1$(Internals.pid(co))/folders/by_path$(path)"; kwargs...)
    return Folder.(json)
end

"""
    Canvas.folder(f::Folder; kwargs...) -> Folder

Return the details of a folder specified by id. Return a [`Folder`](@ref).

**Canvas API documentation**:
[Get folder (`GET /api/v1/folders/:id`)]
(https://canvas.instructure.com/doc/api/files#method.folders.show)
"""
function folder(f::Folder; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(f))"; kwargs...)
    return Folder(json)
end

# TODO: allow "root" as folder id to get the root folder
"""
    Canvas.folder(co::Union{Course,Group,User}, f::Folder; kwargs...) -> Folder

Return the details of a folder specified by id. Return a [`Folder`](@ref).

**Canvas API documentation**:
[Get folder (`GET /api/v1/courses/:course_id/folders/:id`,
`GET /api/v1/groups/:group_id/folders/:id`, `GET /api/v1/users/:user_id/folders/:id`)]
(https://canvas.instructure.com/doc/api/files#method.folders.show)
"""
function folder(co::Union{Course,Group,User}, f::Folder; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(co))$(Internals.pid(f))"; kwargs...)
    return Folder(json)
end

"""
    Canvas.update_folder(f::Folder; kwargs...) -> Folder

Update settings for the folder. Return a [`Folder`](@ref).

**Canvas API documentation**:
[Update folder (`PUT /api/v1/folders/:id`)]
(https://canvas.instructure.com/doc/api/files#method.folders.update)
"""
function update_folder(f::Folder; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(f))"; kwargs...)
    return Folder(json)
end

"""
    Canvas.create_folder(co::Union{Course,Group,User,Folder}; kwargs...) -> Folder

Create a new folder in the course/group/user/folder. Return a [`Folder`](@ref).

**Canvas API documentation**:
[Create folder (`POST /api/v1/courses/:course_id/folders`, `POST /api/v1/groups/:group_id/folders`,
`POST /api/v1/users/:user_id/folders`, `POST /api/v1/folders/:folder_id/folders`)]
(https://canvas.instructure.com/doc/api/files#method.folders.create)
"""
function create_folder(f::Folder; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(f))/folders"; kwargs...)
    return Folder(json)
end

"""
    Canvas.delete_folder(f::Folder; kwargs...) -> Dict

Delete the folder. Return a dictionary.

**Canvas API documentation**:
[Delete folder (`DELETE /api/v1/folders/:id`)]
(https://canvas.instructure.com/doc/api/files#method.folders.api_destroy)
"""
function delete_folder(f::Folder; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(f))"; kwargs...)
    return json # ??
end

"""
    Canvas.upload_file(f::Folder, file::String; kwargs...) -> File

Upload a file to a folder. Return a [`File`](@ref).

**Canvas API documentation**:
[Upload a file (`POST /api/v1/folders/:folder_id/files`)]
(https://canvas.instructure.com/doc/api/files#method.folders.create_file)
"""
function upload_file(f::Folder, file::String; kwargs...)
    return Internals.upload_file("/api/v1$(Internals.pid(f))/files", file; kwargs...)
end

"""
    Canvas.copy_file(f::Folder; kwargs...) -> File

Copy a file into the specified folder. Return a [`File`](@ref).

**Canvas API documentation**:
[Copy a file (`POST /api/v1/folders/:dest_folder_id/copy_file`)]
(https://canvas.instructure.com/doc/api/files#method.folders.copy_file)
"""
function copy_file(f::Folder; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(f))/copy_file"; kwargs...)
    return File(json)
end

"""
    Canvas.copy_folder(f::Folder; kwargs...) -> Folder

Copy a folder and its contents into the specified folder. Return a [`Folder`](@ref).

**Canvas API documentation**:
[Copy a folder (`POST /api/v1/folders/:dest_folder_id/copy_folder`)]
(https://canvas.instructure.com/doc/api/files#method.folders.copy_folder)
"""
function copy_folder(f::Folder; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(f))/copy_folder"; kwargs...)
    return Folder(json)
end

"""
    Canvas.media_folder(co::Union{Course,Group}; kwargs...) -> Folder

Return the details for a designated upload folder that the user has rights to upload to.
Return a [`Folder`](@ref).

**Canvas API documentation**:
[Get uploaded media folder for user (`GET /api/v1/courses/:course_id/folders/media`, `GET /api/v1/groups/:group_id/folders/media`)]
(https://canvas.instructure.com/doc/api/files#method.folders.media_folder)
"""
function media_folder(co::Union{Course,Group}; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(co))/folders/media"; kwargs...)
    return Folder(json)
end

"""
    Canvas.update_usage_rights(co::Union{Course,Group,User}; kwargs...) -> UsageRights

Set copyright and license information for one or more files. Return a [`UsageRights`](@ref).

**Canvas API documentation**:
[Set usage rights (`PUT /api/v1/courses/:course_id/usage_rights`,
`PUT /api/v1/groups/:group_id/usage_rights`, `PUT /api/v1/users/:user_id/usage_rights`)]
(https://canvas.instructure.com/doc/api/files#method.usage_rights.set_usage_rights)
"""
function update_usage_rights(co::Union{Course,Group,User}; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(co))/usage_rights"; kwargs...)
    return UsageRights(json)
end

"""
    Canvas.remove_usage_rights(co::Union{Course,Group,User}; kwargs...) -> Dict

Remove copyright and license information associated with one or more files.

**Canvas API documentation**:
[Remove usage rights (`DELETE /api/v1/courses/:course_id/usage_rights`,
`DELETE /api/v1/groups/:group_id/usage_rights`, `DELETE /api/v1/users/:user_id/usage_rights`)]
(https://canvas.instructure.com/doc/api/files#method.usage_rights.remove_usage_rights)
"""
function remove_usage_rights(co::Union{Course,Group,User}; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(co))/usage_rights"; kwargs...)
    return json # ??
end

"""
    Canvas.licenses(co::Union{Course,Group,User}; kwargs...) -> Vector{License}, page_data

Return a paginated list of licenses that can be applied to files.
Return a vector of [`License`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List licenses (`GET /api/v1/courses/:course_id/content_licenses`,
`GET /api/v1/groups/:group_id/content_licenses`, `GET /api/v1/users/:user_id/content_licenses`)]
(https://canvas.instructure.com/doc/api/files#method.usage_rights.licenses)
"""
function licenses(co::Union{Course,Group,User}; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(co))/content_licenses"; kwargs...)
    return License.(json), page_data
end

# Helper function which is not really part of the upstream API
"""
    Canvas.download_file(f::File, local_path=mktempdir(); update_period=Inf, kwargs...)

Download the [`File`](@ref). If `local_path` is a directory the file will be
placed in this directory with the upstream filename. If `local_path` is a
file path the file will be downloaded to that path directly, ignoring the
upstream filename.
"""
function download_file(f::File, local_path=mktempdir(); update_period=Inf, kwargs...)
    if isdir(local_path) # existing directory, append filename
        local_path = joinpath(local_path, f.filename)
    elseif !isdir(dirname(local_path)) # a filepath, make sure directory exist
        error("destination directory does not exist; $(dirname(local_path))")
    end
    headers = Internals.canvas_headers()
    r = HTTP.download(f.url, local_path, headers; update_period=update_period, kwargs...)
    return r
end


text:
export solve_vf_all!

function solve_vf_all!(evs::dcdp_Emax, t::dcdp_tmpvars, p::dcdp_primitives, θt::AbstractVector, σ::Real, itype::Tuple, dograd::Bool; kwargs...)

    solve_vf_terminal!(evs,    p)
    solve_vf_infill!(  evs, t, p, θt, σ, dograd, itype; kwargs...)
    learningUpdate!(   evs, t, p,     σ, dograd)
    solve_vf_explore!( evs, t, p, θt, σ, dograd, itype; kwargs...)
end

function solve_vf_all!(evs, t, p, θ, itype, dograd; kwargs...)
    solve_vf_all!(evs, t, p, _θt(θ, p), _σv(θ), itype, dograd; kwargs...)
end


text:
# toggle output relay
# KO
export psp_togglerelay

const TOGGLE_RELAY = [0x4b, 0x4F, 0x0d]

"""
    psp_togglerelay(io)

Toggle output relay.
"""
function psp_togglerelay(io_psp::IO)
  write(io_psp, TOGGLE_RELAY)
  return nothing
end


text:
function fill_cpmodel!(optimizer::Optimizer)
    # Adding variables
    bridge_variables!(optimizer)

    # Adding affine functions
    bridge_affines!(optimizer)

    # Adding constraints
    bridge_constraints!(optimizer)

    # Adding objective
    bridge_objective!(optimizer)

    optimizer
end

function bridge_objective!(optimizer::Optimizer)
    if !isnothing(optimizer.moimodel.objective_identifier)
        SeaPearl.addObjective!(optimizer.cpmodel, get_cp_variable(optimizer, optimizer.moimodel.objective_identifier))
    end
end

function get_cp_variable(optimizer::Optimizer, index::MOI.VariableIndex)
    cp_identifier = optimizer.moimodel.variables[index.value].cp_identifier
    return optimizer.cpmodel.variables[cp_identifier]
end

function get_cp_variable(optimizer::Optimizer, index::AffineIndex)
    cp_identifier = optimizer.moimodel.affines[index.value].cp_identifier
    return optimizer.cpmodel.variables[cp_identifier]
end

function bridge_constraints!(optimizer::Optimizer)
    for moiconstraint in optimizer.moimodel.constraints
        constraint = create_CPConstraint(moiconstraint, optimizer)
        if isnothing(constraint)
            continue
        end
        # add constraint to the model
        push!(optimizer.cpmodel.constraints, constraint)
    end
end

function bridge_variables!(optimizer::Optimizer)
    i = 1
    for x in optimizer.moimodel.variables
        @assert !isnothing(x.min) "Every variable must have a lower bound"
        @assert !isnothing(x.max) "Every variable must have an upper bound"
        x.cp_identifier = string(length(keys(optimizer.cpmodel.variables)) + 1)
        newvariable = SeaPearl.IntVar(x.min, x.max, x.cp_identifier, optimizer.cpmodel.trailer)
        SeaPearl.addVariable!(optimizer.cpmodel, newvariable)
        i += 1
    end
end

function build_affine_term!(optimizer::Optimizer, vi::MOI.VariableIndex, coeff::Float64)
    @assert isinteger(coeff) "You can't give a float coefficient."

    intCoeff = convert(Int, coeff)

    @assert intCoeff != 0 "Coefficient cannot be null."

    if intCoeff == 1
        return get_cp_variable(optimizer, vi)
    end

    if intCoeff < 0
        x = build_affine_term!(optimizer, vi, -coeff)
        new_id = string(length(keys(optimizer.cpmodel.variables)) + 1)
        new_var = IntVarViewOpposite(x, new_id)
        addVariable!(optimizer.cpmodel, new_var)
        return new_var
    end

    x = build_affine_term!(optimizer, vi, 1.)
    new_id = string(length(keys(optimizer.cpmodel.variables)) + 1)
    new_var = IntVarViewMul(x, intCoeff, new_id)
    addVariable!(optimizer.cpmodel, new_var)
    return new_var
end

function build_affine!(optimizer::Optimizer, aff_function::MOIAffineFunction)
    if !isnothing(aff_function.cp_identifier)
        return aff_function.cp_identifier
    end

    vars = Vector{AbstractIntVar}(undef, length(aff_function.content.terms) + 1)
    for i in 1:length(aff_function.content.terms)
        vi = aff_function.content.terms[i].variable_index
        coeff = aff_function.content.terms[i].coefficient
        x = build_affine_term!(optimizer, vi, coeff)
        vars[i] = x
    end

    # Dealing with the constant
    c = aff_function.content.constant
    c = convert(Int, c)
    if c != 0
        x = vars[length(aff_function.content.terms)]
        new_id = string(length(keys(optimizer.cpmodel.variables)) + 1)
        new_var = IntVarViewOffset(x, c, new_id)
        addVariable!(optimizer.cpmodel, new_var)
        vars[length(aff_function.content.terms)] = new_var
    end

    # Creating the variable that will be equal to the affine function
    minSum, maxSum = 0, 0
    for i in 1:length(aff_function.content.terms)
        minSum += minimum(vars[i].domain)
        maxSum += maximum(vars[i].domain)
    end
    new_id = string(length(keys(optimizer.cpmodel.variables)) + 1)
    sum = IntVar(minSum, maxSum, new_id, optimizer.cpmodel.trailer)
    addVariable!(optimizer.cpmodel, sum)

    # Creating the constraint that will make it equal
    new_id = string(length(keys(optimizer.cpmodel.variables)) + 1)
    lastForSumToZero = IntVarViewOpposite(sum, new_id)
    addVariable!(optimizer.cpmodel, lastForSumToZero)
    vars[end] = lastForSumToZero
    constraint = SumToZero(vars, optimizer.cpmodel.trailer)
    push!(optimizer.cpmodel.constraints, constraint)

    # Returning the CP identifier of the variable that is equal to the affine function (and storing it, a bit ugly)
    aff_function.cp_identifier = sum.id
    sum.id
end

function bridge_affines!(optimizer::Optimizer)
    for aff in optimizer.moimodel.affines
        build_affine!(optimizer, aff)
    end
end


text:
include("imports.jl")

##

@testset ExtendedTestSet "FastAI.jl" begin
    @testset ExtendedTestSet "encodings/" begin
        @testset ExtendedTestSet "projective.jl" begin
            include("encodings/projective.jl")
        end
        @testset ExtendedTestSet "imagepreprocessing.jl" begin
            include("encodings/imagepreprocessing.jl")
        end
        @testset ExtendedTestSet "keypointpreprocessing.jl" begin
            include("encodings/keypointpreprocessing.jl")
        end
    end

    @testset ExtendedTestSet "methods/" begin
        @testset ExtendedTestSet "imageclassification.jl" begin
            include("methods/imageclassification.jl")
        end
        @testset ExtendedTestSet "imagesegmentation.jl" begin
            include("methods/imagesegmentation.jl")
        end
        @testset ExtendedTestSet "singlekeypointregression.jl" begin
            include("methods/singlekeypointregression.jl")
        end
    end

    @testset ExtendedTestSet "datasets/" begin
        @testset ExtendedTestSet "transformations.jl" begin
            include("datasets/transformations.jl")
        end
        @testset ExtendedTestSet "containers.jl" begin
            include("datasets/containers.jl")
        end
    end

    @testset ExtendedTestSet "training/" begin
        @testset ExtendedTestSet "paramgroups.jl" begin
            include("training/paramgroups.jl")
        end
        @testset ExtendedTestSet "discriminativelrs.jl" begin
            include("training/discriminativelrs.jl")
        end
        @testset ExtendedTestSet "fitonecycle.jl" begin
            include("training/fitonecycle.jl")
        end
        @testset ExtendedTestSet "finetune.jl" begin
            include("training/finetune.jl")
        end
        @testset ExtendedTestSet "lrfind.jl" begin
            include("training/lrfind.jl")
        end
        # TODO: test learning rate finder
    end
end


text:
using CSVFiles
using DataFrames
using FileIO
using FLAC
using WAV

export normalize, normalize!, readwav, readflac, readwav_rfft, writewav, readfrd

function normalize!(a::SampleArray)
    a ./= maximum(abs.(a))
    a
end

normalize(a::SampleArray) = normalize!(similar(a))

function readwav(fname; format="double", names=nothing)
    y, Fs, nbits, opt = wavread(fname; format=format)
    SampleArray(y, Fs, names)
end

function readflac(fname; names=nothing)
    y, Fs = load(fname)
    SampleArray(y, Fs, names)
end

function readwav_rfft(fname; normalized=false, format="double", names=nothing)
    x = readwav(fname; format=format, names=names)
    if normalized
        normalize!(x)
    end
    x, rfft(x)
end

function readflac_rfft(fname; normalized=false, names=nothing)
    x = readflac(fname; names=names)
    if normalized
        normalize!(x)
    end
    x, rfft(x)
end

writewav(x::SampleArray, fname::String; kwargs...) = wavwrite(data(x), fname; Fs=rate(x), kwargs...)

function readfrd(fname, rate; name=nothing, magf=db2amp, phasef=deg2rad)
    df = DataFrame(load(File(format"TSV", fname); header_exists=false, colnames=[:freq, :mag, :phase], spacedelim=true))
    # df = DataFrames.readtable(fname, header=false, colnames=[:freq, :mag, :phase])
    sort!(df, [:freq])
    df = df[df.freq .<= toHz(rate)/2, :]
    SpectrumArray(Complex.(MagPhase.(magf.(df.mag), phasef.(df.phase))), rate, df.freq, isnothing(name) ? nothing : [name])
    # SpectrumArray(MagPhase.(magf.(df.mag), phasef.(df.phase)), rate, df.freq, isnothing(name) ? nothing : [name])
end


text:
using InvariantMeasures
using ValidatedNumerics

@testset "Dynamics" begin

D = Mod1Dynamic(x->2*x)

@test D.T(0.1) == 0.2

@test dfly(Lipschitz, L1, D) == (0.5, 0.0)
@test dfly(TotalVariation, L1, D) == (0.5, 0.0)

D = PwMap([ x-> x^2+0.25, x -> 4*x-2, x -> 4*x-3], [0, 0.5, 0.75, 1])

@test InvariantMeasures.DynamicDefinition.derivative(D, 0.1..0.1) ≈ 0.2
@test InvariantMeasures.DynamicDefinition.derivative(D, 0.2..0.3) ≈ 0.4..0.6
@test InvariantMeasures.DynamicDefinition.derivative(D, 0.4..0.6) ≈ 0.8..4
@test InvariantMeasures.DynamicDefinition.derivative(D, 0.7..0.8) ≈ 4

D = PwMap([ x-> 4*x, x -> 2*x-0.5, x -> 4*x-3], [0, 0.25, 0.75, 1])

@test InvariantMeasures.DynamicDefinition.derivative(D, 0.1..0.3) ≈ 2..4
@test InvariantMeasures.DynamicDefinition.derivative(D, 0.1..0.2) ≈ 4..4

D = PwMap([x->2*x, x->2-2*x], [@interval(0), @interval(0.5), @interval(1)])

@test InvariantMeasures.orientation(D, 1) == 1
@test InvariantMeasures.orientation(D, 2) == -1

D = mod1_dynamic(x -> 3.5x, (0,1))

@test D.endpoints ≈ [0, 2/7, 4/7, 6/7, 1]
@test D.Ts[2](0.5) == 0.75
@test D.y_endpoints == [0 1; 0 1; 0 1; 0 0.5]
@test D.increasing == [1,1,1,1]

D = mod1_dynamic(x -> 3.5x + 0.5)
@test D.endpoints ≈ [0, 1/7, 3/7, 5/7, 1]
@test D.Ts[end](1) == 1
@test D.y_endpoints == [0.5 1; 0 1; 0 1; 0 1]
@test D.increasing == [1,1,1,1]

D = mod1_dynamic(x -> -3.5x + 0.5)
@test D.endpoints ≈ [0, 1/7, 3/7, 5/7, 1]
@test D.Ts[end](5/7) == 1
@test D.y_endpoints == [0.5 0; 1 0; 1 0; 1 0]
@test D.increasing == [0,0,0,0]

D = Iterate(mod1_dynamic(x->2*x), 3)

using InvariantMeasures.DynamicDefinition: derivative, distorsion, nbranches

@test derivative(D, 0.1..0.2) == 8..8
@test distorsion(D, 0.1..0.2) == 0..0
@test map(k->preim(D, k, 0.5), 1:8) == [1//16, 3//16, 5//16, 7//16, 9//16, 11//16, 13//16, 15//16]

# Lanford map, so that we test also something that is not linear
f = x->2*x+0.5*x*(1-x)
D = Iterate(mod1_dynamic(f), 3)
# Interval Newton should converge in 4 steps here, so we set max_iter = 4
# to make sure that we haven't inadvertently implemented something that is not quadratically convergent
# (for instance, by using the wrong Jacobian)
its = map(k->preim(D, k, 0.5, 1e-15, max_iter = 4), 1:nbranches(D))
g(x) = f(x) - floor(f(x))
@test g.(g.(g.(its))) ≈ fill(0.5, 8)

@test endpoints(D) ≈ [0.0, 0.07389363935392047, 0.18200396341631753, 0.28117073603385473, 0.4384471871911697, 0.5287080193012084, 0.6633983486837269, 0.7902622123165944, 1.0]

@test [branch(D, k)(Interval(0.2)) for k in 1:nbranches(D)] ≈ [∅, ∅, g(g(g(Interval(0.2)))), ∅, ∅, ∅, ∅, ∅]

end


text:
module Cards

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export card, card_section, card_actions
export cardsection, cardactions

register_normal_element("q__card", context = @__MODULE__)
register_normal_element("q__card__section", context = @__MODULE__)
register_normal_element("q__card__actions", context = @__MODULE__)

function card(args...; kwargs...)
  q__card(args...; kw(kwargs)...)
end

function card_section(args...; kwargs...)
  q__card__section(args...; kw(kwargs)...)
end

const cardsection = card_section

function card_actions(args...; kwargs...)
  q__card__actions(args...; kw(kwargs)...)
end

const cardactions = card_actions

end


text:
# AWS Lambda runtime API
# https://docs.aws.amazon.com/lambda/latest/dg/runtimes-api.html

module AWSLambdaRuntimeInterfaceClient

import HTTP
import JSON

function _init()
    global AWS_LAMBDA_RUNTIME_API = ENV["AWS_LAMBDA_RUNTIME_API"]
    global BASE_URL = "http://$(AWS_LAMBDA_RUNTIME_API)/2018-06-01/runtime"
end

function next_invocation()
    res = HTTP.get("$(BASE_URL)/invocation/next")
    body = JSON.parse(String(res.body))
    headers = Dict(res.headers)
    return headers, body
end

function invocation_response(request_id, response)
    url = "$(BASE_URL)/invocation/$(request_id)/response"
    headers = ["Content-type" => "application/json"]
    HTTP.post(url, headers, response)
end

function invocation_error(request_id)
    url = "$(BASE_URL)/invocation/$(request_id)/error"
    body = Dict("errorType" => "InvalidEventDataException",
                "errorMessage" => "Error parsing event data.")
    headers = Dict("Content-type" => "application/json",
                   "Lambda-Runtime-Function-Error-Type" => "Unhandled")
    HTTP.post(url, headers, body)
end

function initialization_error()
    body = Dict("errorMessage" => "Failed to load function.",
                "errorType" => "InvalidFunctionException")
    headers = Dict("Content-type" => "application/json",
                   "Lambda-Runtime-Function-Error-Type" => "Unhandled")
    url = "$(BASE_URL)/init/error"
    HTTP.post(url, headers, body)
end

function start(handler)
    _init()
    while true
        local state, request_id
        try
            state = :start
            headers, body = next_invocation()
            request_id = headers["Lambda-Runtime-Aws-Request-Id"]
            state = :received

            response  = handler(body, headers)
            state = :handled

            invocation_response(request_id, response)
            state = :finished
        catch ex
            if state == :start
                @error "Unable to receive request" ex
            elseif state == :received
                @error "Unable to handle event" ex request_id
                invocation_error(request_id)
            elseif state == :handled
                @error "Unable to notify lambda runtime (invocation response)" ex request_id
            else
                @error "Unknown state" ex state
            end
        end
    end
end

end # module