Skip to content

Implement MB02ED #214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bnavigator merged 13 commits into from
Sep 24, 2023
Merged

Implement MB02ED #214

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
84c9552
mb02ed
saasaa Sep 19, 2023
05fd0b7
fix formatting
saasaa Sep 20, 2023
d373407
added numpydoc style docstring
saasaa Sep 21, 2023
15188ef
fixed typos in docstring
saasaa Sep 21, 2023
bad3c74
Added new test, added checking of input variables
saasaa Sep 23, 2023
7c14647
Merge branch 'python-control:master' into master
saasaa Sep 23, 2023
3a97769
remove checks in wrapper
saasaa Sep 23, 2023
da77ca9
changed Warns to Raises in docstring
saasaa Sep 23, 2023
7da2c4c
remove unused import
saasaa Sep 23, 2023
56609fc
Update docstring, fix checks in signature,
saasaa Sep 23, 2023
bf0d133
remove checks in signature
saasaa Sep 23, 2023
eb93764
Apply suggestions from code review
saasaa Sep 23, 2023
05797af
update docstring
saasaa Sep 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions slycot/__init__.py
View file
Open in desktop
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
104 changes: 103 additions & 1 deletion slycot/math.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,113 @@
# MA 02110-1301, USA.

from . import _wrapper
from .exceptions import raise_if_slycot_error
from .exceptions import SlycotParameterError, raise_if_slycot_error
Comment thread
bnavigator marked this conversation as resolved.
Outdated

import numpy as np


def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int):
Comment thread
saasaa marked this conversation as resolved.
Outdated
""" 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'.
Comment thread
saasaa marked this conversation as resolved.
T : ndarray
The leading k-by-n*k or n*k-by-k part of this array must contain the first
block row/column of an s.p.d. block Toeplitz matrix.
B : ndarray
The leading nrhs-by-n*k or n*k-by-nrhs part of this array must contain the
right-hand side matrix B.
Comment thread
saasaa marked this conversation as resolved.
Outdated
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
Comment thread
saasaa marked this conversation as resolved.
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).

Warns
Comment thread
bnavigator marked this conversation as resolved.
Outdated
-----
SlycotArithmeticError
:info = 1:
the reduction algorithm failed. The Toeplitz matrix associated
Comment thread
saasaa marked this conversation as resolved.
Outdated
with T is not numerically positive definite.
SlycotParameterError
:info < 0:
if INFO = -i, the i-th argument had an illegal value.

@bnavigator bnavigator Sep 23, 2023

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

>       with pytest.raises(expected_exception=SlycotParameterError, match="k must be >= 0") as cm:
E       AssertionError: Regex pattern did not match.
E        Regex: 'k must be >= 0'
E        Input: '\nif INFO = -i, the i-th argument had an illegal value.'

If you want to have more meaningful error messages (not mandatory);

Suggested change
SlycotParameterError
:info < 0:
if INFO = -i, the i-th argument had an illegal value.
SlycotParameterError
:info = -2:
k must be >= 0
:info = -3:
n must be >= 0

and so on

@bnavigator bnavigator Sep 23, 2023

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Otherwise just remove the info < 0 zero case and make sure the arg_list is in the correct order, then raise_if_slycot_error will generate a reasonable message. You will have to adjust the regex pattern in the tests then.


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",
]

if k < 0:
raise SlycotParameterError(message="k must be >= 0",info=-2)
if n< 0:
raise SlycotParameterError(message="n must be >= 0", info=-3)
if nrhs < 0:
raise SlycotParameterError(message="nrhs must be >= 0", info=-4)
Comment thread
bnavigator marked this conversation as resolved.
Outdated

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
Open in desktop
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) :: k
integer intent(in),required,check(n>=0) :: n
integer intent(in),required,check(n>=0) :: 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
172 changes: 170 additions & 2 deletions slycot/tests/test_mb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,180 @@
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 negative number of columns
with pytest.raises(expected_exception=SlycotParameterError, match="k must be >= 0") as cm:
result,_ = 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:
result,_ = mb02ed(T=T, B=B, n=-1, k=n, 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:
result,_ = mb02ed(T=T, B=B, n=n, k=n, typet=TYPET, nrhs=-1)
assert cm.value.info == -4



def test_mb02ed_matrix_error():
Comment thread
saasaa marked this conversation as resolved.
"""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 = r"the reduction algorithm failed."
r"The Toeplitz matrix associated with 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():
Comment thread
saasaa marked this conversation as resolved.
""" Test for Schur form reduction.

Expand Down