Skip to content

Commit

Permalink
Implement MB02ED (#214)
Browse files Browse the repository at this point in the history
  • Loading branch information
@saasaa
saasaa committed Sep 24, 2023
1 parent 15d2a69 commit 9485d94
Show file tree
Hide file tree
Showing 4 changed files with 291 additions and 7 deletions.
10 changes: 5 additions & 5 deletions slycot/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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,
Expand Down
101 changes: 101 additions & 0 deletions slycot/math.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down
14 changes: 14 additions & 0 deletions slycot/src/math.pyf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
173 changes: 171 additions & 2 deletions slycot/tests/test_mb.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down

0 comments on commit 9485d94

Please sign in to comment.