python-control · mdclemen · Aug 13, 2016 · Aug 15, 2016 · Aug 16, 2016 · Aug 17, 2016
diff --git a/README.rst b/README.rst
@@ -1,10 +1,10 @@
 Slycot
 =============
 
-.. image:: https://travis-ci.org/jgoppert/Slycot.svg
-        :target: https://travis-ci.org/jgoppert/Slycot
-.. image:: https://coveralls.io/repos/jgoppert/Slycot/badge.png
-        :target: https://coveralls.io/r/jgoppert/Slycot
+.. image:: https://travis-ci.org/python-control/Slycot.svg?branch=master
+        :target: https://travis-ci.org/python-control/Slycot
+.. image:: https://coveralls.io/repos/python-control/Slycot/badge.png
+        :target: https://coveralls.io/r/python-control/Slycot
 
 Python wrapper for selected SLICOT routines, notably including solvers for
 Riccati, Lyapunov and Sylvester equations.

diff --git a/slycot/__init__.py b/slycot/__init__.py
@@ -25,8 +25,8 @@
     # Mathematical routines (3/81 wrapped)
     from .math import mc01td, mb05md, mb05nd
 
-    # Synthesis routines (11/50 wrapped)
-    from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md
+    # Synthesis routines (12/50 wrapped)
+    from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od
     from .synthesis import sb04md,sb04qd,sb10ad,sb10hd,sg03ad
     from .synthesis import sg02ad
 

diff --git a/slycot/analysis.py b/slycot/analysis.py
@@ -876,7 +876,7 @@ def ab09bd(dico,job,equil,n,m,p,A,B,C,D,nr=None,tol1=0,tol2=0,ldwork=None):
         e = ArithmeticError('The computation of Hankel singular values failed')
         e.info = out[-1]
         raise e
-    nr,A,B,C,D,hsv = out[:-2]
+    Nr,A,B,C,D,hsv = out[:-2]
     return nr, A[:Nr,:Nr], B[:Nr,:], C[:,:Nr],D[:,:], hsv
 
 # to be replaced by python wrappers
diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
@@ -231,4 +231,30 @@ subroutine ab09ax(dico,job,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,hsv,t,ldt,ti,ldti,t
     integer optional :: ldwork = max(1,n*(max(n,max(m,p))+5)+n*(n+1)/2)
     integer intent(out) :: iwarn
     integer intent(out) :: info
-end subroutine ab09ax
+end subroutine ab09ax
+subroutine ab09bd(dico,job,equil,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,d,ldd,hsv,tol1,tol2,iwork,dwork,ldwork,iwarn,info) !in :balred:AB09BD.f
+    character intent(in) :: dico
+    character intent(in) :: job
+    character intent(in) :: equil
+    character intent(in) :: ordsel
+    integer check(n>=0) :: n
+    integer check(m>=0) :: m
+    integer check(p>=0) :: p
+    integer intent(in,out) :: nr
+    double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
+    integer intent(hide),depend(a) :: lda=shape(a,0)
+    double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b
+    integer intent(hide),depend(b) :: ldb=shape(b,0)
+    double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c
+    integer intent(hide),depend(c) :: ldc=shape(c,0)
+    double precision intent(in,out,copy),dimension(p,m),depend(m,p) :: d
+    integer intent(hide),depend(d) :: ldd=shape(d,0)
+    double precision intent(out),dimension(n),depend(n) :: hsv
+    double precision :: tol1 =0.0
+    double precision :: tol2 =0.0
+    integer intent(hide,cache),dimension(max(m,p)) :: iwork
+    double precision intent(hide,cache),dimension(ldwork) :: dwork
+    integer optional :: ldwork = max(1,n*(2*n+max(n,max(m,p))+5)+n*(n+1)/2)
+    integer intent(out) :: iwarn
+    integer intent(out) :: info
+end subroutine ab09bd
diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf
@@ -347,6 +347,34 @@ subroutine sb03md(dico,job,fact,trana,n,c,ldc,a,lda,u,ldu,scale,sep,ferr,wr,wi,i
     integer optional,check(ldwork>=max(2*n*n,3*n)),depend(n) :: ldwork=max(2*n*n,3*n)
     integer intent(out) :: info
 end subroutine sb03md
+<<<<<<< HEAD
+<<<<<<< HEAD
+subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info) ! in SB03OD.f
+=======
+subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info)
+>>>>>>> f1b51c5... Fixed sb03od bug in synthesis.pyf.
+=======
+subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info) ! in SB03OD.f
+>>>>>>> 47a3169... Added wrapper to SLICOT routine SB03OD.
+    fortranname sb03od
+    character :: dico
+    character :: fact='N'
+    character :: trans='N'
+    integer check(n>0) :: n
+    integer check(m>0) :: m
+    double precision dimension(n,n),depend(n) :: a
+    integer intent(hide),depend(a) :: lda=shape(a,0)
+    double precision dimension(n,n),depend(n) :: q
+    integer intent(hide),depend(q) :: ldq=shape(q,0)
+    double precision intent(in,out,copy),dimension(n,n),depend(n) :: b
+    integer intent(hide),depend(b) :: ldb=shape(b,0)
+    double precision intent(out) :: scale
+    double precision intent(out),dimension(n),depend(n) :: wr
+    double precision intent(out),dimension(n),depend(n) :: wi
+    double precision intent(hide,cache),dimension(ldwork) :: dwork
+    integer optional,check(ldwork>=max(1,4*n + min(m,n))),depend(n,m) :: ldwork=max(1,4*n + min(m,n))
+    integer intent(out) :: info    
+end subroutine sb03od
 subroutine sb04md(n,m,a,lda,b,ldb,c,ldc,z,ldz,iwork,dwork,ldwork,info) ! in SB04MD.f
     integer check(n>0) :: n
     integer check(m>0) :: m

diff --git a/slycot/synthesis.py b/slycot/synthesis.py
@@ -881,6 +881,252 @@ def sb03md(n,C,A,U,dico,job='X',fact='N',trana='N',ldwork=None):
     w.imag = wi[0:n]
     return X,scale,sep,ferr,w
 
+def sb03od(n,m,A,Q,B,dico,fact='N',trans='N',ldwork=None):
+    """  U,scale,w = sb03od(dico,n,m,A,Q,B,[fact,trans,ldwork])
+
+    To solve for X = op(U)'*op(U) either the stable non-negative
+    definite continuous-time Lyapunov equation
+                                 2
+      op(A)'*X + X*op(A) = -scale *op(B)'*op(B)                   (1)
+
+    or the convergent non-negative definite discrete-time Lyapunov
+    equation
+                                 2
+      op(A)'*X*op(A) - X = -scale *op(B)'*op(B)                   (2)
+
+    where op(K) = K or K' (i.e., the transpose of the matrix K), A is
+    an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper
+    triangular matrix containing the Cholesky factor of the solution
+    matrix X, X = op(U)'*op(U), and scale is an output scale factor,
+    set less than or equal to 1 to avoid overflow in X. If matrix B
+    has full rank then the solution matrix X will be positive-definite
+    and hence the Cholesky factor U will be nonsingular, but if B is
+    rank deficient then X may be only positive semi-definite and U
+    will be singular.
+
+    In the case of equation (1) the matrix A must be stable (that
+    is, all the eigenvalues of A must have negative real parts),
+    and for equation (2) the matrix A must be convergent (that is,
+    all the eigenvalues of A must lie inside the unit circle).
+
+    Required arguments
+    ------------------
+
+        n : input int
+            The order of the matrix A and the number of columns in
+            matrix op(B).  n >= 0.
+        m : input int
+            The number of rows in matrix op(B).  m >= 0.
+        A : input rank-2 array('d'), shape  (n,n)
+            On entry, the leading n-by-n part of this array must
+            contain the matrix A. If fact = 'F', then A contains
+            an upper quasi-triangular matrix S in Schur canonical
+            form; the elements below the upper Hessenberg part of the
+            array A are not referenced.
+            On exit, the leading n-by-n upper Hessenberg part of this
+            array contains the upper quasi-triangular matrix S in
+            Schur canonical form from the Shur factorization of A.
+            The contents of array A is not modified if fact = 'F'.
+        Q : input rank-2 array('d'), shape  (n,n)
+            On entry, if fact = 'F', then the leading n-by-n part of
+            this array must contain the orthogonal matrix Q of the
+            Schur factorization of A.
+            Otherwise, Q need not be set on entry.
+            On exit, the leading n-by-n part of this array contains
+            the orthogonal matrix Q of the Schur factorization of A.
+            The contents of array Q is not modified if fact = 'F'.
+        B : input rank-2 array('d'), shape  (m,n)
+            On entry, if trans = 'N', the leading m-by-n part of this
+            array must contain the coefficient matrix B of the
+            equation.
+            On entry, if trans = 'T', the leading N-by-m part of this
+            array must contain the coefficient matrix B of the
+            equation.
+            On exit, the leading n-by-n part of this array contains
+            the upper triangular Cholesky factor U of the solution
+            matrix X of the problem, X = op(U)'*op(U).
+            If m = 0 and n > 0, then U is set to zero.
+        dico : input string(len=1)
+            Specifies the type of Lyapunov equation to be solved as
+            follows:
+            = 'C':  Equation (1), continuous-time case;
+            = 'D':  Equation (2), discrete-time case.
+
+    Optional arguments
+    ------------------
+
+        fact := 'N' input string(len=1)
+            Specifies whether or not the real Schur factorization
+            of the matrix A is supplied on entry, as follows:
+            = 'F':  On entry, A and Q contain the factors from the
+            real Schur factorization of the matrix A;
+            = 'N':  The Schur factorization of A will be computed
+            and the factors will be stored in A and Q.
+        trans := 'N' input string(len=1)
+            Specifies the form of op(K) to be used, as follows:
+            = 'N':  op(K) = K    (No transpose);
+            = 'T':  op(K) = K**T (Transpose).
+        ldwork := None input int
+            The length of the array DWORK.
+            If m > 0, ldwork >= max(1,4*n + min(m,n));
+            If m = 0, ldwork >= 1.
+            For optimum performance ldwork should sometimes be larger.
+
+    Return objects
+    ______________
+
+        U : rank-2 array('d'), shape  (n,n)
+            The leading n-by-n part of this array contains
+            the upper triangular Cholesky factor U of the solution
+            matrix X of the problem, X = op(U)'*op(U).
+        scale : float
+            The scale factor, scale, set less than or equal to 1 to
+            prevent the solution overflowing.
+        w : rank-1 array('c'), shape (n)
+            If fact = 'N', this array contains the eigenvalues of A.
+
+    Raises
+    ______
+
+        ValueError : e
+            e.info contains information about the exact type of exception
+             < 0:  if info = -i, the i-th argument had an illegal value;
+             = 1:  if the Lyapunov equation is (nearly) singular
+                (warning indicator);
+                if DICO = 'C' this means that while the matrix A
+                (or the factor S) has computed eigenvalues with
+                negative real parts, it is only just stable in the
+                sense that small perturbations in A can make one or
+                more of the eigenvalues have a non-negative real
+                part;
+                if DICO = 'D' this means that while the matrix A
+                (or the factor S) has computed eigenvalues inside
+                the unit circle, it is nevertheless only just
+                convergent, in the sense that small perturbations
+                in A can make one or more of the eigenvalues lie
+                outside the unit circle;
+                perturbed values were used to solve the equation;
+             = 2:  if FACT = 'N' and DICO = 'C', but the matrix A is
+                not stable (that is, one or more of the eigenvalues
+                of A has a non-negative real part), or DICO = 'D',
+                but the matrix A is not convergent (that is, one or
+                more of the eigenvalues of A lies outside the unit
+                circle); however, A will still have been factored
+                and the eigenvalues of A returned in WR and WI.
+             = 3:  if FACT = 'F' and DICO = 'C', but the Schur factor S
+                supplied in the array A is not stable (that is, one
+                or more of the eigenvalues of S has a non-negative
+                real part), or DICO = 'D', but the Schur factor S
+                supplied in the array A is not convergent (that is,
+                one or more of the eigenvalues of S lies outside the
+                unit circle);
+             = 4:  if FACT = 'F' and the Schur factor S supplied in
+                the array A has two or more consecutive non-zero
+                elements on the first sub-diagonal, so that there is
+                a block larger than 2-by-2 on the diagonal;
+             = 5:  if FACT = 'F' and the Schur factor S supplied in
+                the array A has a 2-by-2 diagonal block with real
+                eigenvalues instead of a complex conjugate pair;
+             = 6:  if FACT = 'N' and the LAPACK Library routine DGEES
+                has failed to converge. This failure is not likely
+                to occur. The matrix B will be unaltered but A will
+                be destroyed.
+    """
+    hidden = ' (hidden by the wrapper)'
+    arg_list = ['dico','fact', 'trans', 'n', 'm', 'a', 'lda'+hidden, 'q',
+        'ldq'+hidden, 'b', 'ldb'+hidden, 'scale', 'wr'+hidden,
+        'wi'+hidden, 'dwork'+hidden, 'ldwork', 'info'+hidden]
+    if ldwork is None:
+        if m > 0:
+            ldwork = max(1,4*n + min(m,n))
+        elif m == 0:
+            ldwork = 1
+    if dico != 'C' and dico != 'D':
+        raise ValueError('dico must be either D or C')
+    out = _wrapper.sb03od(dico,n,m,A,Q,B,fact=fact,trans=trans,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 = out[-1]
+        raise e
+    if out[-1] == 1:
+        if dico == 'D':
+            error_text = """this means that while the matrix A
+                (or the factor S) has computed eigenvalues inside
+                the unit circle, it is nevertheless only just
+                convergent, in the sense that small perturbations
+                in A can make one or more of the eigenvalues lie
+                outside the unit circle;
+                perturbed values were used to solve the equation;"""
+        else:
+            error_text = """this means that while the matrix A
+                (or the factor S) has computed eigenvalues with
+                negative real parts, it is only just stable in the
+                sense that small perturbations in A can make one or
+                more of the eigenvalues have a non-negative real
+                part;"""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if out[-1] == 2:
+        if dico == 'D':
+            error_text = """the matrix A is not convergent (that is, one or
+                more of the eigenvalues of A lies outside the unit
+                circle); however, A will still have been factored
+                and the eigenvalues of A returned in WR and WI."""
+        else:
+            error_text = """the matrix A is
+                not stable (that is, one or more of the eigenvalues
+                of A has a non-negative real part)."""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if out[-1] == 3:
+        if dico == 'D':
+            error_text = """the Schur factor S
+                supplied in the array A is not convergent (that is,
+                one or more of the eigenvalues of S lies outside the
+                unit circle)."""
+        else:
+            error_text = """the Schur factor S
+                supplied in the array A is not stable (that is, one
+                or more of the eigenvalues of S has a non-negative
+                real part)."""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if out[-1] == 4:
+        if fact == 'F':
+            error_text = """the Schur factor S supplied in
+                the array A has two or more consecutive non-zero
+                elements on the first sub-diagonal, so that there is
+                a block larger than 2-by-2 on the diagonal."""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if out[-1] == 5:
+        if fact == 'F':
+            error_text = """the Schur factor S supplied in
+                the array A has a 2-by-2 diagonal block with real
+                eigenvalues instead of a complex conjugate pair."""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if out[-1] == 6:
+        if fact == 'N':
+            error_text = """the LAPACK Library routine DGEES
+                has failed to converge. This failure is not likely
+                to occur. The matrix B will be unaltered but A will
+                be destroyed."""
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    U,scale,wr,wi = out[:-1]
+    w = _np.zeros(n,'complex64')
+    w.real = wr[0:n]
+    w.imag = wi[0:n]
+    return U,scale,w
+
 def sb04md(n,m,A,B,C,ldwork=None):
     """X = sb04md(n,m,A,B,C[,ldwork])