- corrected TD04AD.f to handle the case of 'C' and static gain tf

- corrected transform.pyf, to not allocate zero-dim matrices - added tests for td04ad, one static, one checked against Octave with a transfer function set.
repagh · Jun 4, 2018 · Jun 6, 2018 · Jun 14, 2018 · Jun 15, 2018 · Jun 15, 2018
commit f897e7526a0c28c8486dd878f20c78bf50d06bec
diff --git a/slycot/src/TD04AD.f b/slycot/src/TD04AD.f
@@ -245,7 +245,8 @@ SUBROUTINE TD04AD( ROWCOL, M, P, INDEX, DCOEFF, LDDCOE, UCOEFF,
      $                  UCOEFF(LDUCO1,LDUCO2,*)
 C     .. Local Scalars ..
       LOGICAL           LROCOC, LROCOR
-      INTEGER           I, J, JSTOP, K, KDCOEF, MPLIM, MWORK, N, PWORK
+      INTEGER           I, J, JSTOP, K, KDCOEF, MPLIM, MWORK, N, PWORK,
+     $                  KU
 C     .. External Functions ..
       LOGICAL           LSAME
       EXTERNAL          LSAME
@@ -402,17 +403,27 @@ SUBROUTINE TD04AD( ROWCOL, M, P, INDEX, DCOEFF, LDDCOE, UCOEFF,
 C        to get an upper block Hessenberg state dynamics matrix.
 C
          K = IWORK(1)+IWORK(2)-1
-         CALL TB01XD( 'D', NR, MWORK, PWORK, K, NR-1, A, LDA, B, LDB,
+C
+C        RvP 180615 Try to protect against re-working an empty [] A
+C        matrix, failed with K < 0
+C                  
+         IF ( NR.EQ.0 ) THEN
+            K = 0
+            KU = 0
+         ELSE
+            KU = NR - 1
+         END IF
+         CALL TB01XD( 'D', NR, MWORK, PWORK, K, KU, A, LDA, B, LDB,
      $                C, LDC, D, LDD, INFO )
          IF ( MPLIM.NE.1 ) THEN
 C
 C           Also, retranspose U(s) if this is non-scalar.
 C
             DO 70 K = 1, KDCOEF
-C
+C 
                DO 60 J = 1, JSTOP
                   CALL DSWAP( MPLIM-J, UCOEFF(J+1,J,K), 1,
-     $                        UCOEFF(J,J+1,K), LDUCO1 )
+     $                           UCOEFF(J,J+1,K), LDUCO1 )
    60          CONTINUE
 C
    70       CONTINUE

diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf
@@ -324,7 +324,7 @@ subroutine td04ad_r(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
     integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
     integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
     integer intent(in,out) :: nr !=sum(index_bn)
-    double precision intent(out),dimension(max(1,nr),nr),depend(nr) :: a
+    double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
     integer intent(hide),depend(a) :: lda = shape(a,0) 
     double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
     integer intent(hide),depend(b) :: ldb = shape(b,0)
@@ -354,7 +354,7 @@ subroutine td04ad_c(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
     integer intent(hide),depend(a) :: lda = shape(a,0) 
     double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
     integer intent(hide),depend(b) :: ldb = shape(b,0)
-    double precision intent(out),dimension(max(1,max(m,p)),nr),depend(nr,m,p) :: c
+    double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
     integer intent(hide),depend(c) :: ldc = shape(c,0)
     double precision intent(out),dimension(max(1,max(m,p)),max(m,p)),depend(p,m) :: d
     integer intent(hide),depend(d) :: ldd = shape(d,0)

diff --git a/slycot/tests/test_td04ad.py b/slycot/tests/test_td04ad.py
@@ -11,7 +11,7 @@
 
 from numpy.testing import assert_raises, assert_almost_equal
 
-class test_td04ad(unittest.TestCase):
+class TestTf2SS(unittest.TestCase):
 
     def test_td04ad_case1(self):
         # for octave:
@@ -32,8 +32,10 @@ def test_td04ad_case1(self):
             [ [0.0,  0.0, 1.0 ], [ 1.0, 0.0, 0.0 ] ],
             [ [3.0, -1.0, 1.0 ], [ 0.0, 1.0, 0.0 ] ],
             [ [0.0, 0.0, 1.0],   [ 0.0, 2.0, 0.0 ] ] ])
-        numc = np.zeros((3, 3, 3),dtype=float)
-        numc[:,:2,:] = num
+        p, m, d = num.shape
+        numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
+        numc[:p,:m,:] = num
+
         denc = np.array(
             [ [1.0,  0.4, 3.0],  [ 1.0, 1.0, 0.0 ] ])
         indc = np.array(
@@ -44,13 +46,59 @@ def test_td04ad_case1(self):
             [ 2, 1, 0 ], dtype=int)
 
         n, A, B, C, D = transform.td04ad('C', 2, 3, indc, denc, numc)
-        print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
-
+        #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
+        Ac = [ [-1, 0, 0], [ 0, -0.4, -0.3], [ 0, 10, 0]]
+        Bc = [ [0, -1] ,[ 1 , 0], [ 0, 0]]
+        Cc = [ [1, 0, 0.1], [-1, -2.2, -0.8], [ -2, 0, 0.1] ]
+        Dc = [ [0, 1], [ 3, 0], [ 0, 0]]
+        np.testing.assert_array_almost_equal(A, Ac)
+        np.testing.assert_array_almost_equal(B, Bc)
+        np.testing.assert_array_almost_equal(C, Cc)
+        np.testing.assert_array_almost_equal(D, Dc)
+
         resr = transform.td04ad('R', 2, 3, indr, denr, num)
-        print(resr)
+        #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
+
+    def test_staticgain(self):
+
+        # 2 inputs, 3 outputs? columns share a denominator
+        num = np.array([ [ [1.0], [2.0] ],
+                         [ [0.2], [4.3] ],
+                         [ [1.2], [3.2] ] ])
+        p, m, d = num.shape
+        numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
+        numc[:p,:m,:] = num
+
+        # denc, columns share a common denominator
+        denc = np.array([ [ 1.0], [0.5] ])
+        Dc = (num / denc).reshape((3,2))
+        idxc = np.zeros((2,), dtype=int)
+
+        # denr, rows share a common denominator
+        denr = np.array([ [1.0], [0.5], [3.0] ])
+        idxr = np.zeros((3,), dtype=int)
+        Dr = (num / denr[:, np.newaxis]).reshape((3,2))
+
+        # fails with:
+        # On entry to TB01XD parameter number  5 had an illegal value
+
+        n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc)
+        #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
+        self.assertEqual(A.shape, (0,0))
+        self.assertEqual(B.shape, (0,2))
+        self.assertEqual(C.shape, (3,0))
+        np.testing.assert_array_almost_equal(D, Dc)
+
+        n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num)
+        #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
+        self.assertEqual(A.shape, (0,0))
+        self.assertEqual(B.shape, (0,2))
+        self.assertEqual(C.shape, (3,0))
+        np.testing.assert_array_almost_equal(D, Dr)
+
 
 def suite():
-   return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
+   return unittest.TestLoader().loadTestsFromTestCase(TestTF2SS)
 
 
 if __name__ == "__main__":