diff --git a/slycot/src/TB01TD.f b/slycot/src/TB01TD.f index 7c52957a..8f2fc650 100644 --- a/slycot/src/TB01TD.f +++ b/slycot/src/TB01TD.f @@ -245,7 +245,7 @@ SUBROUTINE TB01TD( N, M, P, A, LDA, B, LDB, C, LDC, D, LDD, LOW, C to transform B and C. C DO 10 K = 1, N - KOLD = K + KOLD = N + 1 - K ! RvP, rabraker, slycot #11 IF ( ( LOW.GT.KOLD ) .OR. ( KOLD.GT.IGH ) ) THEN IF ( KOLD.LT.LOW ) KOLD = LOW - KOLD KNEW = INT( SCSTAT(KOLD) ) diff --git a/slycot/src/TB05AD.f b/slycot/src/TB05AD.f index c7b93e91..55490b84 100644 --- a/slycot/src/TB05AD.f +++ b/slycot/src/TB05AD.f @@ -346,10 +346,10 @@ SUBROUTINE TB05AD( BALEIG, INITA, N, M, P, FREQ, A, LDA, B, LDB, C defined in the subroutine DGEBAL. C DO 10 J = 1, N - JJ = J + JJ = N + 1 - J ! RvP, rabraker, slycot #11 IF ( JJ.LT.LOW .OR. JJ.GT.IGH ) THEN IF ( JJ.LT.LOW ) JJ = LOW - JJ - JP = DWORK(JJ) + JP = INT( DWORK(JJ) ) IF ( JP.NE.JJ ) THEN C C Permute rows of B. diff --git a/slycot/tests/test_tb05ad.py b/slycot/tests/test_tb05ad.py index e64f7360..b5330cca 100644 --- a/slycot/tests/test_tb05ad.py +++ b/slycot/tests/test_tb05ad.py @@ -11,6 +11,7 @@ from numpy.testing import assert_almost_equal + # set the random seed so we can get consistent results. np.random.seed(40) CASES = {} @@ -47,9 +48,10 @@ def test_tb05ad_ng(self): for key in CASES: sys = CASES[key] self.check_tb05ad_AG_NG(sys, 10*1j, 'NG') + self.check_tb05ad_AG_NG(sys, 10*1j, 'AG') - @unittest.expectedFailure - def test_tb05ad_ag_failure(self): + #@unittest.expectedFailure + def test_tb05ad_ag_no_longer_failure(self): """ Test tb05ad and job 'AG' (i.e., balancing enabled) fails on certain A matrices. """ @@ -187,6 +189,29 @@ def test_tb05ad_resonance(self): transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH') assert cm.exception.info == 2 + def test_tb05ad_balance(self): + """Specifically check for the balancing output. + + Based on https://rdrr.io/rforge/expm/man/balance.html + """ + A = np.array(((-1, -1, 0, 0), + ( 0, 0, 10, 10), + ( 0, 0, 10, 0), + ( 0, 10, 0, 0)), dtype=float) + B = np.eye(4) + C = np.eye(4) + Ar = np.array(((-1, -1, 0, 0), + ( 0, 0, 10, 10), + ( 0, 10, 0, 0), + ( 0, 0, 0, 10))) + jomega = 1.0 + At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad( + 4, 4, 4, jomega, A, B, C, job='AG') + assert_almost_equal(At, Ar) + + At, Bt, Ct, rcond, g_jwb, ev, hinvb, info = transform.tb05ad( + 4, 4, 4, jomega, A, B, C, job='AG') + if __name__ == "__main__": unittest.main()