Implement MB02ED

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,

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])
 

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

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.