text:
module Test_sklf

using LinearAlgebra
using MatrixPencils
using Test

@testset "Structured Matrix Pencils Utilities" begin

@testset "sreduceBF" begin

fast = true; Ty = Float64; Ty = Complex{Float64}
for fast in (true, false)

for Ty in (Float64, Complex{Float64})

n2 = 3; m2 = 1; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = zeros(Ty,n2,n2); B2 = zeros(Ty,n2); C2 = zeros(Ty,p2,n2); D2 = zeros(Ty,p2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,C,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 2; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = zeros(Ty,n2,n2); B2 = zeros(Ty,n2,m2); C2 = zeros(Ty,p2,n2); D2 = zeros(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,C,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 2; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,C,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 0; m2 = 2; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,C,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

      #
n2 = 0; m2 = 2; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(missing,missing,missing,missing,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

@time M, N, Q, Z, n, m, p = sreduceBF(missing,missing,missing,missing,missing,fast = fast)
@test n == 0 && m == 0 && p == 0

n2 = 3; m2 = 0; p2 = 4;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,missing,C,missing,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 3; p2 = 0;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,missing,missing,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 0; p2 = 0;
A2 = rand(Ty,n2,n2); E2 = triu(rand(Ty,n2,n2),1); B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,missing,missing,missing,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == rank(E) && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 0; p2 = 0;
A2 = rand(Ty,n2,n2); E2 = I; B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,missing,missing,missing,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == n2 && m == n2-n+m2 && p == n2-n+p2

n2 = 3; m2 = 3; p2 = 2;
A2 = rand(Ty,n2,n2); E2 = I; B2 = rand(Ty,n2,m2); C2 = rand(Ty,p2,n2); D2 = rand(Ty,p2,m2);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2); D = copy(D2);
M2 = [A2 B2; C2 D2]
N2 = [E zeros(Ty,n2,m2); zeros(Ty,p2,n2+m2)]

@time M, N, Q, Z, n, m, p = sreduceBF(A,E,B,C,D,fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      n == n2 && m == n2-n+m2 && p == n2-n+p2

end
end
end

@testset "sklf, sklf_right and sklf_left" begin

fast = true; Ty = Float64; Ty = Complex{Float64}
for fast in (true, false)
    # test example for SLICOT subroutine AB08ND
A = [
   1.0   0.0   0.0   0.0   0.0   0.0
   0.0   1.0   0.0   0.0   0.0   0.0
   0.0   0.0   3.0   0.0   0.0   0.0
   0.0   0.0   0.0  -4.0   0.0   0.0
   0.0   0.0   0.0   0.0  -1.0   0.0
   0.0   0.0   0.0   0.0   0.0   3.0];
E = I;
B = [
   0.0  -1.0
  -1.0   0.0
   1.0  -1.0
   0.0   0.0
   0.0   1.0
  -1.0  -1.0];
C = [
   1.0   0.0   0.0   1.0   0.0   0.0
   0.0   1.0   0.0   1.0   0.0   1.0
   0.0   0.0   1.0   0.0   0.0   1.0];
D = [
   0.0   0.0
   0.0   0.0
   0.0   0.0];
M2 = [A B; C D]; N2 = [E zeros(6,2); zeros(3,8)];

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [2, 2] && νl == [1, 1, 1] && μl == [0, 1, 1] && nf == 2

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = false, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [2, 2] && νl == [1, 1, 1] && μl == [0, 1, 1] && nf == 2

M2 = [A' B; C D]; N2 = [E zeros(6,2); zeros(3,8)];

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A', E, B, C, D, fast = fast, finite_infinite = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [2, 2] && νl == [1, 1, 1] && μl == [0, 1, 1] && nf == 2

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [2, 2] && νl == [1, 1, 1] && μl == [0, 1, 1] && nf == 2

M2 = [A; C]
N2 = [E; zeros(3,6)]
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, missing, C, missing, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [2, 3, 3] && μl == [0, 2, 3] && nf == 1

M2 = [A B]
N2 = [E zeros(6,2)]
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, missing, missing, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 2, 1, 0] && μr == [2, 2, 2, 1] && νi == [] && νl == [] && μl == [] && nf == 1

M2 = [A B; C D]; N2 = [E zeros(6,2); zeros(3,8)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = sklf_right(A, E, B, C, D, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && ν == [1, 3, 3] && μ == [0, 3, 3] && nf == 2

M2 = [A B; C D]; N2 = [E zeros(6,2); zeros(3,8)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, ν, μ, nf, νl, μl = sklf_left(A, E, B, C, D, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [2, 2] && μ == [2, 2] && νl == [1, 1, 1] && μl == [0, 1, 1] && nf == 2

# test example for SLICOT subroutine AB08ND
A = [
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E = [
    0     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B = [
    -1     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C = [
    0     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];
D = [
    1     2    -2
    0    -1    -2
    0     0     0];
M2 = [A B; C D]; N2 = [E zeros(9,3); zeros(3,12)];

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1 ] && νi == [1, 1, 5] && νl == [1, 1] && μl == [0, 1] && nf == 1

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = false, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1 ] && νi == [1, 1, 5] && νl == [1, 1] && μl == [0, 1] && nf == 1

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1 ] && νi == [5, 1, 1] && νl == [1, 1] && μl == [0, 1] && nf == 1

M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, C, D, fast = fast, finite_infinite = true, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1 ] && νi == [5, 1, 1] && νl == [1, 1] && μl == [0, 1] && nf == 1

M2 = [A; C]
N2 = [E; zeros(3,9)]
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, missing, C, missing, fast = fast, finite_infinite = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [2, 2, 3] && νl == [2, 3] && μl == [0, 2] && nf == 0

M2 = [A; C]
N2 = [E; zeros(3,9)]
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, missing, C, missing, fast = fast, finite_infinite = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [3, 2, 2] && νl == [2, 3] && μl == [0, 2] && nf == 0

M2 = [A B]
N2 = [E zeros(9,3)]
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = sklf(A, E, B, missing, missing, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [3, 3, 0] && μr == [3, 3, 3] && νi == [3] && νl == [] && μl == [] && nf == 0

M2 = [A B; C D]; N2 = [E zeros(9,3); zeros(3,12)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = sklf_right(A, E, B, C, D, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1] && ν == [1, 2, 6] && μ == [1, 1, 6] && nf == 1

M2 = [A B; C D]; N2 = [E zeros(9,3); zeros(3,12)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = sklf_right(A, E, B, C, D, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 0] && μr == [1, 1, 1] && ν == [1, 2, 6] && μ == [1, 1, 6] && nf == 1

M2 = [A B; C D]; N2 = [E zeros(9,3); zeros(3,12)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, ν, μ, nf, νl, μl = sklf_left(A, E, B, C, D, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [6, 2, 1] && μ == [6, 2, 2] && νl == [1, 1] && μl == [0, 1] && nf == 1

M2 = [A B; C D]; N2 = [E zeros(9,3); zeros(3,12)];
M = copy(M2); N = copy(N2);
M1, N1, Q1, Z1, ν, μ, nf, νl, μl = sklf_left(A, E, B, C, D, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [6, 2, 1] && μ == [6, 2, 2] && νl == [1, 1] && μl == [0, 1] && nf == 1

end
end

@testset "sklf_left! - generalized case" begin

for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nfuo == 0 && niuo == 0

A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(3,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [0] && no == 0 && nfuo == 0 && niuo == 0

A2 = rand(3,3); E2 = zeros(3,3); C2 = zeros(0,3); B2 = rand(3);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nfuo == 0 && niuo == 3

A2 = rand(1,1); E2 = rand(1,1); C2 = rand(1,1); B2 = rand(1,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1] && no == 1 && nfuo == 0 && niuo == 0

A2 = rand(1,1); E2 = rand(1,1); C2 = zeros(1,1); B2 = rand(1,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nfuo == 1 && niuo == 0

# Ty = Complex{Float64}; fast = true
# Ty = Float64; fast = true

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,1,3); B2 = rand(Ty,3);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 1, 1] && no == 3 && nfuo == 0 && niuo == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = zeros(Ty,1,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nfuo == 3 && niuo == 0

 A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,2,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nfuo == 0 && niuo == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,4,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [3] && no == 3 && nfuo == 0 && niuo == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = [rand(Ty,4,2) zeros(Ty,4,1)]; B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nfuo == 0 && niuo == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
B2 = rand(Ty,7,4); C2 = [ zeros(Ty,2,3) rand(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [2, 2] && no == 4 && nfuo == 3 && niuo == 0

A2 = [rand(Ty,3,3) zeros(Ty,3,4) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,3) zeros(Ty,3,4); rand(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = rand(Ty,7,4); C2 = [ rand(Ty,2,3) zeros(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nfuo == 3 && niuo == 1

A2 = [rand(Ty,3,3) zeros(Ty,3,4) ; zeros(Ty,4,3) triu(rand(Ty,4,4)) ];
E2 = [rand(Ty,3,3) zeros(Ty,3,4); rand(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = rand(Ty,7,4); C2 = [ rand(Ty,2,3) zeros(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nfuo, niuo  = sklf_left!(A,E,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nfuo == 0 && niuo == 4

end
end
end

@testset "sklf_right! - generalized case" begin

for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,0); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nfu == 0 && niu == 0

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,3); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [0] && nc == 0 && nfu == 0 && niu == 0

A2 = zeros(3,3); E2 = rand(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nfu == 3 && niu == 0

A2 = rand(3,3); E2 = zeros(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nfu == 0 && niu == 3

A2 = rand(3,3); E2 = rand(3,3); B2 = zeros(3,1); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nfu == 3 && niu == 0

A2 = rand(3,3); E2 = rand(3,3); B2 = rand(3,1); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nfu == 0 && niu == 0

# Ty = Complex{Float64}
# fast = true; Ty = Float64

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = rand(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1] && nc == 1 && nfu == 0 && niu == 0

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = zeros(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nfu == 0 && niu == 1

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,1); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nfu == 0 && niu == 0

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,2); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 1] && nc == 3 && nfu == 0 && niu == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); B2 = rand(Ty,3,4); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nfu == 0 && niu == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); B2 = [zeros(Ty,1,4); rand(Ty,2,4)]; C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 1]  && nc == 3 && nfu == 0 && niu == 0

Ty = Float64; fast = true;
A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nfu == 3 && niu == 1

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) triu(rand(Ty,4,4)) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nfu == 0 && niu == 4

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nfu, niu  = sklf_right!(A,E,B,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nfu == 4 && niu == 0

end
end
end

@testset "sklf_rightfin! - generalized case" begin

for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,0); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 0

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,3); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [0] && nc == 0 && nuc == 0

A2 = zeros(3,3); E2 = rand(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = zeros(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = rand(3,3); B2 = zeros(3,1); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = rand(3,3); B2 = rand(3,1); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nuc == 0

# Ty = Complex{Float64}
# fast = true; Ty = Float64

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = rand(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1] && nc == 1 && nuc == 0

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = zeros(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 1

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,1); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nuc == 0

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,2); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 1] && nc == 3 && nuc == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); B2 = rand(Ty,3,4); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nuc == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); B2 = [zeros(Ty,1,4); rand(Ty,2,4)]; C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 1]  && nc == 3 && nuc == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nuc == 4

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin!(A,E,B,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nuc == 4

end
end
end

@testset "sklf_leftfin! - generalized case" begin

for fast in (true, false)

Ty = Float64; fast = false
A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nuo == 0

A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(3,0); B2 = missing;
A = copy(A2); E = copy(E2); C = copy(C2); B = missing;

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [0] && no == 0 && nuo  == 0

A2 = rand(3,3); E2 = zeros(3,3); C2 = zeros(0,3); B2 = rand(3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nuo == 3

A2 = rand(1,1); E2 = rand(1,1); C2 = rand(1,1); B2 = rand(1,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1] && no == 1 && nuo  == 0

A2 = rand(1,1); E2 = rand(1,1); C2 = zeros(1,1); B2 = rand(1,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nuo == 1

# Ty = Complex{Float64}; fast = true
# Ty = Float64; fast = true

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,1,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 1, 1] && no == 3 && nuo  == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,2,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nuo  == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = rand(Ty,4,3); B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [3] && no == 3 && nuo  == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); C2 = [rand(Ty,4,2) zeros(Ty,4,1)]; B2 = rand(Ty,3,2);
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nuo  == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
B2 = rand(Ty,7,4); C2 = [ zeros(Ty,2,3) rand(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [2, 2] && no == 4 && nuo == 3

A2 = [rand(Ty,3,3) zeros(Ty,3,4) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,3) zeros(Ty,3,4); rand(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = rand(Ty,7,4); C2 = [ rand(Ty,2,3) zeros(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(A,E,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [2, 2] && no == 4 && nuo == 3

A2 = [rand(Ty,3,3) zeros(Ty,3,4) ; zeros(Ty,4,3) triu(rand(Ty,4,4)) ];
E2 = [rand(Ty,3,3) zeros(Ty,3,4); rand(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = rand(Ty,7,4); C2 = [ rand(Ty,2,3) zeros(Ty,2,4)];
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, μl, no, nuo  = sklf_leftfin!(E,A,C,B,fast=fast,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(C2*Z-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nuo == 4

end
end
end

@testset "sklf_left! - standard case" begin

for fast in (true, false)

A2 = zeros(0,0); C2 = zeros(0,0); B2 = missing;
A = copy(A2); C = copy(C2); B = missing;

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nu == 0

A2 = zeros(0,0); C2 = zeros(3,0); B2 = missing;
A = copy(A2); C = copy(C2); B = missing;

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [0] && no == 0 && nu == 0

A2 = zeros(3,3); C2 = zeros(0,3); B2 = rand(3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nu == 3

Ty = Complex{Float64}
for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); C2 = rand(Ty,1,3); B2 = rand(Ty,3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 1, 1] && no == 3 && nu == 0

A2 = rand(Ty,3,3); C2 = zeros(Ty,1,3); B2 = rand(Ty,3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [] && no == 0 && nu == 3

A2 = rand(Ty,3,3); C2 = rand(Ty,2,3); B2 = rand(Ty,3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nu == 0

A2 = rand(Ty,3,3); C2 = rand(Ty,4,3); B2 = rand(Ty,3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [3] && no == 3 && nu == 0

A2 = rand(Ty,3,3); C2 = [rand(Ty,4,2) zeros(Ty,4,1)]; B2 = rand(Ty,3,2);
A = copy(A2); C = copy(C2); B = copy(B2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(C2*Q-C) < sqrt(eps(1.)) &&
      (ismissing(B) || norm(Q'*B2-B) < sqrt(eps(1.))) &&
      μl == [1, 2] && no == 3 && nu == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ]; B2 = rand(Ty,7,4); C2 = [ zeros(Ty,2,3) rand(Ty,2,4)];
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, μl, no, nu  = sklf_left!(A,C,B,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      μl == [2, 2] && no == 4 && nu == 3

end
end
end

@testset "sklf_right! - standard case" begin

for fast in (true, false)

A2 = zeros(0,0); B2 = zeros(0,0); C2 = missing;
A = copy(A2); B = copy(B2); C = missing;

@time Q, νr, nc, nu  = sklf_right!(A,B,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nu == 0

A2 = zeros(0,0); B2 = zeros(0,3);C2 = missing;
A = copy(A2); B = copy(B2); C = missing;

@time Q, νr, nc, nu  = sklf_right!(A,B,missing,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [0] && nc == 0 && nu == 0

A2 = zeros(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nu == 3

Ty = Complex{Float64}
for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); B2 = rand(Ty,3,1); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nu == 0

A2 = rand(Ty,3,3); B2 = rand(Ty,3,2); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [2, 1] && nc == 3 && nu == 0

A2 = rand(Ty,3,3); B2 = rand(Ty,3,4); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nu == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ]; B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [3] && nc == 3 && nu == 4

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ]; B2 = [rand(Ty,3,1); zeros(Ty,4,1)]; C2 = rand(Ty,2,7);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right!(A,B,C)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1] && nc == 3 && nu == 4
end
end

end

@testset "sklf_right! and sklf_left! - pencil case" begin

fast = true; Ty = Float64; Ty = Complex{Float64}
for fast in (true, false)

# ensuring strong observability
# Example 4: Van Dooren, Dewilde, LAA 1983
# P = zeros(3,3,3)
# P[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0]
# P[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2]
# P[:,:,3] = [1 4 2; 0 0 0; 1 4 2]

# use a strongly controllable pencil realization

sys1  = ([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[1.0 3.0 0.0; 1.0 4.0 2.0; 0.0 -1.0 -2.0],
[-1.0 -4.0 -2.0; -0.0 -0.0 -0.0; -1.0 -4.0 -2.0],
[1.0 2.0 -2.0; 0.0 -1.0 -2.0; 0.0 0.0 0.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])

# the initial realization is already strongly controllable
(A, E, B, F, C, G, D, H) = copy.(sys1)
Q, Z, nc = sklf_right!(A, E, B, F, C, G, D, H, fast = fast)
i1 = 1:nc
W = Z[1:3,1:3]
@test Q'*[sys1[3] sys1[1]]*Z ≈ [B A] &&
      Q'*[sys1[4] sys1[2]]*Z ≈ [F E] &&
      [sys1[8] sys1[6]]*Z ≈ [H G] &&
      [sys1[7] sys1[5]]*Z ≈ [D C] && nc == size(A,1) &&
      lpsequal(A[i1,i1], E[i1,i1], B[i1,:]/W, F[i1,:]/W, C[:,i1], G[:,i1], D/W, H/W, sys1...,atol1=1.e-7,atol2=1.e-7)

# determine a strongly controllable and strongly observable realization
(A, E, B, F, C, G, D, H) = copy.(sys1)
Q, Z, no = sklf_left!(A, E, C, G, B, F, D, H, fast = fast)
p,n = size(C)
i1 = n-no+1:n
V = Q[end-p+1:end,end-p+1:end]'
@test Q'*[sys1[1]; sys1[5]]*Z ≈ [A;C] &&
      Q'*[sys1[2]; sys1[6]]*Z ≈ [E; G] &&
      Q'*[sys1[4]; sys1[8]] ≈ [F; H] &&
      Q'*[sys1[3]; sys1[7]] ≈ [B; D] &&
      lpsequal(A[i1,i1], E[i1,i1], B[i1,:], F[i1,:], V\C[:,i1], V\G[:,i1], V\D, V\H, sys1...,atol1=1.e-7,atol2=1.e-7)

# the eigenvalues and Kronecker structure of P(s) are the expected ones
Mo = [A[i1,i1] B[i1,:]; C[:,i1] D]
No = [E[i1,i1] F[i1,:]; G[:,i1] H]
val, info = peigvals(Mo,No)
@test val ≈ [1, Inf] && (info.rki, info.lki,info.id, info.nf) == ([0], [1], [1], 1)

# use a strongly observable pencil realizations

sys2  = ([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 3.0 0.0; 1.0 4.0 2.0; 0.0 -1.0 -2.0],
[-1.0 -4.0 -2.0; -0.0 -0.0 -0.0; -1.0 -4.0 -2.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[1.0 2.0 -2.0; 0.0 -1.0 -2.0; 0.0 0.0 0.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])

# the initial realization is already strongly observable
(A, E, B, F, C, G, D, H) = copy.(sys2)
V, Z, no = sklf_left!(A, E, C, G, B, F, D, H, fast = fast, withQ = false, withZ = false)
n = size(A,1)
i1 = n-no+1:n
@test n == no && lpsequal(A[i1,i1], E[i1,i1], B[i1,:], F[i1,:], V\C[:,i1], V\G[:,i1], V\D, V\H, sys2...,atol1=1.e-7,atol2=1.e-7)

# compute a strongly controllable realization
(A, E, B, F, C, G, D, H) = copy.(sys2)
Q, Z, nc = sklf_right!(A, E, B, F, C, G, D, H, fast = fast)
i1 = 1:nc
W = Z[1:3,1:3]
@test Q'*[sys2[3] sys2[1]]*Z ≈ [B A] &&
      Q'*[sys2[4] sys2[2]]*Z ≈ [F E] &&
      [sys2[8] sys2[6]]*Z ≈ [H G] &&
      [sys2[7] sys2[5]]*Z ≈ [D C] &&
      lpsequal(A[i1,i1], E[i1,i1], B[i1,:]/W, F[i1,:]/W, C[:,i1], G[:,i1], D/W, H/W, sys2...,atol1=1.e-7,atol2=1.e-7)

# the eigenvalues and Kronecker structure of P(s) are the expected ones
Mc = [A[i1,i1] B[i1,:]; C[:,i1] D]
Nc = [E[i1,i1] F[i1,:]; G[:,i1] H]
val, info = peigvals(Mc,Nc,atol1=1.e-7,atol2=1.e-7)
@test val[1] ≈ 1 && val[2] == Inf && (info.rki, info.lki,info.id, info.nf) == ([0], [1], [1], 1)

# use a minimal descriptor realization (which is however not strongly minimal)
(A, E, B, C, D) =
([0.2413793103448278 -1.2772592536173888 -0.18569533817705197; -0.04561640191490692 0.2413793103448277 -0.9826073688810348; -0.9826073688810347 -0.18569533817705186 0.0],
[1.0 -4.5102810375397046e-17 0.0; 0.0 -1.0 0.0; 0.0 0.0 0.0],
[0.9191450300180584 1.8382900600361163 -1.838290060036116; -0.17370208344491272 -0.34740416688982534 0.34740416688982545; -1.8708286933869704 -7.483314773547881 -3.741657386773941],
[-0.0656532164298613 0.34740416688982556 0.7071067811865471; -0.1313064328597228 0.6948083337796513 4.163336342344337e-17; 0.06565321642986145 -0.3474041668898257 0.7071067811865477],
[0.7500000000000001 1.5000000000000004 -1.5; -0.5 -1.9999999999999996 -0.9999999999999996; 0.24999999999999997 0.4999999999999998 -0.5000000000000001])
F = zeros(size(B)...)
G = zeros(size(C)...)
H = zeros(size(D)...)
sys = (A, E, B, F, C, G, D, H)
sys2 = copy.(sys)

# this realization is not strongly minimal
# the computed eigenvalues contains two spurious infinite eigenvalues
val, info = peigvals([A B;C D],[E F;G H])
@test val ≈ [1, Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf) == ([0], [1], [1, 1, 1], 1)

# compute a strongly controllable pencil realization
#(A, E, B, F, C, G, D, H) = copy.(sys2)
Q, Z, nc = sklf_right!(A, E, B, F, C, G, D, H, fast = fast)
i1 = 1:nc
W = Z[1:3,1:3]
@test Q'*[sys2[3] sys2[1]]*Z ≈ [B A] &&
      Q'*[sys2[4] sys2[2]]*Z ≈ [F E] &&
      [sys2[8] sys2[6]]*Z ≈ [H G] &&
      [sys2[7] sys2[5]]*Z ≈ [D C] &&
      lpsequal(A[i1,i1], E[i1,i1], B[i1,:]/W, F[i1,:]/W, C[:,i1], G[:,i1], D/W, H/W, sys2...,atol1=1.e-7,atol2=1.e-7)

# there is still a spurious infinite eigenvalue
Mc = [A[i1,i1] B[i1,:]; C[:,i1] D]
Nc = [E[i1,i1] F[i1,:]; G[:,i1] H]
val, info = peigvals(Mc,Nc)
@test val ≈ [1, Inf, Inf] && (info.rki, info.lki,info.id, info.nf) == ([0], [1], [1, 1], 1)

# compute a strongly observable pencil realization
(A, E, B, F, C, G, D, H) = copy.(sys2)
Q, Z, no = sklf_left!(A, E, C, G, B, F, D, H, fast = fast)
p, n = size(C)
i1 = n-no+1:n
V = Q[end-p+1:end,end-p+1:end]'
@test Q'*[sys2[1]; sys2[5]]*Z ≈ [A;C] &&
      Q'*[sys2[2]; sys2[6]]*Z ≈ [E; G] &&
      Q'*[sys2[4]; sys2[8]] ≈ [F; H] &&
      Q'*[sys2[3]; sys2[7]] ≈ [B; D] &&
      lpsequal(A[i1,i1], E[i1,i1], B[i1,:], F[i1,:], V\C[:,i1], V\G[:,i1], V\D, V\H, sys2...,atol1=1.e-7,atol2=1.e-7)

# there is still a spurious infinite eigenvalue
Mo = [A[i1,i1] B[i1,:]; C[:,i1] D]
No = [E[i1,i1] F[i1,:]; G[:,i1] H]
val, info = peigvals(Mo,No)
@test val ≈ [1, Inf, Inf] && (info.rki, info.lki,info.id, info.nf) == ([0], [1], [1, 1], 1)

end
end

@testset "sklf_right2! - standard case" begin

#fast = true
for fast in (true, false)

A2 = zeros(0,0); B2 = zeros(0,0); C2 = missing;
A = copy(A2); B = copy(B2); C = missing;

@time Q, νr, nc, nu  = sklf_right2!(A,B,0,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == Int[] && nc == 0 && nu == 0

A2 = zeros(0,0); B2 = zeros(0,3);C2 = missing;
A = copy(A2); B = copy(B2); C = missing;

@time Q, νr, nc, nu  = sklf_right2!(A,B,2,missing,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [0, 0] && nc == 0 && nu == 0

A2 = zeros(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,0,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nu == 3

A2 = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B2 = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C2 = Matrix{Float64}(I,5,5);

A = copy(A2); B = copy(B2); C = copy(C2);
@time Q, νr, nc, nu  = sklf_right2!(A,B,0,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr ==  [0, 4, 0, 1] && nc == 5 && nu == 0

A = copy(A2); B = copy(B2); C = copy(C2);
@time Q, νr, nc, nu  = sklf_right2!(A,B,1,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr ==  [1, 3, 1, 0] && nc == 5 && nu == 0

A = copy(A2); B = copy(B2); C = copy(C2);
@time Q, νr, nc, nu  = sklf_right2!(A,B,2,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [2, 2, 1, 0] && nc == 5 && nu == 0

A = copy(A2); B = copy(B2); C = copy(C2);
@time Q, νr, nc, nu  = sklf_right2!(A,B,3,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [3, 1, 1, 0] && nc == 5 && nu == 0

A = copy(A2); B = copy(B2); C = copy(C2);
@time Q, νr, nc, nu  = sklf_right2!(A,B,4,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [4, 0, 1, 0] && nc == 5 && nu == 0

#Ty = Complex{Float64}
for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,3,3); B2 = rand(Ty,3,1); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,1,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 0, 1, 0, 1, 0] && nc == 3 && nu == 0

A2 = rand(Ty,3,3); B2 = rand(Ty,3,2); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,1,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1, 0]  && nc == 3 && nu == 0

A2 = rand(Ty,3,3); B2 = rand(Ty,3,4); C2 = rand(Ty,2,3);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,2,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [2, 1] && nc == 3 && nu == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ]; B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,1,C)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 2] && nc == 3 && nu == 4

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ]; B2 = [rand(Ty,3,2); zeros(Ty,4,2)]; C2 = rand(Ty,2,7);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,1,C)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1, 0] && nc == 3 && nu == 4

A2 = rand(Ty,20,20); B2 = [rand(Ty,2,5);zeros(Ty,18,5)]; C2 = rand(Ty,2,20);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,4,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0] && nc == 20 && nu == 0

A2 = rand(Ty,20,20); B2 = rand(Ty,20,5); C2 = rand(Ty,2,20);
A = copy(A2); B = copy(B2); C = copy(C2);

@time Q, νr, nc, nu  = sklf_right2!(A,B,2,C,fast = fast)
@test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
      νr == [2, 3, 2, 3, 2, 3, 2, 3] && nc == 20 && nu == 0

# A2 = [ rand(Ty,5,20); rand(Ty,2,20); [zeros(Ty,13,5) rand(Ty,13,15)]]; B2 = [rand(Ty,5,5);zeros(Ty,15,5)]; C2 = rand(Ty,2,20);
# A = copy(A2); B = copy(B2); C = copy(C2);

# @time Q, νr, nc, nu  = sklf_right2!(A,B,1,C,fast = fast)
# @test norm(Q'*A2*Q-A) < sqrt(eps(1.)) &&
#       norm(Q'*B2-B) < sqrt(eps(1.)) &&
#       (ismissing(C) || norm(C2*Q-C) < sqrt(eps(1.))) &&
#       νr == [2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0] && nc == 20 && nu == 0

end
end
end

@testset "sklf_rightfin2! - generalized case" begin

#fast = true
for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,0); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,0,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 0

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,3); C2 = missing;
A = copy(A2); E = copy(E2); B = copy(B2); C = missing;

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,2,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [0, 0] && nc == 0 && nuc == 0

A2 = zeros(3,3); E2 = rand(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,0,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = zeros(3,3); B2 = zeros(3,0); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,0,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = rand(3,3); B2 = zeros(3,2); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 3

A2 = rand(3,3); E2 = rand(3,3); B2 = rand(3,2); C2 = rand(2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast=fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1, 0] && nc == 3 && nuc == 0

A2 = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B2 = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C2 = Matrix{Float64}(I,5,5);
E2 = Matrix{Float64}(I,5,5);

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time Q, Z, νr, nc, nu  = sklf_rightfin2!(A,E,B,0,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr ==  [0, 4, 0, 1] && nc == 5 && nu == 0

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time Q, Z, νr, nc, nu  = sklf_rightfin2!(A,E,B,1,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr ==  [1, 3, 1, 0]  && nc == 5 && nu == 0

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time Q, Z, νr, nc, nu  = sklf_rightfin2!(A,E,B,2,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 2, 1, 0]  && nc == 5 && nu == 0

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time Q, Z, νr, nc, nu  = sklf_rightfin2!(A,E,B,3,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [3, 1, 1, 0]  && nc == 5 && nu == 0

A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);
@time Q, Z, νr, nc, nu  = sklf_rightfin2!(A,E,B,4,C,atol1 = 1.e-7, fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [4, 0, 1, 0]  && nc == 5 && nu == 0

#fast = false; Ty = Float64; Ty = Complex{Float64}

for Ty in (Float64, Complex{Float64})

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = rand(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 0] && nc == 1 && nuc == 0

A2 = rand(Ty,1,1); E2 = triu(rand(Ty,1,1),1); B2 = zeros(Ty,1,1); C2 = rand(Ty,2,1);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [] && nc == 0 && nuc == 1

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,1); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 0, 1, 0, 1, 0] && nc == 3 && nuc == 0

A2 = rand(Ty,3,3); E2 = triu(rand(Ty,3,3),1); B2 = rand(Ty,3,2); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2); B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 1, 1, 0] && nc == 3 && nuc == 0

A2 = rand(Ty,3,3); E2 = rand(Ty,3,3); B2 = rand(Ty,3,4); C2 = rand(Ty,2,3);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,2,C,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 1] && nc == 3 && nuc == 0

A2 = [rand(Ty,3,7) ; zeros(Ty,4,3) rand(Ty,4,4) ];
E2 = [rand(Ty,3,7) ; zeros(Ty,4,3) triu(rand(Ty,4,4),1) ];
B2 = [rand(Ty,3,4); zeros(Ty,4,4)]; C2 = rand(Ty,2,7);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,1,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [1, 2] && nc == 3 && nuc == 4

A2 = rand(Ty,20,20); E2 = triu(rand(Ty,20,20)); B2 = [rand(Ty,2,5);zeros(Ty,18,5)]; C2 = rand(Ty,2,20);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,4,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0] && nc == 20 && nu == 0

A2 = rand(Ty,20,20); E2 = rand(Ty,20,20); B2 = rand(Ty,20,5); C2 = rand(Ty,2,20);
A = copy(A2); E = copy(E2);  B = copy(B2); C = copy(C2);

@time Q, Z, νr, nc, nuc  = sklf_rightfin2!(A,E,B,2,C,atol1=1.e-7,atol2=1.e-7,fast = fast)
@test norm(Q'*A2*Z-A) < sqrt(eps(1.)) &&
      norm(Q'*E2*Z-E) < sqrt(eps(1.)) &&
      norm(Q'*B2-B) < sqrt(eps(1.)) &&
      (ismissing(C) || norm(C2*Z-C) < sqrt(eps(1.))) &&
      νr == [2, 3, 2, 3, 2, 3, 2, 3] && nc == 20 && nu == 0

end
end
end

end
end


text:
using RecipeUtils
using Test

@testset "RecipeUtils.jl" begin
    # Write your own tests here.
end

ext: jl


text:
function train_layer_1_rbf_deriv(neurons,vars,obs,n_per_part,inds_all,x,y,cc)
    a_all=Vector{Vector{Float64}}()
    ooops=0
    errs=Array{Float64}(undef,0)
    obs1=50
    itmps=Vector{Vector{Int64}}()
    for i=1:neurons
        if i in 1:1:neurons println(Dates.format(now(), "HH:MM:SS")," Calculating weights w for ",i," neuron, of ",neurons," Total.") end
        # try

            # x1=copy(x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:])
            # obs1=size(x1,1)
            # phi1 = calc_phi(x1,x1,obs1,obs1,cc,vars)
            # atmp=phi1\y[inds_all[n_per_part[i]+1:n_per_part[i+1]]]
            # # phi1_deriv = calc_phi_deriv(x1,x1,obs1,obs1,cc,vars,rand(1:vars))
            # # println(sum(abs.(atmp[1:end])))
            # # atmp=[phi1;phi1_deriv]\[y[inds_all[n_per_part[i]+1:n_per_part[i+1]]];zeros(obs1)]

            itmp=randperm(obs)[1:obs1]
            x1=copy(x[itmp,:])
            phi1 = calc_phi(x1,x1,obs1,obs1,cc,vars)
            phi1_deriv = calc_phi_deriv(x1,x1,obs1,obs1,cc,vars,rand(1:vars))
            atmp=[phi1;phi1_deriv]\[y[itmp];zeros(obs1)]
            err1=mean(abs.(phi1*atmp-y[itmp]))

            if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 && sum(abs.(atmp[1:end]))>0
                push!(a_all,atmp)
                push!(itmps,itmp)
                push!(errs,err1)
            else ooops+=1 end
        # catch ex
        #     ooops+=1;
        #     # println(i," ",ex)
        # end
    end
    println("ooops=",ooops)

    ii1=sortperm(errs)[1:2000]
    a_all=a_all[ii1]
    itmps=itmps[ii1]

    return a_all,itmps
end


text:
# WARNING this file is only loaded if DataFrames.jl is also active
# This all happens thanks to the Requires.jl package

import Base: getindex
import Base: setindex!
import SimpleSDMLayers: clip, latitudes, longitudes

"""
    Base.getindex(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}

Returns the values of a layer at all occurrences in a `DataFrame`. Note that the
function looks for columns named `:latitude` and `:longitude` by default, but
these can be changed using the `latitude` and `longitude` arguments.
"""
function Base.getindex(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}
    lats = df[!, latitude]
    lons = df[!, longitude]
    return [layer[lon, lat] for (lon, lat) in zip(lons, lats)]
end

"""
    clip(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}

Returns a clipped version (with a 10% margin) around all occurences in a
`DataFrame`. Note that the function looks for columns named `:latitude` and
`:longitude` by default, but these can be changed using the `latitude` and
`longitude` arguments.
"""
function SimpleSDMLayers.clip(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}
   occ_latitudes = filter(!ismissing, df[!, latitude])
   occ_longitudes = filter(!ismissing, df[!, longitude])

   lat_min = minimum(occ_latitudes)
   lat_max = maximum(occ_latitudes)
   lon_min = minimum(occ_longitudes)
   lon_max = maximum(occ_longitudes)

   lat_Δ = abs(lat_max - lat_min)
   lon_Δ = abs(lon_max - lon_min)

   scaling = 0.1
   lon_s = scaling*lon_Δ
   lat_s = scaling*lat_Δ

   lat_max = min(layer.top, lat_max+lat_s)
   lat_min = max(layer.bottom, lat_min-lat_s)
   lon_max = min(layer.right, lon_max+lon_s)
   lon_min = max(layer.left, lon_min-lon_s)

   return layer[left=lon_min, right=lon_max, bottom=lat_min, top=lat_max]
end

"""
    DataFrames.DataFrame(layer::T) where {T <: SimpleSDMLayer}

Returns a DataFrame from a `SimpleSDMLayer` element, with columns for latitudes,
longitudes and grid values. `nothing` values in the layer grid are replaced by
`missing` values in the DataFrame.
"""
function DataFrames.DataFrame(layer::T; kw...) where {T <: SimpleSDMLayer}
    lats = repeat(latitudes(layer), outer = size(layer, 2))
    lons = repeat(longitudes(layer), inner = size(layer, 1))
    values = replace(vec(layer.grid), nothing => missing)
    df = DataFrames.DataFrame(longitude = lons, latitude = lats, values = values; kw...)
    return df
end

"""
    DataFrames.DataFrame(layers::Array{SimpleSDMLayer})

Returns a single DataFrame from an `Array` of compatible`SimpleSDMLayer`
elements, with every layer as a column, as well as columns for latitudes and
longitudes. `nothing` values in the layer grids are replaced by
`missing` values in the DataFrame.
"""
function DataFrames.DataFrame(layers::Array{T}; kw...) where {T <: SimpleSDMLayer}
    l1 = layers[1]
    all(x -> SimpleSDMLayers._layers_are_compatible(x, l1), layers)

    lats = repeat(latitudes(l1), outer = size(l1, 2))
    lons = repeat(longitudes(l1), inner = size(l1, 1))
    values = mapreduce(x -> vec(x.grid), hcat, layers)
    values = replace(values, nothing => missing)

    df = DataFrames.DataFrame(values, :auto; kw...)
    DataFrames.insertcols!(df, 1, :longitude => lons)
    DataFrames.insertcols!(df, 2, :latitude => lats)
    return df
end

for ty in (:SimpleSDMResponse, :SimpleSDMPredictor)
    eval(
        quote
            """
                $($ty)(df::DataFrame, col::Symbol, layer::T; latitude::Symbol = :latitude, longitude::Symbol = :longitude) where {T <: SimpleSDMLayer}

            Returns a `$($ty)` from a `DataFrame`. This requires to select a
            column to be returned as a layer, as well as an existing layer from
            which to copy the dimensions. Note that the function looks for
            columns named `:latitude` and `:longitude` by default, but these can
            changed using the `latitude` and `longitude` arguments.
            """
            function SimpleSDMLayers.$ty(df::DataFrames.DataFrame, col::Symbol, layer::SimpleSDMLayer; latitude::Symbol = :latitude, longitude::Symbol = :longitude)
                lats = df[!, latitude]
                lons = df[!, longitude]

                uniquelats = unique(lats)
                uniquelons = unique(lons)

                grid = Array{Any}(nothing, size(layer))

                if uniquelats == latitudes(layer) && uniquelons == longitudes(layer)
                    grid[:] = df[:, col]
                else
                    lats_idx = [SimpleSDMLayers._match_latitude(layer, lat) for lat in lats]
                    lons_idx = [SimpleSDMLayers._match_longitude(layer, lon) for lon in lons]
                    for (lat, lon, value) in zip(lats_idx, lons_idx, df[:, col])
                        grid[lat, lon] = value
                    end
                end

                internal_types = unique(typeof.(grid))
                return SimpleSDMLayers.$ty(Array{Union{internal_types...}}(grid), layer)
            end
        end,
    )
end

"""
    mask!(layer::SimpleSDMResponse{T}, records::DataFrames.DataFrame) where {T <: AbstractBool}

Fills a layer (most likely created with `similar`) so that the values are `true`
if an occurrence is found in the cell, `false` if not. Note that the function
looks for columns named `:latitude` and `:longitude` by default, but these can
be changed using the `latitude` and `longitude` arguments.
"""
function mask!(layer::SimpleSDMResponse{T}, df::DataFrames.DataFrame; latitude::Symbol = :latitude, longitude::Symbol = :longitude) where {T <: Bool}
    uniquedf = unique(df, [longitude, latitude])
    lons = uniquedf[!, longitude]
    lats = uniquedf[!, latitude]
    for (lon, lat) in zip(lons, lats)
        layer[lon, lat] = true
    end
    return layer
end

"""
    mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Number}

Fills a layer (most likely created with `similar`) so that the values reflect
the number of occurrences in the cell. Note that the function looks for columns
named `:latitude` and `:longitude` by default, but these can be changed using
the `latitude` and `longitude` arguments.
"""
function mask!(layer::SimpleSDMResponse{T}, df::DataFrames.DataFrame; latitude::Symbol = :latitude, longitude::Symbol = :longitude) where {T <: Number}
    lons = df[!, longitude]
    lats = df[!, latitude]
    for (lon, lat) in zip(lons, lats)
        layer[lon, lat] = layer[lon, lat] + one(T)
    end
    return layer
end

"""
    mask(layer::SimpleSDMLayer, records::GBIF.GBIFRecords, element_type::Type=Bool)

Create a new layer storing information about the presence of occurrences in the
cells, either counting (numeric types) or presence-absence-ing (boolean types)
them. Note that the function looks for columns named `:latitude` and
`:longitude` by default, but these can be changed using the `latitude` and
`longitude` arguments.
"""
function mask(layer::SimpleSDMLayer, df::DataFrames.DataFrame, element_type::Type=Bool; latitude::Symbol = :latitude, longitude::Symbol = :longitude)
    returnlayer = similar(layer, element_type)
    mask!(returnlayer, df, latitude = latitude, longitude = longitude)
    return returnlayer
end


text:
#############################################################################
# modularity.jl
# This file handles the basic rules on interactions of mathematical expressions
# to create new expressions concerning their modularity.
#
# For example: negative of a SubModularity expression is SuperModularity.
#
# It's not an all-inclusive file, and set functions whose modularity is not determined
# are classified as NotDetermined.
#############################################################################

import Base: +, -, *, sign

export ExprModularity, SubModularity, SuperModularity, Modularity, ConstModularity, NotDetermined
export +, -, *

# ExprModularity subtypes, abstract types
abstract type ExprModularity		        end
struct SubModularity <: ExprModularity		end
struct SuperModularity <: ExprModularity	end
struct Modularity <: ExprModularity       	end
struct ConstModularity <: ExprModularity  	end

struct NotDetermined <: ExprModularity
	function NotDetermined()
		warn("Expression's modularity is not determined. Trying to solve problems with uncertain modularity can lead to unexpected behavior.")
    return new()
	end
end

+(m::ExprModularity) = m

+(m::NotDetermined, n::ExprModularity) = m
+(m::ExprModularity, n::NotDetermined) = n
+(m::NotDetermined, n::NotDetermined) = m

+(m::ConstModularity, n::ConstModularity) = m
+(m::ConstModularity, n::SubModularity) = n
+(m::SubModularity, n::ConstModularity) = m
+(m::ConstModularity, n::SuperModularity) = n
+(m::SuperModularity, n::ConstModularity) = m

+(m::Modularity, n::Modularity) = m
+(m::Modularity, n::SubModularity) = n
+(m::SubModularity, n::Modularity) = m
+(m::Modularity, n::SuperModularity) = NotDetermined()
+(m::SuperModularity, n::Modularity) = NotDetermined()

-(m::ExprModularity) = m
-(m::SubModularity) = SuperModularity()
-(m::SuperModularity) = SubModularity()

-(m::ExprModularity, n::ExprModularity) = m + (-n)


text:
using QMatrices
using Test

include("test_basic.jl")


text:
using Revise
using BosonSampling
using Permanents
using PrettyTables
using ArgCheck
using Test
#
# using DocumenterTools


text:
import FLoops

let sumwith, useinit, smoketest

    function sumwith(f, xs, ex = nothing)
        FLoops.@floop ex for x in xs
            FLoops.@reduce s += f(x)
        end
        return s
    end

    function useinit(xs, ex = nothing)
        T = eltype(xs)
        FLoops.@floop ex for x in xs
            FLoops.@init v = zeros(T, 1)
            v[1] = x
            FLoops.@reduce s += v[1]
        end
        return s
    end

    function smoketest(ex = nothing)
        sumwith(identity, 1:10, ex)
    end

    (
        sumwith = sumwith,
        useinit = useinit,
        smoketest = smoketest,  # just run all functions
    )
end


text:
"""

$(TYPEDEF)

Simple structure to contain the number of samples of each type
of calculation to compute final results

$(TYPEDFIELDS)

"""
@with_kw struct Samples
    md::Float64
    random::Int
end


text:
import Base: ==
using Serialization
import Serialization: serialize, deserialize

###### Array Domains ######

struct ArrayDomain{N}
    indexes::NTuple{N, Any}
end

include("../lib/domain-blocks.jl")

ArrayDomain(xs...) = ArrayDomain(xs)
ArrayDomain(xs::Array) = ArrayDomain((xs...,))

indexes(a::ArrayDomain) = a.indexes
chunks(a::ArrayDomain{N}) where {N} = DomainBlocks(
    ntuple(i->first(indexes(a)[i]), Val(N)), map(x->[length(x)], indexes(a)))

(==)(a::ArrayDomain, b::ArrayDomain) = indexes(a) == indexes(b)
Base.getindex(arr::AbstractArray, d::ArrayDomain) = arr[indexes(d)...]

function intersect(a::ArrayDomain, b::ArrayDomain)
    if a === b
        return a
    end
    ArrayDomain(map((x, y) -> _intersect(x, y), indexes(a), indexes(b)))
end

function project(a::ArrayDomain, b::ArrayDomain)
    map(indexes(a), indexes(b)) do p, q
        q .- (first(p) - 1)
    end |> ArrayDomain
end

function getindex(a::ArrayDomain, b::ArrayDomain)
    ArrayDomain(map(getindex, indexes(a), indexes(b)))
end

"""
    alignfirst(a)

Make a subdomain a standalone domain. For example,

    alignfirst(ArrayDomain(11:25, 21:100))
    # => ArrayDomain((1:15), (1:80))
"""
alignfirst(a::ArrayDomain) =
    ArrayDomain(map(r->1:length(r), indexes(a)))

function size(a::ArrayDomain, dim)
    idxs = indexes(a)
    length(idxs) < dim ? 1 : length(idxs[dim])
end
size(a::ArrayDomain) = map(length, indexes(a))
length(a::ArrayDomain) = prod(size(a))
ndims(a::ArrayDomain) = length(size(a))
isempty(a::ArrayDomain) = length(a) == 0

"The domain of an array is a ArrayDomain"
domain(x::AbstractArray) = ArrayDomain([1:l for l in size(x)])

abstract type ArrayOp{T, N} <: AbstractArray{T, N} end
Base.IndexStyle(::Type{<:ArrayOp}) = IndexCartesian()

compute(ctx, x::ArrayOp) =
    compute(ctx, cached_stage(ctx, x)::DArray)

collect(ctx::Context, x::ArrayOp) =
    collect(ctx, compute(ctx, x))

collect(x::ArrayOp) = collect(Context(), x)

function Base.show(io::IO, ::MIME"text/plain", x::ArrayOp)
    write(io, string(typeof(x)))
    write(io, string(size(x)))
end

function Base.show(io::IO, x::ArrayOp)
    m = MIME"text/plain"()
    show(io, m, x)
end

"""
`DArray{T,N,F}(domain, subdomains, chunks, concat)`

An N-dimensional distributed array of element type T.

- `domain`: the whole ArrayDomain of the array
- `subdomains`: a `DomainBlocks` of the same dimensions as the array
- `chunks`: an array of chunks of dimension N
- `concat`: a function of type `F`. `concat(x, y; dims=d)` takes two chunks `x` and `y`
            and concatenates them along dimension `d`. `cat` is used by default.
"""
mutable struct DArray{T,N,F} <: ArrayOp{T, N}
    domain::ArrayDomain{N}
    subdomains::AbstractArray{ArrayDomain{N}, N}
    chunks::AbstractArray{Union{Chunk,Thunk}, N}
    concat::F
    freed::Threads.Atomic{UInt8}
    function DArray{T,N,F}(domain, subdomains, chunks, concat) where {T, N,F}
        A = new(domain, subdomains, chunks, concat, Threads.Atomic{UInt8}(0))
        refcount_chunks(A.chunks)
        finalizer(free!, A)
        A
    end
end

function serialize(io::AbstractSerializer, A::DArray)
    @async refcount_chunks(A)
    invoke(serialize, Tuple{AbstractSerializer,Any}, io, A)
end

function deserialize(io::AbstractSerializer, dt::Type{DArray{T,N,F}}) where {T,N,F}
    nf = fieldcount(dt)
    A = ccall(:jl_new_struct_uninit, Any, (Any,), dt)
    Serialization.deserialize_cycle(io, A)
    for i in 1:nf
        tag = Int32(read(io.io, UInt8)::UInt8)
        if tag != Serialization.UNDEFREF_TAG
            ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), A, i-1, Serialization.handle_deserialize(io, tag))
        end
    end
    finalizer(free!, A)
    A
end

refcount_chunks(A::DArray) = refcount_chunks(A.chunks)
function refcount_chunks(chunks)
    for c in chunks
        if c isa Chunk{<:Any, DRef}
            # increment refcount on the master node
            addrefcount(c.handle, 1)
        elseif c isa Thunk
            refcount_chunks(c.inputs)
        end
    end
end

function free_chunks(chunks)
    @sync for c in chunks
        if c isa Chunk{<:Any, DRef}
            # increment refcount on the master node
            cnt = addrefcount(c.handle, -1)
            cnt <= 0 && @async free!(c.handle)
        elseif c isa Thunk
            free_chunks(c.inputs)
        end
    end
end

function free!(x::DArray)
    freed = Bool(Threads.atomic_cas!(x.freed, UInt8(0), UInt8(1)))
    !freed && @async Dagger.free_chunks(x.chunks)
    nothing
end

# mainly for backwards-compatibility
DArray{T, N}(domain, subdomains, chunks) where {T,N} = DArray(T, domain, subdomains, chunks)

"""
`DArray(T, domain, subdomains, chunks, [concat=cat])`

Creates a distributed array of element type T.

- `T`: element type

rest of the arguments are the same as the DArray constructor.
"""
function DArray(T, domain::ArrayDomain{N},
             subdomains::AbstractArray{ArrayDomain{N}, N},
             chunks::AbstractArray{<:Any, N}, concat=cat) where N
    DArray{T, N, typeof(concat)}(domain, subdomains, chunks, concat)
end

domain(d::DArray) = d.domain
chunks(d::DArray) = d.chunks
domainchunks(d::DArray) = d.subdomains
size(x::DArray) = size(domain(x))
stage(ctx, c::DArray) = c

function collect(ctx::Context, d::DArray; tree=false)
    a = compute(ctx, d)

    if isempty(d.chunks)
        return Array{eltype(d)}(undef, size(d)...)
    end

    dimcatfuncs = [(x...) -> d.concat(x..., dims=i) for i in 1:ndims(d)]
    if tree
        collect(treereduce_nd(delayed.(dimcatfuncs), a.chunks))
    else
        treereduce_nd(dimcatfuncs, asyncmap(collect, a.chunks))
    end
end

function (==)(x::ArrayOp, y::ArrayOp)
    x === y || reduce((a,b)->a&&b, map(==, x, y))
end

function Base.hash(x::ArrayOp, i::UInt)
    7*objectid(x)-2
end

function Base.isequal(x::ArrayOp, y::ArrayOp)
    x === y
end

"""
`view` of a `DArray` chunk returns a `DArray` of thunks
"""
function Base.view(c::DArray, d)
    subchunks, subdomains = lookup_parts(chunks(c), domainchunks(c), d)
    d1 = alignfirst(d)
    DArray(eltype(c), d1, subdomains, subchunks)
end

function group_indices(cumlength, idxs,at=1, acc=Any[])
    at > length(idxs) && return acc
    f = idxs[at]
    fidx = searchsortedfirst(cumlength, f)
    current_block = (get(cumlength, fidx-1,0)+1):cumlength[fidx]
    start_at = at
    end_at = at
    for i=(at+1):length(idxs)
        if idxs[i] in current_block
            end_at += 1
            at += 1
        else
            break
        end
    end
    push!(acc, fidx=>idxs[start_at:end_at])
    group_indices(cumlength, idxs, at+1, acc)
end

function group_indices(cumlength, idx::Int)
    group_indices(cumlength, [idx])
end

function group_indices(cumlength, idxs::AbstractRange)
    f = searchsortedfirst(cumlength, first(idxs))
    l = searchsortedfirst(cumlength, last(idxs))
    out = cumlength[f:l]
    isempty(out) && return []
    out[end] = last(idxs)
    map(=>, f:l, map(UnitRange, vcat(first(idxs), out[1:end-1].+1), out))
end

_cumsum(x::AbstractArray) = length(x) == 0 ? Int[] : cumsum(x)
function lookup_parts(ps::AbstractArray, subdmns::DomainBlocks{N}, d::ArrayDomain{N}) where N
    groups = map(group_indices, subdmns.cumlength, indexes(d))
    sz = map(length, groups)
    pieces = Array{Union{Chunk,Thunk}}(undef, sz)
    for i = CartesianIndices(sz)
        idx_and_dmn = map(getindex, groups, i.I)
        idx = map(x->x[1], idx_and_dmn)
        dmn = ArrayDomain(map(x->x[2], idx_and_dmn))
        pieces[i] = delayed(getindex)(ps[idx...], project(subdmns[idx...], dmn))
    end
    out_cumlength = map(g->_cumsum(map(x->length(x[2]), g)), groups)
    out_dmn = DomainBlocks(ntuple(x->1,Val(N)), out_cumlength)
    pieces, out_dmn
end

"""
A DArray object may contain a thunk in it, in which case
we first turn it into a Thunk object and then compute it.
"""
function compute(ctx, x::DArray; persist=true)
    thunk = thunkize(ctx, x, persist=persist)
    if isa(thunk, Thunk)
        compute(ctx, thunk)
    else
        x
    end
end

"""
If a DArray tree has a Thunk in it, make the whole thing a big thunk
"""
function thunkize(ctx, c::DArray; persist=true)
    if any(istask, chunks(c))
        thunks = chunks(c)
        sz = size(thunks)
        dmn = domain(c)
        dmnchunks = domainchunks(c)
        if persist
            foreach(persist!, thunks)
        end
        Thunk(thunks...; meta=true) do results...
            t = eltype(results[1])
            DArray(t, dmn, dmnchunks,
                                  reshape(Union{Chunk,Thunk}[results...], sz))
        end
    else
        c
    end
end

global _stage_cache = WeakKeyDict{Context, Dict}()
"""
A memoized version of stage. It is important that the
tasks generated for the same DArray have the same
identity, for example:

    A = rand(Blocks(100,100), Float64, 1000, 1000)
    compute(A+A')

must not result in computation of A twice.
"""
function cached_stage(ctx, x)
    cache = if !haskey(_stage_cache, ctx)
        _stage_cache[ctx] = Dict()
    else
        _stage_cache[ctx]
    end

    if haskey(cache, x)
        cache[x]
    else
        cache[x] = stage(ctx, x)
    end
end

Base.@deprecate_binding Cat DArray
Base.@deprecate_binding ComputedArray DArray

export Distribute, distribute

struct Distribute{T, N} <: ArrayOp{T, N}
    domainchunks
    data::AbstractArray{T,N}
end

size(x::Distribute) = size(domain(x.data))

export BlockPartition, Blocks

struct Blocks{N}
    blocksize::NTuple{N, Int}
end
Blocks(xs::Int...) = Blocks(xs)

Base.@deprecate BlockPartition Blocks

Distribute(p::Blocks, data::AbstractArray) =
    Distribute(partition(p, domain(data)), data)

function stage(ctx, d::Distribute)
    if isa(d.data, ArrayOp)
        # distributing a dsitributed array
        x = cached_stage(ctx, d.data)
        if d.domainchunks == domainchunks(x)
            return x # already properly distributed
        end
        Nd = ndims(x)
        T = eltype(d.data)
        concat = x.concat
        cs = map(d.domainchunks) do idx
            chunks = cached_stage(ctx, x[idx]).chunks
            shape = size(chunks)
            (delayed() do shape, parts...
                if prod(shape) == 0
                    return Array{T}(undef, shape)
                end
                dimcatfuncs = [(x...) -> concat(x..., dims=i) for i in 1:length(shape)]
                ps = reshape(Any[parts...], shape)
                collect(treereduce_nd(dimcatfuncs, ps))
            end)(shape, chunks...)
        end
    else
        cs = map(c -> delayed(identity)(d.data[c]), d.domainchunks)
    end

    DArray(
           eltype(d.data),
           domain(d.data),
           d.domainchunks,
           cs
    )
end

function distribute(x::AbstractArray, dist)
    compute(Distribute(dist, x))
end

function distribute(x::AbstractArray{T,N}, n::NTuple{N}) where {T,N}
    p = map((d, dn)->ceil(Int, d / dn), size(x), n)
    distribute(x, Blocks(p))
end

function distribute(x::AbstractVector, n::Int)
    distribute(x, (n,))
end

function distribute(x::AbstractVector, n::Vector{<:Integer})
    distribute(x, DomainBlocks((1,), (cumsum(n),)))
end

function Base.:(==)(x::ArrayOp{T,N}, y::AbstractArray{S,N}) where {T,S,N}
    collect(x) == y
end

function Base.:(==)(x::AbstractArray{T,N}, y::ArrayOp{S,N}) where {T,S,N}
    return collect(x) == y
end


text:
##################################################
# Compute a map of the NR diffusion coefficients
# In log-log action space
##################################################
using HDF5
########################################
include("../sources/julia/Main.jl") # Loading the main code
########################################

m_field  = _M/nbGlobularCluster # Mass of field star

########################################
# Action space parameter
########################################

JrminMeasure, JrmaxMeasure = 0.00001*_L0,100.0*_L0 # Jr range
LminMeasure, LmaxMeasure = 0.00001*_L0,100.0*_L0 # L range
nbJrMeasure = 150 # Number of Jr sampling points
nbLMeasure = 150 # Number of L sampling points

########################################

tabJrMeasure = exp.(range(log(JrminMeasure),length=nbJrMeasure,log(JrmaxMeasure)))
tabLMeasure = exp.(range(log(LminMeasure),length=nbLMeasure,log(LmaxMeasure)))
nbJrLGrid = nbJrMeasure*nbLMeasure

const tabLJrGrid  = zeros(Float64,2,nbJrLGrid) # Location (L,Jr) of the grid points where the diffusion coefficients are computed
const tabDNRJr  = zeros(Float64,nbJrLGrid) # Values of the DNR_Jr coefficients on the (L,Jr)-grid
const tabDNRJrJr = zeros(Float64,nbJrLGrid) # Values of the DNR_JrJr coefficients on the (L,Jr)-grid
const tabDNRL  = zeros(Float64,nbJrLGrid) # Values of the DNR_L coefficients on the (L,Jr)-grid
const tabDNRLL = zeros(Float64,nbJrLGrid) # Values of the DNR_LL coefficients on the (L,Jr)-grid
const tabDNRJrL = zeros(Float64,nbJrLGrid) # Values of the DNR_JrL coefficients on the (L,Jr)-grid

########################################
# Functions to fill the arrays of coefficients
########################################

function tabLJrGrid!()
    index = 1
    for iJr=1:nbJrMeasure
    JrMeasure = tabJrMeasure[iJr]
        for iL=1:nbLMeasure
            LMeasure = tabLMeasure[iL]
            tabLJrGrid[1,index], tabLJrGrid[2,index] = LMeasure, JrMeasure
            index += 1
        end
    end
end

function tabDNR!()
    if (PARALLEL == "yes")
        println("parallel")
        Threads.@threads for iGrid=1:nbJrLGrid
            LMeasure, JrMeasure = tabLJrGrid[1,iGrid], tabLJrGrid[2,iGrid]

            avrDJr, avrDL, avrDJrJr, avrDJrL, avrDLL = avr_action_coefficients(JrMeasure,LMeasure,m_field)

            tabDNRJr[iGrid] = avrDJr
            tabDNRJrJr[iGrid] = avrDJrJr
            tabDNRL[iGrid] = avrDL
            tabDNRLL[iGrid] = avrDLL
            tabDNRJrL[iGrid] = avrDJrL

        end
    else # Computation is not made in parallel
        for iGrid=1:nbJrLGrid
            LMeasure, JrMeasure = tabLJrGrid[1,iGrid], tabLJrGrid[2,iGrid]

            avrDJr, avrDL, avrDJrJr, avrDJrL, avrDLL = avr_action_coefficients(JrMeasure,LMeasure,m_field)

            tabDNRJr[iGrid] = avrDJr
            tabDNRJrJr[iGrid] = avrDJrJr
            tabDNRL[iGrid] = avrDL
            tabDNRLL[iGrid] = avrDLL
            tabDNRJrL[iGrid] = avrDJrL

        end
    end
end

########################################
namefile = "../data/Dump_Diffusion_Coefficients_Action_Orbital_Map_q_"*string(qCalc)*".hf5"
########################################
# Function that writes to .hf5 files
########################################
function writedump!(namefile)
    file = h5open(namefile,"w") # Opening the file
    write(file,"tabL",tabLMeasure)
    write(file,"q",qCalc)
    write(file,"tabJr",tabJrMeasure)
    write(file,"tabLJr",tabLJrGrid)
    write(file,"tabDNRJr",tabDNRJr)
    write(file,"tabDNRJrJr",tabDNRJrJr)
    write(file,"tabDNRL",tabDNRL)
    write(file,"tabDNRLL",tabDNRLL)
    write(file,"tabDNRJrL",tabDNRJrL)

    write(file,"nbJrMeasure",nbJrMeasure)
    write(file,"nbLMeasure",nbLMeasure)

    close(file) # Closing the file
end
########################################

@time tabLJrGrid!()
@time tabDNR!()

########################################
writedump!(namefile) # Dumping the computed table

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
using Missings

#=
Dictionary
=#

"""
Space complexity: O(n)
"""
mutable struct SimpleDict{K,V} <: AbstractDict{K,V}
    keys::Vector{Union{K, Missing}}
    values::Vector{Union{V, Missing}}
    count::Int

    function SimpleDict{K,V}(; capacity::Integer = 8) where {K,V}
        new{K,V}(missings(K, capacity), missings(V, capacity), 0)
    end
end

#=
Core methods
=#

# TODO collision policy

"""
Time complexity: O(1)
"""
function Base.setindex!(sd::SimpleDict{K,V}, value::V, key::K) where {K,V}
    index = hash_index(key, length(sd.keys))
    filled = !ismissing(sd.keys[index])
    sd.keys[index] = key
    sd.values[index] = value
    if !filled
        sd.count += 1
    end
    return sd
end

# TODO collision policy

"""
Time complexity: O(1)
"""
function Base.getindex(sd::SimpleDict{K,V}, key::K) where {K,V}
    haskey(sd, key) || throw(KeyError(key))
    index = hash_index(key, length(sd.keys))
    return sd.values[index]
end

hash_index(key, capacity::Integer) = ((hash(key) % Int) & (capacity - 1)) + 1

"""
Time complexity: O(1)
"""
function Base.delete!(sd::SimpleDict{K,V}, key::K) where {K,V}
    haskey(sd, key) || throw(KeyError(key))
    index = hash_index(key, length(sd.keys))
    sd.keys[index] = missing
    sd.values[index] = missing
    sd.count -= 1
    return sd
end

#=
Interface
=#

Base.length(sd::SimpleDict) = sd.count

function Base.haskey(sd::SimpleDict{K,V}, key::K) where {K,V}
    index = hash_index(key, length(sd.keys))
    k = sd.keys[index]
    if !ismissing(k)
        return k == key
    end
    return false
end

function Base.iterate(sd::SimpleDict{K,V}, i::Integer = 1) where {K,V}
    i = findnext(!ismissing, sd.keys, i) # find index of next filled, else return nothing
    if !isnothing(i)
        return (
            Pair{K,V}(sd.keys[i], sd.values[i]),
            i == typemax(Int) ? 0 : i + 1,
        )
    end
    return nothing
end


text:
## ------------------------------------------------------------------
# Meta

getdoc(s::RBSection) = get_parent(s)
setdoc!(s::RBSection, doc::RBDoc) = set_parent!(s, doc)

get_doi(s::RBSection) = getproperty!(s, :doi, "")
set_doi!(s::RBSection, doi::String) = setproperty!(s, :doi, _doi_to_url(doi))

get_title(s::RBSection) = getproperty!(s, :title, "")
set_title!(s::RBSection, title::String) = setproperty!(s, :title, strip(title))

## ------------------------------------------------------------------
# Data
items(s::RBSection) = getproperty!(() -> OrderedDict{String, RBObject}(), s, :items)
hasobj(s::RBSection, label::String) = haskey(items(s), label)

## ------------------------------------------------------------------
# Abstract Dict
Base.length(s::RBSection) = length(items(s))
Base.getindex(s::RBSection, label::String) = getindex(items(s), label)
Base.collect(s::RBSection) = collect(_values(items(s)))
Base.collect(T::Type, s::RBSection) = filter((i) -> i isa T, _values(items(s)))

## ------------------------------------------------------------------
# Array
Base.getindex(s::RBSection, idx) = (objs = items(s); getindex(_values(objs), idx))
Base.lastindex(s::RBSection) = lastindex(items(s))
Base.firstindex(s::RBSection) = firstindex(items(s))
Base.iterate(s::RBSection) = iterate(_values(items(s)))
Base.iterate(s::RBSection, state) = iterate(_values(items(s)), state)

## ------------------------------------------------------------------
# Show
function Base.show(io::IO, s::RBSection)
    nobjs = length(s)
    println(io, "RBSection(\"", get_label(s), "\") with ", nobjs, " object(s)")
    println(io, "path: '", localrelpath(s), "'")

    # meta
    for meta in [:title]
        str = getproperty(s, meta, "")
        if !isempty(str)
            println(io, meta, ": \"", _preview(io, str), "\"")
        end
    end

    # data
    if nobjs > 0
        print(io, "object(s):")
        _show_data_preview(io, s) do obj
            string("[", get_label(obj), "] ", typeof(obj))
        end
    end
    print(io, "\n", _preview(io, "-"^70))
    return nothing
end

## ------------------------------------------------------------------
# API
function eachobj(sec::RBSection)
    ch = Channel{RBObject}(1) do ch_
        for obj in sec
            put!(ch_, obj)
        end
    end

    return (obj for obj in ch)
end

# This should be the only method to add nre documents
# It updates all sort of stuff
function add_obj!(sec::RBSection, obj::RBObject)

    objlabel = get_label(obj)

    # Update stuff
    # obj labels
    add_objlabel!(obj)

    # add doc
    objs = items(sec)
    push!(objs, objlabel => obj)

    return sec
end


text:
mutable struct NilakanthaSeries
  p::Float64
  a::Float64
  s::Float64
  iterate::Function
  function NilakanthaSeries(_p, _a, _s)
    this = new(_p, _a, _s)
    this.iterate = function()
      s = this.s
      a = this.a
      this.p += s*(4.0/(a*(a*(a+3.0)+2.0)))
      this.a += 2.0
      this.s = -s
    end
    this
  end
end

approximator =  NilakanthaSeries(3.0, 2.0, 1.0)

for _ = 0:150000
  approximator.iterate()
end

println("pi = ",approximator.p)


text:
@enum ObjectiveType minimization=0 maximization=1 allDiff=2
abstract type ObjectiveFunction end

struct Range{T<:Real}
    max     ::T
    min     ::T
end

include("./sequence.jl")
include("./geometry/geometryMain.jl")
include("./finalModels/finalModelsMain.jl")

abstract type Node end
include("./nodes/nodesMain.jl")


text:
using BenchmarkTools

using DiffEqBase: remake
using Setfield: @lens

using Bifurcations
using Bifurcations: LimitCycleProblem
using Bifurcations.Continuations: pre_solve!, sweep!
using Bifurcations.Examples: DuffingVanDerPol

function make_van_der_pol_lc_solver()
    param = DuffingVanDerPol.DuffingVanDerPolParam(
        d = 0.1,
    )
    ode = remake(
        DuffingVanDerPol.ode,
        p = param,
    )

    num_mesh = 20
    degree = 3

    # Start continuation of the limit cycle at the small damping limit.
    xs0 = let
        n = num_mesh * degree
        dt = 2π / n
        ts = ((1:n) - 1) * dt
        hcat(
            2 .* cos.(ts),
            -2 .* sin.(ts),
        )'
    end

    prob = LimitCycleProblem(
        de_prob = ode;
        param_axis = DuffingVanDerPol.param_axis,
        num_mesh = num_mesh,
        degree = degree,
        xs0 = xs0,
        l0 = 2π,
        t_domain = (0.01, 1.5),
        t0 = get(ode.p, DuffingVanDerPol.param_axis),
    )

    @assert size(prob.xs0) == (2, num_mesh * degree)

    return init(
        prob;
        start_from_nearest_root = true,
        max_samples = 3,
    )
end

suite = BenchmarkGroup()
suite["sweep!(lc_solver)"] = @benchmarkable(
    sweep!(lc_solver),
    setup = (lc_solver = pre_solve!(make_van_der_pol_lc_solver())))

if !isdefined(:SUITE)
    SUITE = BenchmarkGroup()
end
SUITE["van_der_pol"] = suite