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

repagh · May 11, 2020 · fca28a3 · fca28a3
1 parent 429989b
commit fca28a3
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 5 deletions.
diff --git 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
@@ -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
@@ -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()