Defaults to None (no sorting).
    overwrite_a : bool, optional
        Whether to overwrite data in a (may improve performance)
    overwrite_b : bool, optional
        Whether to overwrite data in b (may improve performance)
    check_finite : bool, optional
        If true checks the elements of `A` and `B` are finite numbers. If
        false does no checking and passes matrix through to
        underlying algorithm.

    Returns
    -------
    AA : (N, N) ndarray
        Generalized Schur form of A.
    BB : (N, N) ndarray
        Generalized Schur form of B.
    Q : (N, N) ndarray
        The left Schur vectors.
    Z : (N, N) ndarray
        The right Schur vectors.

    See Also
    --------
    ordqz

    Notes
    -----
    Q is transposed versus the equivalent function in Matlab.

    .. versionadded:: 0.11.0

    Examples
    --------
    >>> import numpy as np
    >>> from scipy.linalg import qz

    >>> A = np.array([[1, 2, -1], [5, 5, 5], [2, 4, -8]])
    >>> B = np.array([[1, 1, -3], [3, 1, -1], [5, 6, -2]])

    Compute the decomposition.  The QZ decomposition is not unique, so
    depending on the underlying library that is used, there may be
    differences in the signs of coefficients in the following output.

    >>> AA, BB, Q, Z = qz(A, B)
    >>> AA
    array([[-1.36949157, -4.05459025,  7.44389431],
           [ 0.        ,  7.65653432,  5.13476017],
           [ 0.        , -0.65978437,  2.4186015 ]])  # may vary
    >>> BB
    array([[ 1.71890633, -1.64723705, -0.72696385],
           [ 0.        ,  8.6965692 , -0.        ],
           [ 0.        ,  0.        ,  2.27446233]])  # may vary
    >>> Q
    array([[-0.37048362,  0.1903278 ,  0.90912992],
           [-0.90073232,  0.16534124, -0.40167593],
           [ 0.22676676,  0.96769706, -0.11017818]])  # may vary
    >>> Z
    array([[-0.67660785,  0.63528924, -0.37230283],
           [ 0.70243299,  0.70853819, -0.06753907],
           [ 0.22088393, -0.30721526, -0.92565062]])  # may vary

    Verify the QZ decomposition.  With real output, we only need the
    transpose of ``Z`` in the following expressions.

    >>> Q @ AA @ Z.T  # Should be A
    array([[ 1.,  2., -1.],
           [ 5.,  5.,  5.],
           [ 2.,  4., -8.]])
    >>> Q @ BB @ Z.T  # Should be B
    array([[ 1.,  1., -3.],
           [ 3.,  1., -1.],
           [ 5.,  6., -2.]])

    Repeat the decomposition, but with ``output='complex'``.

    >>> AA, BB, Q, Z = qz(A, B, output='complex')

    For conciseness in the output, we use ``np.set_printoptions()`` to set
    the output precision of NumPy arrays to 3 and display tiny values as 0.

    >>> np.set_printoptions(precision=3, suppress=True)
    >>> AA
    array([[-1.369+0.j   ,  2.248+4.237j,  4.861-5.022j],
           [ 0.   +0.j   ,  7.037+2.922j,  0.794+4.932j],
           [ 0.   +0.j   ,  0.   +0.j   ,  2.655-1.103j]])  # may vary
    >>> BB
    array([[ 1.719+0.j   , -1.115+1.j   , -0.763-0.646j],
           [ 0.   +0.j   ,  7.24 +0.j   , -3.144+3.322j],
           [ 0.   +0.j   ,  0.   +0.j   ,  2.732+0.j   ]])  # may vary
    >>> Q
    array([[ 0.326+0.175j, -0.273-0.029j, -0.886-0.052j],
           [ 0.794+0.426j, -0.093+0.134j,  0.402-0.02j ],
           [-0.2  -0.107j, -0.816+0.482j,  0.151-0.167j]])  # may vary
    >>> Z
    array([[ 0.596+0.32j , -0.31 +0.414j,  0.393-0.347j],
           [-0.619-0.332j, -0.479+0.314j,  0.154-0.393j],
           [-0.195-0.104j,  0.576+0.27j ,  0.715+0.187j]])  # may vary

    With complex arrays, we must use ``Z.conj().T`` in the following
    expressions to verify the decomposition.

    >>> Q @ AA @ Z.conj().T  # Should be A
    array([[ 1.-0.j,  2.-0.j, -1.-0.j],
           [ 5.+0.j,  5.+0.j,  5.-0.j],
           [ 2.+0.j,  4.+0.j, -8.+0.j]])
    >>> Q @ BB @ Z.conj().T  # Should be B
    array([[ 1.+0.j,  1.+0.j, -3.+0.j],
           [ 3.-0.j,  1.-0.j, -1.+0.j],
           [ 5.+0.j,  6.+0.j, -2.+0.j]])

    )