tgsen uses an output argument for computing a default argument

[Sbte]

tgsen uses lwork=max(2*m*(n-m),1) as default argument (for the complex case) where integer intent(out) :: m. I assume this means that lwork will always be 1, which is usually insufficient. m is the amount of eigenpairs that will be selected (which may be larger than the amount of nonzeros in the select array in the non-complex case).

Note that the same issue exists for the non-complex case, but this is less visible, since the default is set to lwork=max(MAX(4*n+16,2*m*(n-m)),1), so this includes 4*n+16.

Reproducing code example:

@pytest.mark.parametrize('dtype', DTYPES)
def test_gges_tgsen(dtype):
    seed(1234)
    atol = np.finfo(dtype).eps*100

    n = 10
    a = generate_random_dtype_array([n, n], dtype=dtype)
    b = generate_random_dtype_array([n, n], dtype=dtype)

    gges, tgsen = get_lapack_funcs(('gges', 'tgsen'), dtype=dtype)

    result = gges(lambda x: None, a, b, overwrite_a=False, overwrite_b=False)
    assert_equal(result[-1], 0)

    s = result[0]
    t = result[1]
    q = result[-4]
    z = result[-3]

    if dtype in COMPLEX_DTYPES:
        assert_allclose(s, np.triu(s), rtol=0, atol=atol)
        assert_allclose(t, np.triu(t), rtol=0, atol=atol)

    assert_allclose(q @ s @ z.conj().T, a, rtol=0, atol=atol)
    assert_allclose(q @ t @ z.conj().T, b, rtol=0, atol=atol)

    result = tgsen(select, s, t, q, z)
    assert_equal(result[-1], 0)

    s = result[0]
    t = result[1]
    q = result[-9]
    z = result[-8]

    if dtype in COMPLEX_DTYPES:
        assert_allclose(s, np.triu(s), rtol=0, atol=atol)
        assert_allclose(t, np.triu(t), rtol=0, atol=atol)

    assert_allclose(q @ s @ z.conj().T, a, rtol=0, atol=atol)
    assert_allclose(q @ t @ z.conj().T, b, rtol=0, atol=atol)

Error message:

 ** On entry to CTGSEN parameter number 21 had an illegal value
 ** On entry to ZTGSEN parameter number 21 had an illegal value

Scipy/Numpy/Python version information:

Scipy commit 55cf5c301e58f24092639c01e48c2bd139807906 Numpy 1.19.4 Python 3.8.5

[Sbte]

Also, it needs to be +1 the actual default value. See also

https://github.com/scipy/scipy/blob/6e15c52a40370d629dfa4951eb58e7f2b267ed61/scipy/linalg/_decomp_qz.py#L372-L374

[ilayn]

Well indeed we clearly made a mistake in that wrapper.

[ilayn]

Correction: I am surprised that it is actually working. Most settings should have lead to segfaults. I'll fix this

[ilayn]

Well apparently, I've opened a can of worms since when I touched these wrappers, linalg.qz and linalg.ordqz stopped working. So it will get a bit complicated as a fix.

However, there is basically no way to get a straightforward formula to get the minimal lwork and liwork values hence there will be an additional ?tgsen_lwork functions wrapped to get the right values. Then the regular ?tgsen call can be requested.