python-control · bnavigator · Sep 24, 2023 · Sep 19, 2023 · Sep 20, 2023 · Sep 21, 2023
diff --git a/slycot/__init__.py b/slycot/__init__.py
@@ -35,7 +35,7 @@
                            ab08nd, ab08nz,
                            ab09ad, ab09ax, ab09bd, ab09md, ab09nd,
                            ab13bd, ab13dd, ab13ed, ab13fd, ab13md)
-    
+
     # Benchmark routines (0/6 wrapped)
 
     # Adaptive control routines (0/0 wrapped)
@@ -46,22 +46,22 @@
 
     # Identification routines (0/15 wrapped)
 
-    # Mathematical routines (7/281 wrapped)
-    from .math import (mb03rd, mb03vd, mb03vy, mb03wd,
+    # Mathematical routines (8/281 wrapped)
+    from .math import (mb02ed, mb03rd, mb03vd, mb03vy, mb03wd,
                        mb05md, mb05nd,
                        mc01td)
 
     # Nonlinear Systems (0/16 wrapped)
 
     # Synthesis routines ((16+1)/131 wrapped), sb03md57 is not part of slicot
     from .synthesis import (sb01bd,
-                            sb02md, sb02mt, sb02od, 
+                            sb02md, sb02mt, sb02od,
                             sb03md, sb03md57, sb03od,
                             sb04md, sb04qd,
                             sb10ad, sb10dd, sb10fd, sb10hd, sb10yd,
                             sg02ad,
                             sg03ad, sg03bd)
-                            
+
     # Transformation routines (10/77 wrapped)
     from .transform import (tb01id, tb01pd,
                             tb03ad,

diff --git a/slycot/math.py b/slycot/math.py
@@ -23,6 +23,107 @@
 import numpy as np
 
 
+def mb02ed(typet: str, T: np.ndarray, B: np.ndarray, n: int, k: int,  nrhs: int):
+    """ X, T = mb02ed(typet, T, B, n, k, nrhs)
+
+    Solve a system of linear equations T*X = B or X*T = B with a positive
+    definite block Toeplitz matrix T.
+
+    Parameters
+    ----------
+    typet: str
+        Specifies the type of T:
+            - 'R': T contains the first block row of an s.p.d. block Toeplitz matrix,
+                and the system X*T = B is solved.
+            - 'C': T contains the first block column of an s.p.d. block Toeplitz matrix,
+                and the system T*X = B is solved.
+        Note: the notation x / y means that x corresponds to
+                typet = 'R' and y corresponds to typet = 'C'.
+    T : array_like
+            The leading k-by-n*k / n*k-by-k part of this array must contain the first
+            block row/column of an s.p.d. block Toeplitz matrix.
+    B : array_like
+            The leading nrhs-by-n*k / n*k-by-nrhs part of this array must contain the
+            right-hand side matrix B.
+    n : int
+            The number of blocks in T. n >= 0.
+    k : int
+       The number of rows/columns in T, equal to the blocksize. k >= 0.
+    nrhs : int
+            The number of right-hand sides. nrhs >= 0.
+
+    Returns
+    -------
+    X : ndarray
+        Leading nrhs-by-n*k / n*k-by-nrhs part of
+        this array contains the solution matrix X.
+    T: ndarray
+        If no error is thrown  and  nrhs > 0,  then the leading
+        k-by-n*k / n*k-by-k part of this array contains the last
+        row / column of the Cholesky factor of inv(T).
+
+    Raises
+    ------
+    SlycotArithmeticError
+        :info = 1:
+            The reduction algorithm failed. The Toeplitz matrix associated
+            with T is not numerically positive definite.
+    SlycotParameterError
+        :info = -1:
+            typet must be either "R" or "C"
+        :info = -2:
+            k must be >= 0
+        :info = -3:
+            n must be >= 0
+        :info = -4:
+            nrhs must be >= 0
+
+    Notes
+    -----
+    The algorithm uses Householder transformations, modified hyperbolic rotations,
+    and block Gaussian eliminations in the Schur algorithm [1], [2].
+
+    References
+    ----------
+    [1] Kailath, T. and Sayed, A.
+        Fast Reliable Algorithms for Matrices with Structure.
+        SIAM Publications, Philadelphia, 1999.
+
+    [2] Kressner, D. and Van Dooren, P.
+        Factorizations and linear system solvers for matrices with Toeplitz structure.
+        SLICOT Working Note 2000-2, 2000.
+
+    Numerical Aspects
+    -----------------
+    The implemented method is numerically equivalent to forming the Cholesky factor R and the
+    inverse Cholesky factor of T using the generalized Schur algorithm and solving the systems
+    of equations R*X = L*B or X*R = B*L by a blocked backward substitution algorithm.
+    The algorithm requires O(K * N^2 + K * N * NRHS) floating-point operations.
+
+    """
+
+    hidden = " (hidden by the wrapper)"
+    arg_list = [
+        "typet",
+        "k",
+        "n",
+        "nrhs",
+        "t",
+        "ldt" + hidden,
+        "b",
+        "ldb" + hidden,
+        "ldwork" + hidden,
+        "dwork" + hidden,
+        "info",
+    ]
+
+    T, X, info = _wrapper.mb02ed(typet=typet, k=k, n=n, nrhs=nrhs, t=T, b=B)
+
+    raise_if_slycot_error(info, arg_list, docstring=mb02ed.__doc__, checkvars=locals())
+
+    return X, T
+
+
 def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0):
     """Ar, Xr, blsize, W = mb03rd(n, A, [X, jobx, sort, pmax, tol])
 

diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf
@@ -12,6 +12,20 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f
 	integer intent(out) :: info
 end subroutine mc01td
 
+subroutine mb02ed(typet,k,n,nrhs,t,ldt,b,ldb,dwork,ldwork,info) ! in MB02ED.f
+  character :: typet
+  integer intent(in),required  :: k
+  integer intent(in),required :: n
+  integer intent(in),required :: nrhs
+  double precision intent(in,out,copy),dimension(ldt,*) :: t
+  integer, intent(hide),optional,check(shape(t, 0) == ldt),depend(t) :: ldt=shape(t, 0)
+  double precision intent(in,out,copy),dimension(ldb,*) :: b
+  integer, intent(hide),optional,check(shape(b, 0) == ldb),depend(b) :: ldb=shape(b, 0)
+  double precision intent(cache,hide),dimension(ldwork) :: dwork
+  integer optional,check(ldwork>=n*k*k+(n+2)*k), depend(n,k) :: ldwork=max(1,n*k*k+(n+2)*k)
+  integer intent(out):: info
+end subroutine mb02ed
+
 subroutine mb03rd(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
   character intent(in) :: jobx
   character intent(in),required :: sort

diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py
@@ -9,12 +9,181 @@
 from numpy.testing import assert_allclose
 from scipy.linalg import schur
 
-from slycot import math, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
-from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning
+from slycot import math, mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
+from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning, SlycotParameterError
 
 from .test_exceptions import assert_docstring_parse
 
 
+def test_mb02ed_ex():
+    """Test MB02ED using the example given in the MB02ED SLICOT Documentation"""
+    n = 3
+    k = 3
+    nrhs = 2
+    TYPET = "C"
+    T = np.array(
+        [
+            [3.0000, 1.0000, 0.2000],
+            [1.0000, 4.0000, 0.4000],
+            [0.2000, 0.4000, 5.0000],
+            [0.1000, 0.1000, 0.2000],
+            [0.2000, 0.0400, 0.0300],
+            [0.0500, 0.2000, 0.1000],
+            [0.1000, 0.0300, 0.1000],
+            [0.0400, 0.0200, 0.2000],
+            [0.0100, 0.0300, 0.0200],
+        ]
+    )
+    B = np.array(
+        [
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+        ]
+    )
+    X = np.array(
+        [
+            [0.2408, 0.4816],
+            [0.1558, 0.3116],
+            [0.1534, 0.3068],
+            [0.2302, 0.4603],
+            [0.1467, 0.2934],
+            [0.1537, 0.3075],
+            [0.2349, 0.4698],
+            [0.1498, 0.2995],
+            [0.1653, 0.3307],
+        ]
+    )
+
+    result,_ = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs)
+    np.testing.assert_almost_equal(result, X, decimal=4)
+
+def test_mb02ed_parameter_errors():
+    """Test for errors in the input parameters of MB02ED"""
+    n = 3
+    k = 3
+    nrhs = 2
+    TYPET = "C"
+    T = np.array(
+        [
+            [3.0000, 1.0000, 0.2000],
+            [1.0000, 4.0000, 0.4000],
+            [0.2000, 0.4000, 5.0000],
+            [0.1000, 0.1000, 0.2000],
+            [0.2000, 0.0400, 0.0300],
+            [0.0500, 0.2000, 0.1000],
+            [0.1000, 0.0300, 0.1000],
+            [0.0400, 0.0200, 0.2000],
+            [0.0100, 0.0300, 0.0200],
+        ]
+    )
+    B = np.array(
+        [
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+        ]
+    )
+    X = np.array(
+        [
+            [0.2408, 0.4816],
+            [0.1558, 0.3116],
+            [0.1534, 0.3068],
+            [0.2302, 0.4603],
+            [0.1467, 0.2934],
+            [0.1537, 0.3075],
+            [0.2349, 0.4698],
+            [0.1498, 0.2995],
+            [0.1653, 0.3307],
+        ]
+    )
+
+    # Test for wrong parameter typet
+    with pytest.raises(expected_exception=SlycotParameterError, match='typet must be either "R" or "C"') as cm:
+        mb02ed(T=T, B=B, n=n, k=k, typet='U', nrhs=nrhs)
+        assert cm.value.info == -1
+    #Test for negative number of columns
+    with pytest.raises(expected_exception=SlycotParameterError, match="k must be >= 0") as cm:
+        mb02ed(T=T, B=B, n=n, k=-1, typet=TYPET, nrhs=nrhs)
+        assert cm.value.info == -2
+    #Test for negative number of blocks
+    with pytest.raises(expected_exception=SlycotParameterError, match="n must be >= 0") as cm:
+        mb02ed(T=T, B=B, n=-1, k=k, typet=TYPET, nrhs=nrhs)
+        assert cm.value.info == -3
+    #Test for negative number of right hand sides
+    with pytest.raises(expected_exception=SlycotParameterError, match="nrhs must be >= 0") as cm:
+        mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=-1)
+        assert cm.value.info == -4
+
+
+def test_mb02ed_matrix_error():
+    """Test for a negative definite input matrix in MB02ED"""
+    n = 3
+    k = 3
+    nrhs = 2
+    TYPET = "C"
+    T = np.array(
+        [
+            [3.0000, 1.0000, 0.2000],
+            [1.0000, 4.0000, 0.4000],
+            [0.2000, 0.4000, 5.0000],
+            [0.1000, 0.1000, 0.2000],
+            [0.2000, 0.0400, 0.0300],
+            [0.0500, 0.2000, 0.1000],
+            [0.1000, 0.0300, 0.1000],
+            [0.0400, 0.0200, 0.2000],
+            [0.0100, 0.0300, 0.0200],
+        ]
+    )
+    # Create a negative definite matrix
+    T = -1 * T
+    B = np.array(
+        [
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+            [1.0000, 2.0000],
+        ]
+    )
+    X = np.array(
+        [
+            [0.2408, 0.4816],
+            [0.1558, 0.3116],
+            [0.1534, 0.3068],
+            [0.2302, 0.4603],
+            [0.1467, 0.2934],
+            [0.1537, 0.3075],
+            [0.2349, 0.4698],
+            [0.1498, 0.2995],
+            [0.1653, 0.3307],
+        ]
+    )
+
+    with pytest.raises(SlycotArithmeticError,
+                       match = "The reduction algorithm failed. "
+                       "The Toeplitz matrix associated\nwith T "
+                       r"is not numerically positive definite.") as cm:
+        mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs)
+        assert cm.value.info == 1
+
+
 def test_mb03rd():
     """ Test for Schur form reduction.