Skip to content

Commit 9485d94

Browse files
saasaasaasaa
authored
Implement MB02ED (#214)
1 parent 15d2a69 commit 9485d94

4 files changed

Lines changed: 291 additions & 7 deletions

File tree

‎slycot/__init__.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
ab08nd, ab08nz,
3636
ab09ad, ab09ax, ab09bd, ab09md, ab09nd,
3737
ab13bd, ab13dd, ab13ed, ab13fd, ab13md)
38-
38+
3939
# Benchmark routines (0/6 wrapped)
4040

4141
# Adaptive control routines (0/0 wrapped)
@@ -46,22 +46,22 @@
4646

4747
# Identification routines (0/15 wrapped)
4848

49-
# Mathematical routines (7/281 wrapped)
50-
from .math import (mb03rd, mb03vd, mb03vy, mb03wd,
49+
# Mathematical routines (8/281 wrapped)
50+
from .math import (mb02ed, mb03rd, mb03vd, mb03vy, mb03wd,
5151
mb05md, mb05nd,
5252
mc01td)
5353

5454
# Nonlinear Systems (0/16 wrapped)
5555

5656
# Synthesis routines ((16+1)/131 wrapped), sb03md57 is not part of slicot
5757
from .synthesis import (sb01bd,
58-
sb02md, sb02mt, sb02od,
58+
sb02md, sb02mt, sb02od,
5959
sb03md, sb03md57, sb03od,
6060
sb04md, sb04qd,
6161
sb10ad, sb10dd, sb10fd, sb10hd, sb10yd,
6262
sg02ad,
6363
sg03ad, sg03bd)
64-
64+
6565
# Transformation routines (10/77 wrapped)
6666
from .transform import (tb01id, tb01pd,
6767
tb03ad,

‎slycot/math.py‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,107 @@
2323
import numpy as np
2424

2525

26+
def mb02ed(typet: str, T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int):
27+
""" X, T = mb02ed(typet, T, B, n, k, nrhs)
28+
29+
Solve a system of linear equations T*X = B or X*T = B with a positive
30+
definite block Toeplitz matrix T.
31+
32+
Parameters
33+
----------
34+
typet: str
35+
Specifies the type of T:
36+
- 'R': T contains the first block row of an s.p.d. block Toeplitz matrix,
37+
and the system X*T = B is solved.
38+
- 'C': T contains the first block column of an s.p.d. block Toeplitz matrix,
39+
and the system T*X = B is solved.
40+
Note: the notation x / y means that x corresponds to
41+
typet = 'R' and y corresponds to typet = 'C'.
42+
T : array_like
43+
The leading k-by-n*k / n*k-by-k part of this array must contain the first
44+
block row/column of an s.p.d. block Toeplitz matrix.
45+
B : array_like
46+
The leading nrhs-by-n*k / n*k-by-nrhs part of this array must contain the
47+
right-hand side matrix B.
48+
n : int
49+
The number of blocks in T. n >= 0.
50+
k : int
51+
The number of rows/columns in T, equal to the blocksize. k >= 0.
52+
nrhs : int
53+
The number of right-hand sides. nrhs >= 0.
54+
55+
Returns
56+
-------
57+
X : ndarray
58+
Leading nrhs-by-n*k / n*k-by-nrhs part of
59+
this array contains the solution matrix X.
60+
T: ndarray
61+
If no error is thrown and nrhs > 0, then the leading
62+
k-by-n*k / n*k-by-k part of this array contains the last
63+
row / column of the Cholesky factor of inv(T).
64+
65+
Raises
66+
------
67+
SlycotArithmeticError
68+
:info = 1:
69+
The reduction algorithm failed. The Toeplitz matrix associated
70+
with T is not numerically positive definite.
71+
SlycotParameterError
72+
:info = -1:
73+
typet must be either "R" or "C"
74+
:info = -2:
75+
k must be >= 0
76+
:info = -3:
77+
n must be >= 0
78+
:info = -4:
79+
nrhs must be >= 0
80+
81+
Notes
82+
-----
83+
The algorithm uses Householder transformations, modified hyperbolic rotations,
84+
and block Gaussian eliminations in the Schur algorithm [1], [2].
85+
86+
References
87+
----------
88+
[1] Kailath, T. and Sayed, A.
89+
Fast Reliable Algorithms for Matrices with Structure.
90+
SIAM Publications, Philadelphia, 1999.
91+
92+
[2] Kressner, D. and Van Dooren, P.
93+
Factorizations and linear system solvers for matrices with Toeplitz structure.
94+
SLICOT Working Note 2000-2, 2000.
95+
96+
Numerical Aspects
97+
-----------------
98+
The implemented method is numerically equivalent to forming the Cholesky factor R and the
99+
inverse Cholesky factor of T using the generalized Schur algorithm and solving the systems
100+
of equations R*X = L*B or X*R = B*L by a blocked backward substitution algorithm.
101+
The algorithm requires O(K * N^2 + K * N * NRHS) floating-point operations.
102+
103+
"""
104+
105+
hidden = " (hidden by the wrapper)"
106+
arg_list = [
107+
"typet",
108+
"k",
109+
"n",
110+
"nrhs",
111+
"t",
112+
"ldt" + hidden,
113+
"b",
114+
"ldb" + hidden,
115+
"ldwork" + hidden,
116+
"dwork" + hidden,
117+
"info",
118+
]
119+
120+
T, X, info = _wrapper.mb02ed(typet=typet, k=k, n=n, nrhs=nrhs, t=T, b=B)
121+
122+
raise_if_slycot_error(info, arg_list, docstring=mb02ed.__doc__, checkvars=locals())
123+
124+
return X, T
125+
126+
26127
def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0):
27128
"""Ar, Xr, blsize, W = mb03rd(n, A, [X, jobx, sort, pmax, tol])
28129

‎slycot/src/math.pyf‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,20 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f
1212
integer intent(out) :: info
1313
end subroutine mc01td
1414

15+
subroutine mb02ed(typet,k,n,nrhs,t,ldt,b,ldb,dwork,ldwork,info) ! in MB02ED.f
16+
character :: typet
17+
integer intent(in),required :: k
18+
integer intent(in),required :: n
19+
integer intent(in),required :: nrhs
20+
double precision intent(in,out,copy),dimension(ldt,*) :: t
21+
integer, intent(hide),optional,check(shape(t, 0) == ldt),depend(t) :: ldt=shape(t, 0)
22+
double precision intent(in,out,copy),dimension(ldb,*) :: b
23+
integer, intent(hide),optional,check(shape(b, 0) == ldb),depend(b) :: ldb=shape(b, 0)
24+
double precision intent(cache,hide),dimension(ldwork) :: dwork
25+
integer optional,check(ldwork>=n*k*k+(n+2)*k), depend(n,k) :: ldwork=max(1,n*k*k+(n+2)*k)
26+
integer intent(out):: info
27+
end subroutine mb02ed
28+
1529
subroutine mb03rd(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
1630
character intent(in) :: jobx
1731
character intent(in),required :: sort

‎slycot/tests/test_mb.py‎

Lines changed: 171 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,181 @@
99
from numpy.testing import assert_allclose
1010
from scipy.linalg import schur
1111

12-
from slycot import math, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
13-
from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning
12+
from slycot import math, mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
13+
from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning, SlycotParameterError
1414

1515
from .test_exceptions import assert_docstring_parse
1616

1717

18+
def test_mb02ed_ex():
19+
"""Test MB02ED using the example given in the MB02ED SLICOT Documentation"""
20+
n = 3
21+
k = 3
22+
nrhs = 2
23+
TYPET = "C"
24+
T = np.array(
25+
[
26+
[3.0000, 1.0000, 0.2000],
27+
[1.0000, 4.0000, 0.4000],
28+
[0.2000, 0.4000, 5.0000],
29+
[0.1000, 0.1000, 0.2000],
30+
[0.2000, 0.0400, 0.0300],
31+
[0.0500, 0.2000, 0.1000],
32+
[0.1000, 0.0300, 0.1000],
33+
[0.0400, 0.0200, 0.2000],
34+
[0.0100, 0.0300, 0.0200],
35+
]
36+
)
37+
B = np.array(
38+
[
39+
[1.0000, 2.0000],
40+
[1.0000, 2.0000],
41+
[1.0000, 2.0000],
42+
[1.0000, 2.0000],
43+
[1.0000, 2.0000],
44+
[1.0000, 2.0000],
45+
[1.0000, 2.0000],
46+
[1.0000, 2.0000],
47+
[1.0000, 2.0000],
48+
]
49+
)
50+
X = np.array(
51+
[
52+
[0.2408, 0.4816],
53+
[0.1558, 0.3116],
54+
[0.1534, 0.3068],
55+
[0.2302, 0.4603],
56+
[0.1467, 0.2934],
57+
[0.1537, 0.3075],
58+
[0.2349, 0.4698],
59+
[0.1498, 0.2995],
60+
[0.1653, 0.3307],
61+
]
62+
)
63+
64+
result,_ = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs)
65+
np.testing.assert_almost_equal(result, X, decimal=4)
66+
67+
def test_mb02ed_parameter_errors():
68+
"""Test for errors in the input parameters of MB02ED"""
69+
n = 3
70+
k = 3
71+
nrhs = 2
72+
TYPET = "C"
73+
T = np.array(
74+
[
75+
[3.0000, 1.0000, 0.2000],
76+
[1.0000, 4.0000, 0.4000],
77+
[0.2000, 0.4000, 5.0000],
78+
[0.1000, 0.1000, 0.2000],
79+
[0.2000, 0.0400, 0.0300],
80+
[0.0500, 0.2000, 0.1000],
81+
[0.1000, 0.0300, 0.1000],
82+
[0.0400, 0.0200, 0.2000],
83+
[0.0100, 0.0300, 0.0200],
84+
]
85+
)
86+
B = np.array(
87+
[
88+
[1.0000, 2.0000],
89+
[1.0000, 2.0000],
90+
[1.0000, 2.0000],
91+
[1.0000, 2.0000],
92+
[1.0000, 2.0000],
93+
[1.0000, 2.0000],
94+
[1.0000, 2.0000],
95+
[1.0000, 2.0000],
96+
[1.0000, 2.0000],
97+
]
98+
)
99+
X = np.array(
100+
[
101+
[0.2408, 0.4816],
102+
[0.1558, 0.3116],
103+
[0.1534, 0.3068],
104+
[0.2302, 0.4603],
105+
[0.1467, 0.2934],
106+
[0.1537, 0.3075],
107+
[0.2349, 0.4698],
108+
[0.1498, 0.2995],
109+
[0.1653, 0.3307],
110+
]
111+
)
112+
113+
# Test for wrong parameter typet
114+
with pytest.raises(expected_exception=SlycotParameterError, match='typet must be either "R" or "C"') as cm:
115+
mb02ed(T=T, B=B, n=n, k=k, typet='U', nrhs=nrhs)
116+
assert cm.value.info == -1
117+
#Test for negative number of columns
118+
with pytest.raises(expected_exception=SlycotParameterError, match="k must be >= 0") as cm:
119+
mb02ed(T=T, B=B, n=n, k=-1, typet=TYPET, nrhs=nrhs)
120+
assert cm.value.info == -2
121+
#Test for negative number of blocks
122+
with pytest.raises(expected_exception=SlycotParameterError, match="n must be >= 0") as cm:
123+
mb02ed(T=T, B=B, n=-1, k=k, typet=TYPET, nrhs=nrhs)
124+
assert cm.value.info == -3
125+
#Test for negative number of right hand sides
126+
with pytest.raises(expected_exception=SlycotParameterError, match="nrhs must be >= 0") as cm:
127+
mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=-1)
128+
assert cm.value.info == -4
129+
130+
131+
def test_mb02ed_matrix_error():
132+
"""Test for a negative definite input matrix in MB02ED"""
133+
n = 3
134+
k = 3
135+
nrhs = 2
136+
TYPET = "C"
137+
T = np.array(
138+
[
139+
[3.0000, 1.0000, 0.2000],
140+
[1.0000, 4.0000, 0.4000],
141+
[0.2000, 0.4000, 5.0000],
142+
[0.1000, 0.1000, 0.2000],
143+
[0.2000, 0.0400, 0.0300],
144+
[0.0500, 0.2000, 0.1000],
145+
[0.1000, 0.0300, 0.1000],
146+
[0.0400, 0.0200, 0.2000],
147+
[0.0100, 0.0300, 0.0200],
148+
]
149+
)
150+
# Create a negative definite matrix
151+
T = -1 * T
152+
B = np.array(
153+
[
154+
[1.0000, 2.0000],
155+
[1.0000, 2.0000],
156+
[1.0000, 2.0000],
157+
[1.0000, 2.0000],
158+
[1.0000, 2.0000],
159+
[1.0000, 2.0000],
160+
[1.0000, 2.0000],
161+
[1.0000, 2.0000],
162+
[1.0000, 2.0000],
163+
]
164+
)
165+
X = np.array(
166+
[
167+
[0.2408, 0.4816],
168+
[0.1558, 0.3116],
169+
[0.1534, 0.3068],
170+
[0.2302, 0.4603],
171+
[0.1467, 0.2934],
172+
[0.1537, 0.3075],
173+
[0.2349, 0.4698],
174+
[0.1498, 0.2995],
175+
[0.1653, 0.3307],
176+
]
177+
)
178+
179+
with pytest.raises(SlycotArithmeticError,
180+
match = "The reduction algorithm failed. "
181+
"The Toeplitz matrix associated\nwith T "
182+
r"is not numerically positive definite.") as cm:
183+
mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs)
184+
assert cm.value.info == 1
185+
186+
18187
def test_mb03rd():
19188
""" Test for Schur form reduction.
20189

0 commit comments

Comments
 (0)