Description
I recently submitted a pull request for a wrapper for TB05AD. The TB05AD routine is used to efficiently compute the frequency response of a state space system at many frequencies. The first time you call TB05AD.f, the routine will transform the system matrices such that the A matrix is in upper Hessenberg form. These transformed system matrices are returned and should then be used in subsequent calls to TB05AD.f.
One of the job options for TB05AD.f is BALEIG='N' or 'C' or 'B' or 'A'. My wrapper currently only uses options 'N' and 'A'.
When using option 'A', TB05AD.f will balance the A matrix before transforming A to upper Hessenberg form. For certain A matrices, this results in an incorrect frequency response calculation. For example:
import numpy as np
from scipy import linalg
from slycot import transform
A = np.array([[-0.5, 0., 0., 0. ],
[ 0., -1., 0. , 0. ],
[ 1., 0., -0.5, 0. ],
[ 0., 1., 0., -1. ]])
B = np.array([[ 1.],
[ 0.],
[ 0.],
[ 0.]])
C = np.array([[ 0., 1., 1., 0.]])
freq = 10*1j
n, m = B.shape
p = C.shape[0]
result = transform.tb05ad(n, m, p, freq, A, B, C, job='AG')
g= result[3]
g_true = C.dot(np.linalg.solve(np.eye(n) * freq - A, B))
print g - g_true
I believe this is a bug in TB05AD.f circa line 354. When doing the balancing, TB05AD calls the lapack routine DGEBAL. From what I can tell, TB05AD misinterprets the permutation information provided by DGEBAL. This is similar to the scipy matrix_balance bug which I reported here. Before that bug was fixed, linalg.matrix_balance + linalg.hessenberg and TB05AD would yield the same transformed set of (Abar, Bbar, Cbar), but both would yield the wrong frequency response.
The documentation for DGEBAL says that rows/cols are permuted going from N to IGH+1, then from 1 to LO-1. I think TB05AD does the backwards. In some preliminary tests, replacing line 349 with
JJ = N + 1 - J
seems to fix the issue, but this needs further testing.