def METHOD_NAME(g_s, g_x, g_e, f_s, f_x, f_S, f_X):
    n_s = g_s.shape[0]  # number of controls
    n_x = g_x.shape[1]  # number of states
    n_e = g_e.shape[1]
    n_v = n_s + n_x
    A = row_stack(
        [column_stack([eye(n_s), zeros((n_s, n_x))]), column_stack([-f_S, -f_X])]
    )
    B = row_stack([column_stack([g_s, g_x]), column_stack([f_s, f_x])])
    [S, T, Q, Z, eigval] = qzordered(A, B, 1.0 - 1e-8)
    Z = Z.real
    diag_S = np.diag(S)
    diag_T = np.diag(T)
    tol_geneigvals = 1e-10
    try:
        ok = sum((abs(diag_S) < tol_geneigvals) * (abs(diag_T) < tol_geneigvals)) == 0
        assert ok
    except Exception as e:
        raise GeneralizedEigenvaluesError(diag_S=diag_S, diag_T=diag_T)
    eigval_s = sorted(eigval, reverse=False)
    if max(eigval[:n_s]) >= 1 and min(eigval[n_s:]) < 1:
        # BK conditions are met
        pass
    else:
        ev_a = eigval_s[n_s - 1]
        ev_b = eigval_s[n_s]
        cutoff = (ev_a - ev_b) / 2
        if not ev_a > ev_b:
            raise GeneralizedEigenvaluesSelectionError(
                A=A,
                B=B,
                eigval=eigval,
                cutoff=cutoff,
                diag_S=diag_S,
                diag_T=diag_T,
                n_states=n_s,
            )
        import warnings
        if cutoff > 1:
            warnings.warn("Solution is not convergent.")
        else:
            warnings.warn(
                "There are multiple convergent solutions. The one with the smaller eigenvalues was selected."
            )
        [S, T, Q, Z, eigval] = qzordered(A, B, cutoff)
    Z11 = Z[:n_s, :n_s]
    # Z12 = Z[:n_s, n_s:]
    Z21 = Z[n_s:, :n_s]
    # Z22 = Z[n_s:, n_s:]
    # S11 = S[:n_s, :n_s]
    # T11 = T[:n_s, :n_s]
    # first order solution
    # P = (solve(S11.T, Z11.T).T @ solve(Z11.T, T11.T).T)
    C = solve(Z11.T, Z21.T).T
    A = g_s + g_x @ C
    B = g_e
    return C, eigval_s