Skip to content

Bug in slicot TB05AD #11

Closed
Closed
Bug in slicot TB05AD#11
Labels
bug

Description

@rabraker
rabraker
opened on May 19, 2017

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.

Activity

  1. murrayrm commented on Jan 13, 2018

    @murrayrm
    murrayrm
    on Jan 13, 2018
    Member

    @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.

  2. murrayrm
    added
    bug
    on Jan 13, 2018

  3. rabraker commented on Jan 14, 2018

    @rabraker
    rabraker
    on Jan 14, 2018
    ContributorAuthor

    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?

  4. repagh commented on May 11, 2020

    @repagh
    repagh
    on May 11, 2020
    Member

    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.

  5. repagh
    added a commit that references this issue on May 11, 2020

    fix for DGEBAL transformation application, issue python-control#11

    fca28a3

  6. repagh commented on May 21, 2020

    @repagh
    repagh
    on May 21, 2020
    Member

    PR #122 should have solved this

    @rabraker, can you confirm?

  7. ilayn commented on Oct 27, 2025

    @ilayn
    ilayn
    on Oct 27, 2025

    This was probably due to the bug that caused some issues present all the way to LAPACK version 3.4.2 and fixed in LAPACK 3.5. It might not be an issue anymore.

  8. bnavigator commented on Jun 16, 2026

    @bnavigator
    bnavigator
    on Jun 16, 2026
    Collaborator

    Closing as stale, feel free to reopen, if you encounter it with a more recent version.

  9. bnavigator
    closed this as completedon Jun 16, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bug

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions