fixing submodules

mccullerlp · Dec 16, 2023 · Dec 18, 2023 · Dec 18, 2023 · Dec 18, 2023 · Dec 18, 2023
commit 30b052eaecef7a96d6d5314b89d57bd30c825c2d
diff --git a/.gitmodules b/.gitmodules
@@ -1,3 +1,4 @@
 [submodule "slycot/src/SLICOT-Reference"]
 	path = slycot/src/SLICOT-Reference
 	url = https://github.com/python-control/SLICOT-Reference
+	branch = balancers
diff --git a/slycot/src/SLICOT-Reference b/slycot/src/SLICOT-Reference
diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf
@@ -421,6 +421,29 @@ subroutine tb01pd(job,equil,n,m,p,a,lda,b,ldb,c,ldc,nr,tol,iwork,dwork,ldwork,in
     integer optional,depend(n,m,p) :: ldwork=MAX(1,n+MAX(n,MAX(3*m,3*p)))
     integer intent(out) :: info
 end subroutine tb01pd
+subroutine tg01jd(job,systyp,equil,n,m,p,a,lda,e,lde,b,ldb,c,ldc,nr,infred,tol,iwork,dwork,ldwork,info) ! in TG01JD.f
+  character intent(in) :: job = 'M'
+  character intent(in) :: systyp = 'R'
+  character intent(in) :: equil = 'S'
+  integer required,check(n>=0) :: n
+  integer required,check(m>=0) :: m
+  integer required,check(p>=0) :: p
+  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,n),depend(n) :: e
+  integer, intent(hide), depend(e) :: lde=shape(e,0)
+  double precision intent(in,out,copy),dimension(n,MAX(m,p)),depend(n,m,p) :: b
+  integer, intent(hide), depend(b) :: ldb=shape(b,0)
+  double precision intent(in,out,copy),dimension(MAX(p,m),n),depend(p,n,m) :: c
+  integer, intent(hide), depend(c) :: ldc=shape(c,0)
+  integer, intent(out) :: nr
+  integer, dimension(7) :: infred
+  double precision, intent(in) :: tol = 0.0
+  integer, intent(hide,cache),dimension(2*n+MAX(m,p)),depend(n,m,p) :: iwork
+  double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
+  integer optional,depend(n,m,p) :: ldwork=(2*n*n+n*m+n*p+MAX(n,MAX(2*p,2*p)))
+  integer intent(out) :: info
+end subroutine tg01jd
 subroutine tg01ad(job,l,n,m,p,thresh,a,lda,e,lde,b,ldb,c,ldc,lscale,rscale,dwork,info) ! in TG01AD.f
     character :: job
     integer intent(in),required,check(l>=0) :: l

diff --git a/slycot/transform.py b/slycot/transform.py
@@ -192,6 +192,119 @@ def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None):
     raise_if_slycot_error(out[-1], arg_list)
     return out[:-1]
 
+def tb01jd(n, m, p, A, E, B, C, job='I', systyp='R', equil='S', tol=1e-13, ldwork=None):
+    """Ar, Br, Cr, nr = tb01jd(n,m,p,A,B,C,[job,equil,tol,ldwork])
+
+    To find a reduced (controllable, observable, or minimal) state-
+    space representation (Ar,Er,Br,Cr) for any original state-space
+    representation (A,E,B,C). The matrix Ar,Er is in upper block
+    Hessenberg form. Also removes unobservable/controllable states
+    at infinity.
+
+    Parameters
+    ----------
+    n : int
+        Order of the State-space representation.
+    m : int
+        Number of inputs.
+    p : int
+        Number of outputs.
+    A : (n, n) array_like
+        State dynamics matrix.
+    E : (n, n) array_like
+        State descriptor dynamics Matrix
+    B : (n, max(m,p)) array_like
+        The leading n-by-m part of this array must contain the original
+        input/state matrix B; the remainder of the leading n-by-max(m,p)
+        part is used as internal workspace.
+    C : (p, n) array_like
+        The leading p-by-n part of this array must contain the original
+        state/output matrix C; the remainder of the leading max(1,m,p)-by-n
+        part is used as internal workspace.
+    job : {'I', 'C', 'O'}, optional
+        Indicates whether the user wishes to remove the
+        uncontrollable and/or unobservable parts as follows:
+        = 'I':  Remove both the uncontrollable and unobservable
+                parts to get a minimal state-space representation;
+        = 'C':  Remove the uncontrollable part only to get a
+                controllable state-space representation;
+        = 'O':  Remove the unobservable part only to get an
+                observable state-space representation.
+        Default is 'I'.
+    systyp : {'R', 'S', 'P'}, defaults to 'R'
+        = 'R':  Rational transfer-function matrix; (most general)
+        = 'S':  Proper (standard) transfer-function matrix;
+        = 'P':  Polynomial transfer-function matrix.
+    equil : {'S', 'N'}, optional
+        Specifies whether the user wishes to preliminarily balance
+        the triplet (A,B,C) as follows:
+        = 'S':  Perform balancing (scaling);
+        = 'N':  Do not perform balancing.
+    tol : float, optional
+        The tolerance to be used in rank determination when
+        transforming (A, B, C). If the user sets tol > 0, then
+        the given value of tol is used as a lower bound for the
+        reciprocal condition number.
+        Default is `1e-13`.
+    ldwork : int, optional
+        The length of the cache array. Default is None to use a
+        default, preferred size 2*n*n+n*m+n*p+MAX(n,2*p,2*p)
+
+    Returns
+    -------
+    Ar : (nr, nr) ndarray
+        Contains the upper block Hessenberg state dynamics matrix
+        Ar of a minimal, controllable, or observable realization
+        for the original system, depending on the value of JOB,
+        JOB = 'I', JOB = 'C', or JOB = 'O', respectively.
+    Er : (nr, nr) ndarray
+        Contains the upper block Hessenberg state descriptor
+        dynamics matrix Er of a minimal, controllable, or
+        observable realization for the original system, depending
+         on the value of JOB,
+        JOB = 'I', JOB = 'C', or JOB = 'O', respectively.
+    Br : (nr, m) ndarray
+        Contains the transformed input/state matrix Br of a
+        minimal, controllable, or observable realization for the
+        original system, depending on the value of JOB, JOB = 'I',
+        JOB = 'C', or JOB = 'O', respectively.  If JOB = 'C', only
+        the first IWORK(1) rows of B are nonzero.
+    Cr : (p, nr) ndarray
+        Contains the transformed state/output matrix Cr of a
+        minimal, C controllable, or observable realization for the
+        original C system, depending on the value of JOB, JOB =
+        'I', C JOB = 'C', or JOB = 'O', respectively.  C If JOB =
+        'I', or JOB = 'O', only the last IWORK(1) columns C (in
+        the first NR columns) of C are nonzero.
+    nr : int
+        The order of the reduced state-space representation
+        (Ar,Br,Cr) of a minimal, controllable, or observable
+        realization for the original system, depending on
+        JOB = 'I', JOB = 'C', or JOB = 'O'.
+    """
+    hidden = ' (hidden by the wrapper)'
+    arg_list = ['job', 'equil', 'n','m','p','A','lda'+hidden,'B','ldb'+hidden,
+                'C','ldc'+hidden,'nr','tol','iwork'+hidden,'dwork'+hidden,
+                'ldwork','info'+hidden]
+    if ldwork is None:
+        # the wrapper knows good ldwork
+        ldwork = None
+    elif EQUIL == 'S':
+        assert(ldwork >= max(8*n,2*m,2*p))
+        raise SlycotParameterError("ldwork is too small", -15)
+    elif EQUIL == 'N':
+        assert(ldwork >= max(n,2*m,2*p))
+        raise SlycotParameterError("ldwork is too small", -15)
+    out = _wrapper.tb01jd(n=n,m=m,p=p,a=A,e=E,b=B,c=C,
+                          job=job,systyp=systyp,equil=equil,
+                          tol=tol,ldwork=ldwork)
+
+    # extra information about block sizes returned by the wrapped
+    # function
+    infred = out[5]
+    raise_if_slycot_error(out[-1], arg_list)
+    return out[:5]
+
 def tb03ad(n,m,p,A,B,C,D,leri,equil='N',tol=0.0,ldwork=None):
     """ A_min,b_min,C_min,nr,index,pcoeff,qcoeff,vcoeff = tb03ad_l(n,m,p,A,B,C,D,leri,[equil,tol,ldwork])
 
@@ -1266,6 +1379,7 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld
 
     return A,E,B,C,ranke,rnka22,Q,Z
 
+#TODO clean documentation
 def tg01gd(
         jobs,
         l, n, m, p,
+3 −1		src/MB04DL.f
+3 −1		src/MB04DP.f
+3 −1		src/MB4DLZ.f
+3 −1		src/MB4DPZ.f
+2 −1		src/TB01ID.f
+3 −1		src/TG01AD.f
+3 −1		src/TG01AZ.f