Bug in slicot TB05AD

[rabraker]

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.

[murrayrm]

@rabraker Do we need to fix before doing the next release of slycot or is this a minor issue that can wait for later? I'm assuming the latter.

[rabraker]

I'm not terribly familiar with how the releases work...

The TB05AD function should return correct results if one uses the 'NG' and 'NH' flags (which are the flags used in the the recent PR 177 on python-control).

The bug is "documented" with a known failure in the unit tests. But perhaps a comment about this bug using the balancing option 'AG' could be added to the doc-string?

[repagh]

Copy & paste issues?

I grep'd DGEBAL in the Fortran sources. The results

• TB01TD found the same presumed erroneous structure.
• TB05AD reported above. I think, from reading the docs, that you are correct. These are permutations that need to be carried over to the B and C matrices, and thus in the same order as those on the A.
• MB04MD is an alternative version for the DGEBAL routine. It is used internally in SLICOT, in the routine MB05OD.
• MB05OY is to reverse the effects of DGEBAL, and it should thus work in the reverse order, which it does. It also is only used internally by MB05OD
• AB13DD has the correct order, down from N and then up from 1.
• AB13FD, MB03SD and MB04DY only use scaling or don't do back transformation, thus no issues here.

The curious thing is that variable names are different on all the pieces where permutation is applied or reversed, whether (assumed) correct or wrong. This may be due to the habit of re-using integer variables with Fortran.

I guess the TB01TD and TB05AD need to be corrected. Now trying to find some test cases.

[repagh]

PR #122 should have solved this

@rabraker, can you confirm?