Merge pull request #79 from bnavigator/fix-unused-vars

  - fixes and tests for mb05md, mb05nd, mc01td
  - removed unnecessary mathematical.pyf

slycot/CMakeLists.txt

@@ -104,7 +104,7 @@ set(FSOURCES
 
 set(F2PYSOURCE src/_wrapper.pyf)
 set(F2PYSOURCE_DEPS
-  src/analysis.pyf src/math.pyf src/mathematical.pyf
+  src/analysis.pyf src/math.pyf
   src/transform.pyf src/synthesis.pyf)
 
 configure_file(version.py.in version.py @ONLY)

slycot/math.py

@@ -18,6 +18,7 @@
 #       MA 02110-1301, USA.
 
 from . import _wrapper
+import warnings
 
 def mc01td(dico,dp,p):
     """ dp,stable,nz = mc01td(dico,dp,p)
@@ -61,16 +62,16 @@ def mc01td(dico,dp,p):
         e.info = out[-1]
         raise e
     if out[-1] == 1:
-        warings.warn('entry P(x) is the zero polynomial.')
+        warnings.warn('entry P(x) is the zero polynomial.')
     if out[-1] == 2:
-        warings.warn('P(x) may have zeros very close to stability boundary.')
+        warnings.warn('P(x) may have zeros very close to stability boundary.')
     if out[-2] > 0:
-        warnings.warn('The degree of P(x) has been reduced to %i' %(dp-k))
+        warnings.warn('The degree of P(x) has been reduced to %i' %(dp-out[-2]))
     return out[:-2]
 
 
 def mb05md(a, delta, balanc='N'):
-    """Ar, Vr, Yr, VALRr, VALDr = mb05md(a, delta, balanc='N')
+    """Ar, Vr, Yr, VAL = mb05md(a, delta, balanc='N')
 
     Matrix exponential for a real non-defective matrix
 
@@ -88,7 +89,7 @@ def mb05md(a, delta, balanc='N'):
             Square matrix
         delta : input 'd'
             The scalar value delta of the problem.
-     
+
     Optional arguments:
         balanc : input char*1
             Indicates how the input matrix should be diagonally scaled
@@ -114,7 +115,7 @@ def mb05md(a, delta, balanc='N'):
             (k+1)-th columns of the eigenvector matrix, respectively,
             then the eigenvector corresponding to the complex
             eigenvalue with positive (negative) imaginary value is
-            given by 
+            given by
               p + q*j (p - q*j), where j^2  = -1.
         Yr : output rank-2 array('d') with bounds (n,n)
             contains an intermediate result for computing the matrix
@@ -126,29 +127,39 @@ def mb05md(a, delta, balanc='N'):
             the (right) eigenvector matrix of A, where Lambda is the
             diagonal matrix of eigenvalues.
 
-        VALr : output rank-1 array('c') with bounds (n)
+        VAL : output rank-1 array('c') with bounds (n)
             Contains the eigenvalues of the matrix A. The eigenvalues
             are unordered except that complex conjugate pairs of values
             appear consecutively with the eigenvalue having positive
             imaginary part first.
     """
     hidden = ' (hidden by the wrapper)'
-    arg_list = [ 'balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden,
-                 'y','ldy'+hidden,'valr','vali',
-                 'iwork'+hidden,'dwork'+hidden,'ldwork'+hidden,'info'+hidden]
-    out = _wrapper.mb05md(balanc=balanc,n=min(a.shape),delta=delta,a=a)
-    if out[-1] == 0:
-        return out[:-1]
-    elif out[-1] < 0:
-        error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
-    elif out[-1] > 0 and out[-1] <= n:
-        error_text = "Incomplete eigenvalue calculation, missing %i eigenvalues" % out[-1]
-    elif out[-1] == n+1:
+    arg_list = ['balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden,
+                'y', 'ldy'+hidden, 'valr', 'vali',
+                'iwork'+hidden, 'dwork'+hidden, 'ldwork'+hidden,
+                'info'+hidden]
+    n = min(a.shape)
+    (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md(balanc=balanc,
+                                                     n=n,
+                                                     delta=delta,
+                                                     a=a)
+    if INFO == 0:
+        if not all(VALi == 0):
+            VAL = VALr + 1J*VALi
+        else:
+            VAL = VALr
+        return (Ar, Vr, Yr, VAL)
+    elif INFO < 0:
+        error_text = "The following argument had an illegal value: " \
+                     + arg_list[-INFO-1]
+    elif INFO > 0 and INFO <= n:
+        error_text = "Incomplete eigenvalue calculation, missing %i eigenvalues" % INFO
+    elif INFO == n+1:
         error_text = "Eigenvector matrix singular"
-    elif out[-1] == n+2:
+    elif INFO == n+2:
         error_text = "A matrix defective"
     e = ValueError(error_text)
-    e.info = out[-1]
+    e.info = INFO
     raise e
 
 """
@@ -182,21 +193,21 @@ def mb05nd(a, delta, tol=1e-7):
         H : Int[F(s) ds] from s = 0 to s = delta,
     """
     hidden = ' (hidden by the wrapper)'
-    arg_list = [ 'n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden,
-                 'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden,
-                 'dwork'+hidden, 'ldwork'+hidden]
-    out = _wrapper.mb05nd(n=min(a.shape), delta=delta, a=a, tol=tol)
+    arg_list = ['n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden,
+                'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden,
+                'dwork'+hidden, 'ldwork'+hidden]
+    n = min(a.shape)
+    out = _wrapper.mb05nd(n=n, delta=delta, a=a, tol=tol)
     if out[-1] == 0:
         return out[:-1]
     elif out[-1] < 0:
         error_text = "The following argument had an illegal value: " \
-                     +arg_list[-out[-1]-1]
+                     + arg_list[-out[-1]-1]
     elif out[-1] == n+1:
         error_text = "Delta too large"
     e = ValueError(error_text)
     e.info = out[-1]
     raise e
-
 
 
 # to be replaced by python wrappers

slycot/src/math.pyf

@@ -3,7 +3,7 @@
 
 subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f
 	character :: dico
-	integer intent(in,out),check(dp>0) :: dp
+	integer intent(in,out),check(dp>=0) :: dp
 	double precision intent(in),check(shape(p,0)==dp+1),dimension(dp+1),depend(dp) :: p
 	logical intent(out) :: stable
 	integer intent(out) :: nz

slycot/synthesis.py

@@ -689,20 +689,20 @@ def sb02od(n,m,A,B,Q,R,dico,p=None,L=None,fact='N',uplo='U',sort='S',tol=0.0,ldw
         out = _wrapper.sb02od_n(dico,n,m,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork)
     if fact == 'C':
         if p is None:
-            p = shape(Q)[0]
+            p = _np.shape(Q)[0]
         out = _wrapper.sb02od_c(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork)
     if fact == 'D':
         if p is None:
-            p = shape(R)[0]
+            p = _np.shape(R)[0]
         out = _wrapper.sb02od_d(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork)
     if fact == 'B':
         if p is None:
-            p = shape(Q)[0]
+            p = _np.shape(Q)[0]
         out = _wrapper.sb02od_b(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork)
     if out[-1] < 0:
         error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
         e = ValueError(error_text)
-        e.info = info
+        e.info = out[-1]
         raise e
     if out[-1] == 1:
         e = ValueError('the computed extended matrix pencil is singular, possibly due to rounding errors')

slycot/tests/CMakeLists.txt

@@ -3,6 +3,8 @@ set(PYSOURCE
   __init__.py
   test.py
   test_ag08bd.py
+  test_mb.py
+  test_mc.py
   test_sb10jd.py
   test_sg02ad.py
   test_sg03ad.py
@@ -11,4 +13,8 @@ set(PYSOURCE
   test_tg01ad.py
   test_tg01fd.py )
 
-install(FILES ${PYSOURCE} DESTINATION slycot/tests)
+install(FILES ${PYSOURCE}
+        PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE
+                    GROUP_READ GROUP_EXECUTE
+                    WORLD_READ WORLD_EXECUTE
+        DESTINATION slycot/tests)

slycot/tests/test.py

@@ -1,6 +1,5 @@
 import unittest
 from slycot import synthesis
-from slycot import math
 from slycot import transform
 
 class Test(unittest.TestCase):
@@ -11,11 +10,6 @@ def setUp(self):
     def test_1(self):
         synthesis.sb02mt(1,1,1,1)
 
-    def test_2(self):
-        from numpy import array
-        a = array([[-2, 0.5], [-1.6, -5]])
-        Ar, Vr, Yr, VALRr, VALDr = math.mb05md(a, 0.1)
-
     def test_sb02ad(self):
         "Test sb10ad, Hinf synthesis"
         import numpy as np
@@ -64,9 +58,5 @@ def test_td04ad_static(self):
             nr,a,b,c,d = transform.td04ad(rc,nin,nout,index,den,num)
 
 
-def suite():
-   return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
-
-
 if __name__ == "__main__":
     unittest.main()

slycot/tests/test_mb.py

@@ -0,0 +1,82 @@
+#!/usr/bin/env python
+#
+# test_mb.py - test suite for linear algebra commands
+# bnavigator <[email protected]>, Aug 2019
+
+import unittest
+import numpy as np
+
+from slycot import mb05md, mb05nd
+
+from numpy.testing import assert_allclose
+
+
+class test_mb(unittest.TestCase):
+
+    def test_mb05md(self):
+        """ test_mb05md: verify Matrix exponential with slicot doc example
+        data from http://slicot.org/objects/software/shared/doc/MB05MD.html
+        """
+        A = np.array([[ 0.5,  0.,   2.3, -2.6],
+                      [ 0.,   0.5, -1.4, -0.7],
+                      [ 2.3, -1.4,  0.5,  0.0],
+                      [-2.6, -0.7,  0.0,  0.5]])
+        delta = 1.0
+        Ar_ref = np.array([[  26.8551,  -3.2824,  18.7409, -19.4430],
+                           [  -3.2824,   4.3474,  -5.1848,   0.2700],
+                           [  18.7409,  -5.1848,  15.6012, -11.7228],
+                           [ -19.4430,   0.2700, -11.7228,  15.6012]])
+        Vr_ref = np.array([[-0.7,  0.7,  0.1, -0.1],
+                           [ 0.1, -0.1,  0.7, -0.7],
+                           [ 0.5,  0.5,  0.5,  0.5],
+                           [-0.5, -0.5,  0.5,  0.5]])
+        Yr_ref = np.array([[  -0.0349,   0.0050,   0.0249,  -0.0249],
+                           [  38.2187,  -5.4598,  27.2991, -27.2991],
+                           [   0.0368,   0.2575,   0.1839,   0.1839],
+                           [  -0.7389,  -5.1723,   3.6945,   3.6945]])
+        VAL_ref = np.array([-3., 4., -1., 2.])
+        (Ar, Vr, Yr, VAL) = mb05md(A, delta)
+
+        assert_allclose(Ar, Ar_ref, atol=0.0001)
+
+        # Order of eigenvalues is not guaranteed, so we check them one by one.
+        for i, e in enumerate(VAL):
+            erow = np.ones(VAL.shape)*e
+            i_ref = np.isclose(erow, VAL_ref)
+            self.assertTrue(any(i_ref),
+                            msg="eigenvalue {} not expected".format(e))
+            # Eigenvectors can have different scaling.
+            vr_ref = Vr_ref[:, i_ref]*Vr[0, i]/Vr_ref[0, i_ref][0]
+            assert_allclose(Vr[:, (i,)], vr_ref, atol=0.0001)
+
+        assert_allclose(np.dot(Vr, Yr), np.dot(Vr_ref, Yr_ref), atol=0.0001)
+
+    def test_mb05nd(self):
+        """ test_mb05nd: verify Matrix exponential and integral
+        data from http://slicot.org/objects/software/shared/doc/MB05ND.html
+        """
+        A = np.array([[5.0,   4.0,   3.0,   2.0,   1.0],
+                      [1.0,   6.0,   0.0,   4.0,   3.0],
+                      [2.0,   0.0,   7.0,   6.0,   5.0],
+                      [1.0,   3.0,   1.0,   8.0,   7.0],
+                      [2.0,   5.0,   7.0,   1.0,   9.0]])
+        delta = 0.1
+        F_ref = np.array([[1.8391,  0.9476,  0.7920,  0.8216,  0.7811],
+                          [0.3359,  2.2262,  0.4013,  1.0078,  1.0957],
+                          [0.6335,  0.6776,  2.6933,  1.6155,  1.8502],
+                          [0.4804,  1.1561,  0.9110,  2.7461,  2.0854],
+                          [0.7105,  1.4244,  1.8835,  1.0966,  3.4134]])
+        H_ref = np.array([[0.1347,  0.0352,  0.0284,  0.0272,  0.0231],
+                          [0.0114,  0.1477,  0.0104,  0.0369,  0.0368],
+                          [0.0218,  0.0178,  0.1624,  0.0580,  0.0619],
+                          [0.0152,  0.0385,  0.0267,  0.1660,  0.0732],
+                          [0.0240,  0.0503,  0.0679,  0.0317,  0.1863]])
+
+        (F, H) = mb05nd(A, delta)
+
+        assert_allclose(F, F_ref, atol=0.0001)
+        assert_allclose(H, H_ref, atol=0.0001)
+
+
+if __name__ == "__main__":
+    unittest.main()

slycot/tests/test_mc.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+#
+# test_mc.py - test suite for polynomial and rational function manipulation
+# bnavigator <[email protected]>, Aug 2019
+
+import unittest
+import warnings
+
+from slycot import mc01td
+
+
+class test_mc(unittest.TestCase):
+
+    def test_mc01td(self):
+        """ test_mc01td: doc example
+        data from http://slicot.org/objects/software/shared/doc/MC01TD.html
+        """
+        (dp, stable, nz) = mc01td('C', 4, [2, 0, 1, -1, 1])
+        self.assertEqual(dp, 4)
+        self.assertEqual(stable, 0)
+        self.assertEqual(nz, 2)
+
+    def test_mc01td_D(self):
+        """ test_mc01td_D: test discrete option """
+        (dp, stable, nz) = mc01td('D', 3, [1, 2, 3, 4])
+        self.assertEqual(dp, 3)
+        self.assertEqual(stable, 1)
+        self.assertEqual(nz, 0)
+        (dp, stable, nz) = mc01td('D', 3, [4, 3, 2, 1])
+        self.assertEqual(dp, 3)
+        self.assertEqual(stable, 0)
+        self.assertEqual(nz, 3)
+
+    def test_mc01td_warnings(self):
+        """ test_mc01td_warnings: Test warnings """
+        T = [([0, 0], "entry P(x) is the zero polynomial."),
+             ([0, 1], "P(x) may have zeros very close to stability boundary."),
+             ([1, 0], "The degree of P(x) has been reduced to 0")]
+        for P, m in T:
+            with warnings.catch_warnings(record=True) as w:
+                (dp, stable, nz) = mc01td('C', len(P)-1, P)
+                self.assertEqual(str(w[0].message), m)
+
+
+if __name__ == "__main__":
+    unittest.main()

slycot/transform.py

@@ -649,7 +649,6 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None):
 
     kdcoef = max(index)+1
     if rowcol == 'R':
-        porm = p
         if ucoeff.ndim != 3:
             e = ValueError("The numerator is not a 3D array!")
             e.info = -7
@@ -664,7 +663,6 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None):
             raise e
         out = _wrapper.td04ad_r(m,p,index,dcoeff,ucoeff,n,tol,ldwork)
     elif rowcol == 'C':
-        porm = m
         if ucoeff.ndim != 3:
             e = ValueError("The numerator is not a 3D array!")
             e.info = -7