From a60cee415b096ba09f93338388bbda064581805a Mon Sep 17 00:00:00 2001
From: Johannes Kaisinger <johannes.kaisinger@gmail.com>
Date: Thu, 20 Jul 2023 22:29:00 +0200
Subject: [PATCH] Add ab04md wrapper

---
 slycot/__init__.py      |  2 +-
 slycot/analysis.py      | 17 +++++++++++++++++
 slycot/src/analysis.pyf | 20 ++++++++++++++++++++
 3 files changed, 38 insertions(+), 1 deletion(-)

diff --git a/slycot/__init__.py b/slycot/__init__.py
index 6bd59014..fb1b413b 100644
--- 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
index 32050236..d2edb67e 100644
--- 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
index 633ff6ec..28629e38 100644
--- 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