scipy/scipy — pull request #3107 by jseabold (2013-12-02 17:33:48+00:00)
Title: WIP: Add ordered QZ decomposition

This currently segfaults. I don't know why. I spent some time checking the wrapper definitions and stepping through with gdb. I didn't notice anything that was obvious to me. If someone else wants to give it a go, this was the working script I was using. If someone can figure out the segfault, I'll take it the rest of the way, but I just don't have time to mess with this.

```
A = [[-21.10-22.50j, 53.5-50.5j, -34.5+127.5j, 7.5+ 0.5j],
[ -0.46- 7.78j, -3.5-37.5j, -15.5+ 58.5j, -10.5- 1.5j],
[ 4.30- 5.50j, 39.7-17.1j, -68.5+ 12.5j, -7.5- 3.5j],
[ 5.50+ 4.40j, 14.4+43.3j, -32.5- 46.0j, -19.0-32.5j]]
B = [[1.0-5.0j, 1.6+1.2j, -3+0j, 0.0-1.0j],
[0.8-0.6j, 3.0-5.0j, -4+3j, -2.4-3.2j],
[1.0+0.0j, 2.4+1.8j, -4-5j, 0.0-3.0j],
[0.0+1.0j, -1.8+2.4j, 0-4j, 4.0-5.0j]]

A = [[3.9, 12.5, -34.5, -0.5],
[4.3, 21.5, -47.5, 7.5],
[4.3, 21.5, -43.5, 3.5],
[4.4, 26.0, -46.0, 6.0]]

B = [[ 1, 2, -3, 1],
[ 1, 3, -5, 4],
[ 1, 3, -4, 3],
[ 1, 3, -4, 4]]

from scipy.linalg._decomp_qz import ordqz, _qz, _select_function, get_lapack_funcs
output = 'real'
lwork = None
overwrite_a = False
overwrite_b = False
check_finite = True
sort = 'lhp'
result, typ = _qz(A, B, output=output, lwork=lwork, sort=None,
overwrite_a=overwrite_a, overwrite_b=overwrite_b,
check_finite=check_finite)
AA, BB, Q, Z = result[0], result[1], result[-4], result[-3]
alpha, beta = result[3], result[4]
sfunction = _select_function(sort, typ)
select = sfunction(alpha, beta)
tgsen, = get_lapack_funcs(('tgsen',), (AA,BB))

n = len(A)
m = 0.
lwork = 2*m*(n-m)
liwork = n+2
#result = tgsen(select, AA, BB, Q, Z, lwork=-1, liwork=AA.shape[1]*2)
result = tgsen(select, AA, BB, Q, Z)
```