Add ab04md wrapper

python-control · Jul 20, 2023 · a60cee4 · a60cee4
1 parent 294eae0
commit a60cee4
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 1 deletion.
diff --git a/slycot/__init__.py b/slycot/__init__.py
@@ -13,7 +13,7 @@
     # import slycot.examples
 
     # Analysis routines (15/40 wrapped)
-    from .analysis import ab01nd, ab05md, ab05nd, ab07nd, ab08nd, ab08nz
+    from .analysis import ab01nd, ab04md, ab05md, ab05nd, ab07nd, ab08nd, ab08nz
     from .analysis import ab09ad, ab09ax, ab09bd, ab09md, ab09nd
     from .analysis import ab13bd, ab13dd, ab13ed, ab13fd, ab13md
 

diff --git a/slycot/analysis.py b/slycot/analysis.py
@@ -139,6 +139,23 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None):
         Z = None
     return Ac, Bc, ncont, indcon, nblk, Z, tau
 
+def ab04md(type_bn, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None):
+    """ a,b,c,d = ab04md()
+    """
+
+    hidden = ' (hidden by the wrapper)'  
+    arg_list = ['type_bn', 'n', 'm', 'p', 'alpha', 'beta', 
+                'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden,
+                'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden]
+
+    if ldwork is None:
+        ldwork = max(n, 3*m)
+
+    out = _wrapper.ab04md(type_bn, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork)
+    info=out[-1]
+    raise_if_slycot_error(info, arg_list)
+
+    return out[:-1]
 
 def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'):
     """ n,a,b,c,d = ab05md(n1,m1,p1,n2,p2,a1,b1,c1,d1,a2,b2,c2,d2,[uplo])

diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
@@ -19,6 +19,26 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo
     integer :: ldwork = max(n,3*m)
     integer intent(out) :: info
 end subroutine ab01nd
+subroutine ab04md(type_bn,n,m,p,alpha,beta,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in AB04MD.f
+    character :: type_bn
+    integer check(n>=0) :: n
+    integer check(m>=0) :: m
+    integer check(p>=0) :: p
+    double precision intent(in) :: alpha
+    double precision intent(in) :: beta
+    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)
+    integer intent(hide,cache),dimension(n),depend(n) :: iwork
+    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
+    integer :: ldwork = max(1,n)
+    integer intent(out) :: info
+end subroutine ab04md
 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