python-control · repagh · May 20, 2020 · May 19, 2020 · May 19, 2020 · May 19, 2020
diff --git a/slycot/analysis.py b/slycot/analysis.py
@@ -21,19 +21,23 @@
 from .exceptions import raise_if_slycot_error, SlycotParameterError
 
 
-def ab01nd(n,m,A,B,jobz='N',tol=0,ldwork=None):
+def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None):
     """ Ac,Bc,ncont,indcon,nblk,Z,tau = ab01nd_i(n,m,A,B,[jobz,tol,ldwork])
 
     To find a controllable realization for the linear time-invariant
     multi-input system
 
+    ::
+
           dX/dt = A * X + B * U,
 
     where A and B are N-by-N and N-by-M matrices, respectively,
     which are reduced by this routine to orthogonal canonical form
     using (and optionally accumulating) orthogonal similarity
-    transformations.  Specifically, the pair (A, B) is reduced to
-    the pair (Ac, Bc),  Ac = Z' * A * Z,  Bc = Z' * B,  given by
+    transformations.  Specifically, the pair ``(A, B)`` is reduced to
+    the pair ``(Ac, Bc),  Ac = Z' * A * Z,  Bc = Z' * B``,  given by
+
+    ::
 
           [ Acont     *    ]         [ Bcont ]
      Ac = [                ],   Bc = [       ],
@@ -49,88 +53,89 @@ def ab01nd(n,m,A,B,jobz='N',tol=0,ldwork=None):
              [  .   .       .    .     .  ]         [ .  ]
              [  0   0   . . .  Ap,p-1 App ]         [ 0  ]
 
-    where the blocks  B1, A21, ..., Ap,p-1  have full row ranks and
-    p is the controllability index of the pair.  The size of the
-    block  Auncont is equal to the dimension of the uncontrollable
-    subspace of the pair (A, B).
+    where the blocks ``B1, A21, ..., Ap,p-1`` have full row ranks and
+    `p` is the controllability index of the pair.  The size of the
+    block `Auncont` is equal to the dimension of the uncontrollable
+    subspace of the pair ``(A, B)``.
 
-    Required arguments:
-        n : input int
-            The order of the original state-space representation, i.e.
-            the order of the matrix A.  N > 0.
-        m : input int
-            The number of system inputs, or of columns of B.  M > 0.
-        A : rank-2 array('d') with bounds (n,n)
-            The leading n-by-n part of this array must contain the original
-            state dynamics matrix A.
-        B : rank-2 array('d') with bounds (n,m)
-            The leading n-by-m part of this array must contain the input
-            matrix B.
-    Optional arguments:
-        jobz := 'N' input string(len=1)
-            Indicates whether the user wishes to accumulate in a matrix Z
-            the orthogonal similarity transformations for reducing the system,
-            as follows:
-            = 'N':  Do not form Z and do not store the orthogonal transformations;
-            = 'F':  Do not form Z, but store the orthogonal transformations in
-                    the factored form;
-            = 'I':  Z is initialized to the unit matrix and the orthogonal
-                    transformation matrix Z is returned.
-        tol := 0 input float
-            The tolerance to be used in rank determination when transforming
-            (A, B). If tol <= 0 a default value is used.
-        ldwork := max(n,3*m) input int
-            The length of the cache array. ldwork >= max(n, 3*m).
-            For optimum performance it should be larger.
-    Return objects:
-        Ac : rank-2 array('d') with bounds (n,n)
-            The leading ncont-by-ncont part contains the upper block
-            Hessenberg state dynamics matrix Acont in Ac, given by Z'*A*Z,
-            of a controllable realization for the original system. The
-            elements below the first block-subdiagonal are set to zero.
-        Bc : rank-2 array('d') with bounds (n,m)
-            The leading ncont-by-m part of this array contains the transformed
-            input matrix Bcont in Bc, given by Z'*B, with all elements but the
-            first block set to zero.
-        ncont : int
-            The order of the controllable state-space representation.
-        indcon : int
-            The controllability index of the controllable part of the system
-            representation.
-        nblk : rank-1 array('i') with bounds (n)
-            The leading indcon elements of this array contain the the orders of
-            the diagonal blocks of Acont.
-        Z : rank-2 array('d') with bounds (n,n)
-            If jobz = 'I', then the leading N-by-N part of this array contains
-            the matrix of accumulated orthogonal similarity transformations
-            which reduces the given system to orthogonal canonical form.
-            If jobz = 'F', the elements below the diagonal, with the array tau,
-            represent the orthogonal transformation matrix as a product of
-            elementary reflectors. The transformation matrix can then be
-            obtained by calling the LAPACK Library routine DORGQR.
-            If jobz = 'N', the array Z is None.
-        tau : rank-1 array('d') with bounds (n)
-            The elements of tau contain the scalar factors of the
-            elementary reflectors used in the reduction of B and A."""
+    Parameters
+    ----------
+    n : int
+        The order of the original state-space representation, i.e.
+        the order of the matrix A.  ``n > 0``.
+    m : int
+        The number of system inputs, or of columns of B.  ``m > 0``.
+    A : (n,n) array_like
+        The original state dynamics matrix A.
+    B : (n,m) array_like
+        The input matrix B.
+    jobz : {'N', 'F', 'I'}, optional
+        Indicates whether the user wishes to accumulate in a matrix Z
+        the orthogonal similarity transformations for reducing the system,
+        as follows:
+
+        := 'N':  Do not form Z and do not store the orthogonal transformations;
+                 (default)
+        := 'F':  Do not form Z, but store the orthogonal transformations in
+                 the factored form;
+        := 'I':  Z is initialized to the unit matrix and the orthogonal
+                 transformation matrix Z is returned.
+    tol : float, optional
+        The tolerance to be used in rank determination when transforming
+        ``(A, B)``. If ``tol <= 0`` a default value is used.
+    ldwork : int, optional
+        The length of the cache array. ``ldwork >= max(n, 3*m)``.
+        For optimum performance it should be larger.
+        default: ``ldwork = max(n, 3*m)``
+
+    Returns
+    -------
+    Ac : (n,n) ndarray
+        The leading ncont-by-ncont part contains the upper block
+        Hessenberg state dynamics matrix Acont in Ac, given by Z'*A*Z,
+        of a controllable realization for the original system. The
+        elements below the first block-subdiagonal are set to zero.
+    Bc : (n,m) ndarray
+        The leading ncont-by-m part of this array contains the transformed
+        input matrix Bcont in Bc, given by ``Z'*B``, with all elements but the
+        first block set to zero.
+    ncont : int
+        The order of the controllable state-space representation.
+    indcon : int
+        The controllability index of the controllable part of the system
+        representation.
+    nblk : (n,) int ndarray
+        The leading indcon elements of this array contain the the orders of
+        the diagonal blocks of Acont.
+    Z : (n,n) ndarray
+        - If jobz = 'I', then the leading N-by-N part of this array contains
+          the matrix of accumulated orthogonal similarity transformations
+          which reduces the given system to orthogonal canonical form.
+        - If jobz = 'F', the elements below the diagonal, with the array tau,
+          represent the orthogonal transformation matrix as a product of
+          elementary reflectors. The transformation matrix can then be
+          obtained by calling the LAPACK Library routine DORGQR.
+        - If jobz = 'N', the array Z is `None`.
+    tau : (n, ) ndarray
+        The elements of tau contain the scalar factors of the
+        elementary reflectors used in the reduction of B and A.
+    """
 
     hidden = ' (hidden by the wrapper)'
     arg_list = ['jobz', 'n', 'm', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden,
                 'ncont', 'indcon', 'nblk', 'Z', 'LDZ'+hidden, 'tau', 'tol',
                 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden]
 
-    wrappermap = {"N": _wrapper.ab01nd_n,
-                  "I": _wrapper.ab01nd_i,
-                  "F": _wrapper.ab01nd_f}
-
     if ldwork is None:
         ldwork = max(n, 3*m)
 
-    out = wrappermap[jobz](n, m, A, B, tol=tol, ldwork=ldwork)
-    raise_if_slycot_error(out[-1], arg_list)
-    # sets Z to None
+    Ac, Bc, ncont, indcon, nblk, Z, tau, info = _wrapper.ab01nd(
+        jobz, n, m, A, B, tol=tol, ldwork=ldwork)
+    raise_if_slycot_error(info, arg_list)
+
     if jobz == "N":
-        out[5] = None
-    return out[:-1]
+        Z = None
+    return Ac, Bc, ncont, indcon, nblk, Z, tau
 
 
 def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'):

diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
@@ -1,67 +1,24 @@
 !    -*- f90 -*-
-subroutine ab01nd_i(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f
-    fortranname ab01nd
-    character intent(hide) :: jobz = 'I'
+subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f
+    character :: jobz
     integer check(n>0) :: n
-    integer check(n>0) :: m
+    integer check(m>0) :: m
     double precision intent(in,out),dimension(n,n),depend(n) :: a
     integer intent(hide),depend(a) :: lda=shape(a,0)
     double precision intent(in,out),dimension(n,m),depend(n,m) :: b
     integer intent(hide),depend(b) :: ldb=shape(b,0)
     integer intent(out) :: ncont
     integer intent(out) :: indcon
     integer intent(out),dimension(n),depend(n) :: nblk
-    double precision intent(out),dimension(n,n),depend(n) :: z
-    integer intent(hide),depend(z) :: ldz=shape(z,0)
+    double precision intent(out),dimension(ldz,n),depend(n,ldz) :: z
+    integer intent(hide),depend(n, jobz) :: ldz = (*jobz == 'N' ? 1 : n)
     double precision intent(out),dimension(n),depend(n) :: tau
     double precision :: tol = 0
     integer intent(hide,cache),dimension(m),depend(m) :: iwork
     double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
     integer :: ldwork = max(n,3*m)
     integer intent(out) :: info
-end subroutine ab01nd_i
-subroutine ab01nd_f(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f
-    fortranname ab01nd
-    character intent(hide) :: jobz = 'F'
-    integer check(n>0) :: n
-    integer check(n>0) :: m
-    double precision intent(in,out),dimension(n,n),depend(n) :: a
-    integer intent(hide),depend(a) :: lda=shape(a,0)
-    double precision intent(in,out),dimension(n,m),depend(n,m) :: b
-    integer intent(hide),depend(b) :: ldb=shape(b,0)
-    integer intent(out) :: ncont
-    integer intent(out) :: indcon
-    integer intent(out),dimension(n),depend(n) :: nblk
-    double precision intent(out),dimension(n,n),depend(n) :: z
-    integer intent(hide),depend(z) :: ldz=shape(z,0)
-    double precision intent(out),dimension(n),depend(n) :: tau
-    double precision :: tol = 0
-    integer intent(hide,cache),dimension(m),depend(m) :: iwork
-    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
-    integer :: ldwork = max(n,3*m)
-    integer intent(out) :: info
-end subroutine ab01nd_f
-subroutine ab01nd_n(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f
-    fortranname ab01nd
-    character intent(hide) :: jobz = 'N'
-    integer check(n>0) :: n
-    integer check(n>0) :: m
-    double precision intent(in,out),dimension(n,n),depend(n) :: a
-    integer intent(hide),depend(a) :: lda=shape(a,0)
-    double precision intent(in,out),dimension(n,m),depend(n,m) :: b
-    integer intent(hide),depend(b) :: ldb=shape(b,0)
-    integer intent(out) :: ncont
-    integer intent(out) :: indcon
-    integer intent(out),dimension(n),depend(n) :: nblk
-    double precision intent(out),dimension(1,1),depend(n) :: z
-    integer intent(hide),depend(z) :: ldz=shape(z,0)
-    double precision intent(out),dimension(n),depend(n) :: tau
-    double precision :: tol = 0
-    integer intent(hide,cache),dimension(m),depend(m) :: iwork
-    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
-    integer :: ldwork = max(n,3*m)
-    integer intent(out) :: info
-end subroutine ab01nd_n
+end subroutine ab01nd
 subroutine ab05md(uplo,over,n1,m1,p1,n2,p2,a1,lda1,b1,ldb1,c1,ldc1,d1,ldd1,a2,lda2,b2,ldb2,c2,ldc2,d2,ldd2,n,a,lda,b,ldb,c,ldc,d,ldd,dwork,ldwork,info) ! in AB05MD.f
     character :: uplo = 'U'
     character intent(hide) :: over = 'N' ! not sure how the overlap works

diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt
@@ -1,6 +1,7 @@
 set(PYSOURCE
 
   __init__.py
+  test_ab01.py
   test_ab08n.py
   test_ag08bd.py
   test_examples.py

diff --git a/slycot/tests/test_ab01.py b/slycot/tests/test_ab01.py
@@ -0,0 +1,66 @@
+"""
+Test ab01 wrappers
+
+@author: bnavigator
+"""
+
+from numpy import array
+from numpy.testing import assert_allclose, assert_equal
+
+from scipy.linalg.lapack import dorgqr
+
+from slycot.analysis import ab01nd
+
+
+def test_ab01nd():
+    """SLICOT doc example
+
+    http://slicot.org/objects/software/shared/doc/AB01ND.html"""
+
+    # Example program data
+    n = 3
+    m = 2
+    tol = 0.0
+
+    A = array([[-1., 0., 0.],
+               [-2., -2., -2.],
+               [-1., 0., -3.]])
+    B = array([[1., 0.],
+               [0,  2.],
+               [0., 1.]])
+
+    for jobz in ['N', 'I', 'F']:
+        Ac, Bc, ncont, indcon, nblk, Z, tau = ab01nd(n, m, A, B,
+                                                     jobz=jobz, tol=tol)
+
+        # The transformed state dynamics matrix of a controllable realization
+        assert_allclose(Ac[:ncont, :ncont], array([[-3.0000,  2.2361],
+                                                   [ 0.0000, -1.0000]]),
+                        atol=0.0001)
+
+        # and the dimensions of its diagonal blocks are
+        assert_equal(nblk[:indcon], array([2]))
+
+        # The transformed input/state matrix B of a controllable realization
+        assert_allclose(Bc[:ncont, :],array([[ 0.0000, -2.2361],
+                                             [ 1.0000,  0.0000]]),
+                        atol=0.0001)
+
+        # The controllability index of the transformed system representation
+        assert indcon == 1
+
+        if jobz == 'N':
+            assert Z is None
+            continue
+        elif jobz == 'I':
+            Z_ = Z
+        elif jobz == 'F':
+            Z_, _, info = dorgqr(Z, tau)
+            assert info == 0
+
+        # The similarity transformation matrix Z
+        assert_allclose(Z_, array([[ 0.0000,  1.0000,  0.0000],
+                                  [-0.8944,  0.0000, -0.4472],
+                                  [-0.4472,  0.0000,  0.8944]]),
+                        atol=0.0001)
+
diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py
@@ -1,5 +1,5 @@
 # ===================================================
-# ag08bd tests
+# ab08n* tests
 
 import unittest
 from slycot import analysis
@@ -9,7 +9,7 @@
 from numpy.testing import assert_equal, assert_allclose
 
 
-class test_ab08n(unittest.TestCase):
+class test_ab08nX(unittest.TestCase):
     """ Test regular pencil construction ab08nX with input parameters
     according to example in documentation """