From 84c9552c30a73e4a889d78fc379b67864a4667e9 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Tue, 19 Sep 2023 20:27:03 +0200 Subject: [PATCH 01/12] mb02ed --- slycot/__init__.py | 90 +++++---- slycot/math.py | 410 ++++++++++++++++++++++++---------------- slycot/src/math.pyf | 14 ++ slycot/tests/test_mb.py | 136 +++++++++---- 4 files changed, 422 insertions(+), 228 deletions(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index 9d8704a8..4c3b52b8 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -6,10 +6,10 @@ if __SLYCOT_SETUP__: import sys as _sys - _sys.stderr.write('Running from Slycot source directory.\n') + + _sys.stderr.write("Running from Slycot source directory.\n") del _sys else: - # import slycot.examples # The Slycot library is organised by 11-chapters. Each chapter can be identified by a single letter. @@ -26,16 +26,27 @@ # T : Transformation Routines (included) # U : Utility Routines - # Analysis routines (17/60 wrapped) - from .analysis import (ab01nd, - ab04md, - ab05md, ab05nd, - ab07nd, - ab08nd, ab08nz, - ab09ad, ab09ax, ab09bd, ab09md, ab09nd, - ab13bd, ab13dd, ab13ed, ab13fd, ab13md) - + from .analysis import ( + ab01nd, + ab04md, + ab05md, + ab05nd, + ab07nd, + ab08nd, + ab08nz, + ab09ad, + ab09ax, + ab09bd, + ab09md, + ab09nd, + ab13bd, + ab13dd, + ab13ed, + ab13fd, + ab13md, + ) + # Benchmark routines (0/6 wrapped) # Adaptive control routines (0/0 wrapped) @@ -46,37 +57,52 @@ # Identification routines (0/15 wrapped) - # Mathematical routines (7/281 wrapped) - from .math import (mb03rd, mb03vd, mb03vy, mb03wd, - mb05md, mb05nd, - mc01td) + # 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, - sb03md, sb03md57, sb03od, - sb04md, sb04qd, - sb10ad, sb10dd, sb10fd, sb10hd, sb10yd, - sg02ad, - sg03ad, sg03bd) - + from .synthesis import ( + sb01bd, + 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, - tb04ad, - tb05ad, - tc01od, tc04ad, - td04ad, - tf01md, tf01rd) + from .transform import ( + tb01id, + tb01pd, + tb03ad, + tb04ad, + tb05ad, + tc01od, + tc04ad, + td04ad, + tf01md, + tf01rd, + ) # Utility routines (0/7 wrapped) - from .version import __version__ def test(): import pytest - pytest.main(['--pyargs', 'slycot']) + + pytest.main(["--pyargs", "slycot"]) diff --git a/slycot/math.py b/slycot/math.py index ae9aab18..f7121926 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -17,13 +17,35 @@ # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, # MA 02110-1301, USA. +import numpy as np from . import _wrapper from .exceptions import raise_if_slycot_error -import numpy as np + +def mb02ed(n: int, k: int, T: np.ndarray, B: np.ndarray, typet: str, nrhs: int): + hidden = " (hidden by the wrapper)" + arg_list = [ + "typet", + "k", + "n", + "nrhs", + "t", + "ldt" + hidden, + "b", + "ldb" + hidden, + "ldwork" + hidden, + "dwork" + hidden, + "info", + ] + + t, b, info = _wrapper.mb02ed(typet="C", k=k, n=n, nrhs=nrhs, t=T, b=B) + + raise_if_slycot_error(info, arg_list) + + return b -def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0): +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]) To reduce a matrix `A` in real Schur form to a block-diagonal form @@ -223,30 +245,44 @@ def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0): SIAM J. Numer. Anal., 20, pp. 599-610, 1983. """ - hidden = ' (hidden by the wrapper)' - arg_list = ['jobx', 'sort', 'n', 'pmax', - 'A', 'lda' + hidden, 'X', 'ldx' + hidden, - 'nblcks', 'blsize', 'wr', 'wi', 'tol', - 'dwork' + hidden, 'info'] + hidden = " (hidden by the wrapper)" + arg_list = [ + "jobx", + "sort", + "n", + "pmax", + "A", + "lda" + hidden, + "X", + "ldx" + hidden, + "nblcks", + "blsize", + "wr", + "wi", + "tol", + "dwork" + hidden, + "info", + ] if X is None: X = np.eye(n) Ar, Xr, nblcks, blsize, wr, wi, info = _wrapper.mb03rd( - jobx, sort, n, pmax, A, X, tol) + jobx, sort, n, pmax, A, X, tol + ) raise_if_slycot_error(info, arg_list) - if jobx == 'N': + if jobx == "N": Xr = None else: Xr = Xr[:n, :n] Ar = Ar[:n, :n] - W = wr + 1J*wi + W = wr + 1j * wi return Ar, Xr, blsize[:nblcks], W def mb03vd(n, ilo, ihi, A): - """ HQ, Tau = mb03vd(n, ilo, ihi, A) + """HQ, Tau = mb03vd(n, ilo, ihi, A) To reduce a product of p real general matrices A = A_1*A_2*...*A_p to upper Hessenberg form, H = H_1*H_2*...*H_p, where H_1 is @@ -354,11 +390,20 @@ def mb03vd(n, ilo, ihi, A): (1,2) in A_p is also unchanged for this example.) """ - hidden = ' (hidden by the wrapper)' - arg_list = ['n', 'p' + hidden, - 'ilo', 'ihi', 'a', - 'lda1' + hidden, 'lda2' + hidden, 'tau', - 'ldtau' + hidden, 'dwork' + hidden, 'info'] + hidden = " (hidden by the wrapper)" + arg_list = [ + "n", + "p" + hidden, + "ilo", + "ihi", + "a", + "lda1" + hidden, + "lda2" + hidden, + "tau", + "ldtau" + hidden, + "dwork" + hidden, + "info", + ] HQ, Tau, info = _wrapper.mb03vd(n, ilo, ihi, A) @@ -367,7 +412,7 @@ def mb03vd(n, ilo, ihi, A): def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): - """ Q = mb03vy(n, ilo, ihi, A, Tau, [ldwork]) + """Q = mb03vy(n, ilo, ihi, A, Tau, [ldwork]) To generate the real orthogonal matrices Q_1, Q_2, ..., Q_p, which are defined as the product of ihi-ilo elementary reflectors @@ -411,11 +456,21 @@ def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): """ - hidden = ' (hidden by the wrapper)' - arg_list = ['n', 'p' + hidden, - 'ilo', 'ihi', 'a', - 'lda1' + hidden, 'lda2' + hidden, 'tau', - 'ldtau' + hidden, 'dwork' + hidden, 'ldwork', 'info' + hidden] + hidden = " (hidden by the wrapper)" + arg_list = [ + "n", + "p" + hidden, + "ilo", + "ihi", + "a", + "lda1" + hidden, + "lda2" + hidden, + "tau", + "ldtau" + hidden, + "dwork" + hidden, + "ldwork", + "info" + hidden, + ] if not ldwork: ldwork = max(1, 2 * n) @@ -428,139 +483,155 @@ def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): def mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork=None): - """ T, Z, Wr = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, [ldwork]) - - To compute the Schur decomposition and the eigenvalues of a - product of matrices, H = H_1*H_2*...*H_p, with H_1 an upper - Hessenberg matrix and H_2, ..., H_p upper triangular matrices, - without evaluating the product. Specifically, the matrices Z_i - are computed, such that - - :: - - Z_1' * H_1 * Z_2 = T_1, - Z_2' * H_2 * Z_3 = T_2, - ... - Z_p' * H_p * Z_1 = T_p, - - where T_1 is in real Schur form, and T_2, ..., T_p are upper - triangular. - - The routine works primarily with the Hessenberg and triangular - submatrices in rows and columns ilo to ihi, but optionally applies - the transformations to all the rows and columns of the matrices - H_i, i = 1,...,p. The transformations can be optionally - accumulated. - - Parameters - ---------- - job : {'E', 'S'} - Indicates whether the user wishes to compute the full - Schur form or the eigenvalues only, as follows: - = 'E': Compute the eigenvalues only; - = 'S': Compute the factors T_1, ..., T_p of the full - Schur form, T = T_1*T_2*...*T_p. - compz : {'N', 'I', 'V'} - Indicates whether or not the user wishes to accumulate - the matrices Z_1, ..., Z_p, as follows: - = 'N': The matrices Z_1, ..., Z_p are not required; - = 'I': Z_i is initialized to the unit matrix and the - orthogonal transformation matrix Z_i is returned, - i = 1, ..., p; - = 'V': Z_i must contain an orthogonal matrix Q_i on - entry, and the product Q_i*Z_i is returned, - i = 1, ..., p. - n : int - The order of the matrix H. n >= 0 - ilo, ihi : int - It is assumed that all matrices H_j, j = 2, ..., p, are - already upper triangular in rows and columns [:ilo-1] and - [ihi:n], and H_1 is upper quasi-triangular in rows and - columns [:ilo-1] and [ihi:n], with H_1[ilo-1,ilo-2] = 0 - (unless ilo = 1), and H_1[ihi,ihi-1] = 0 (unless ihi = n). - The routine works primarily with the Hessenberg submatrix - in rows and columns ilo to ihi, but applies the - transformations to all the rows and columns of the - matrices H_i, i = 1,...,p, if JOB = 'S'. - 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. - iloz, ihiz : int - Specify the rows of Z to which the transformations must be - applied if compz = 'I' or compz = 'V'. - 1 <= iloz <= ilo; ihi <= ihiz <= n. - H : array_like - H[:n,:n,0] must contain the upper Hessenberg matrix H_1 and - H[:n,:n,j-1] for j > 1 must contain the upper triangular matrix - H_j, j = 2, ..., p. - Q : array_like - If compz = 'V', Q[:n,:n,:p] must contain the current matrix Q of - transformations accumulated by SLICOT Library routine - MB03VY. - If compz = 'I', Q is ignored - ldwork : int, optional - The length of the cache array. The default value is - ihi-ilo+p-1 - - Returns - ------- - T : ndarray - 3D array with the same shape as H. - If JOB = 'S', T[:n,:n,0] is upper quasi-triangular in rows - and columns [ilo-1:ihi], with any 2-by-2 diagonal blocks - corresponding to a pair of complex conjugated eigenvalues, and - T[:n,:n,j-1] for j > 1 contains the resulting upper - triangular matrix T_j. - If job = 'E', T is None - Z : ndarray - 3D array with the same shape as Q. - If compz = 'V', or compz = 'I', the leading - N-by-N-by-P part of this array contains the transformation - matrices which produced the Schur form; the - transformations are applied only to the submatrices - Z[iloz-1:ihiz,ilo-1:ihi,j-1], j = 1, ..., p. - If compz = 'N', Z is None - W : (n,) complex ndarray - The computed eigenvalues ilo to ihi. If two eigenvalues - are computed as a complex conjugate pair, they are stored - in consecutive elements of W say the i-th and - (i+1)th, with imag(W][i]) > 0 and imag(W[i+1]) < 0. - If JOB = 'S', the eigenvalues are stored in the same order - as on the diagonal of the Schur form returned in H. - - Warns - ----- - SlycotResultWarning - :info > 0: - failed to compute all the eigenvalues {ilo} to {ihi} - in a total of 30*({ihi}-{ilo}+1) iterations - the elements Wr[{info}:{ihi}] contains those - eigenvalues which have been successfully computed. + """T, Z, Wr = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, [ldwork]) + + To compute the Schur decomposition and the eigenvalues of a + product of matrices, H = H_1*H_2*...*H_p, with H_1 an upper + Hessenberg matrix and H_2, ..., H_p upper triangular matrices, + without evaluating the product. Specifically, the matrices Z_i + are computed, such that + + :: + + Z_1' * H_1 * Z_2 = T_1, + Z_2' * H_2 * Z_3 = T_2, + ... + Z_p' * H_p * Z_1 = T_p, + + where T_1 is in real Schur form, and T_2, ..., T_p are upper + triangular. + + The routine works primarily with the Hessenberg and triangular + submatrices in rows and columns ilo to ihi, but optionally applies + the transformations to all the rows and columns of the matrices + H_i, i = 1,...,p. The transformations can be optionally + accumulated. + + Parameters + ---------- + job : {'E', 'S'} + Indicates whether the user wishes to compute the full + Schur form or the eigenvalues only, as follows: + = 'E': Compute the eigenvalues only; + = 'S': Compute the factors T_1, ..., T_p of the full + Schur form, T = T_1*T_2*...*T_p. + compz : {'N', 'I', 'V'} + Indicates whether or not the user wishes to accumulate + the matrices Z_1, ..., Z_p, as follows: + = 'N': The matrices Z_1, ..., Z_p are not required; + = 'I': Z_i is initialized to the unit matrix and the + orthogonal transformation matrix Z_i is returned, + i = 1, ..., p; + = 'V': Z_i must contain an orthogonal matrix Q_i on + entry, and the product Q_i*Z_i is returned, + i = 1, ..., p. + n : int + The order of the matrix H. n >= 0 + ilo, ihi : int + It is assumed that all matrices H_j, j = 2, ..., p, are + already upper triangular in rows and columns [:ilo-1] and + [ihi:n], and H_1 is upper quasi-triangular in rows and + columns [:ilo-1] and [ihi:n], with H_1[ilo-1,ilo-2] = 0 + (unless ilo = 1), and H_1[ihi,ihi-1] = 0 (unless ihi = n). + The routine works primarily with the Hessenberg submatrix + in rows and columns ilo to ihi, but applies the + transformations to all the rows and columns of the + matrices H_i, i = 1,...,p, if JOB = 'S'. + 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. + iloz, ihiz : int + Specify the rows of Z to which the transformations must be + applied if compz = 'I' or compz = 'V'. + 1 <= iloz <= ilo; ihi <= ihiz <= n. + H : array_like + H[:n,:n,0] must contain the upper Hessenberg matrix H_1 and + H[:n,:n,j-1] for j > 1 must contain the upper triangular matrix + H_j, j = 2, ..., p. + Q : array_like + If compz = 'V', Q[:n,:n,:p] must contain the current matrix Q of + transformations accumulated by SLICOT Library routine + MB03VY. + If compz = 'I', Q is ignored + ldwork : int, optional + The length of the cache array. The default value is + ihi-ilo+p-1 + + Returns + ------- + T : ndarray + 3D array with the same shape as H. + If JOB = 'S', T[:n,:n,0] is upper quasi-triangular in rows + and columns [ilo-1:ihi], with any 2-by-2 diagonal blocks + corresponding to a pair of complex conjugated eigenvalues, and + T[:n,:n,j-1] for j > 1 contains the resulting upper + triangular matrix T_j. + If job = 'E', T is None + Z : ndarray + 3D array with the same shape as Q. + If compz = 'V', or compz = 'I', the leading + N-by-N-by-P part of this array contains the transformation + matrices which produced the Schur form; the + transformations are applied only to the submatrices + Z[iloz-1:ihiz,ilo-1:ihi,j-1], j = 1, ..., p. + If compz = 'N', Z is None + W : (n,) complex ndarray + The computed eigenvalues ilo to ihi. If two eigenvalues + are computed as a complex conjugate pair, they are stored + in consecutive elements of W say the i-th and + (i+1)th, with imag(W][i]) > 0 and imag(W[i+1]) < 0. + If JOB = 'S', the eigenvalues are stored in the same order + as on the diagonal of the Schur form returned in H. + + Warns + ----- + SlycotResultWarning + :info > 0: + failed to compute all the eigenvalues {ilo} to {ihi} + in a total of 30*({ihi}-{ilo}+1) iterations + the elements Wr[{info}:{ihi}] contains those + eigenvalues which have been successfully computed. """ - hidden = ' (hidden by the wrapper)' - arg_list = ['job', 'compz', 'n', 'p' + hidden, - 'ilo', 'ihi', 'iloz', 'ihiz', - 'h', 'ldh1' + hidden, 'ldh2' + hidden, - 'z', 'ldz1' + hidden, 'ldz2' + hidden, - 'wr', 'wi', - 'dwork' + hidden, 'ldwork', 'info' + hidden] + hidden = " (hidden by the wrapper)" + arg_list = [ + "job", + "compz", + "n", + "p" + hidden, + "ilo", + "ihi", + "iloz", + "ihiz", + "h", + "ldh1" + hidden, + "ldh2" + hidden, + "z", + "ldz1" + hidden, + "ldz2" + hidden, + "wr", + "wi", + "dwork" + hidden, + "ldwork", + "info" + hidden, + ] if not ldwork: p = H.shape[2] ldwork = max(1, ihi - ilo + p - 1) T, Z, Wr, Wi, info = _wrapper.mb03wd( - job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork) + job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork + ) raise_if_slycot_error(info, arg_list, mb03rd.__doc__, locals()) - if job == 'E': + if job == "E": T = None - if compz == 'N': + if compz == "N": Z = None - W = Wr + Wi*1J + W = Wr + Wi * 1j return (T, Z, W) -def mb05md(a, delta, balanc='N'): +def mb05md(a, delta, balanc="N"): """Ar, Vr, Yr, w = mb05md(a, delta, balanc='N') Matrix exponential for a real non-defective matrix @@ -637,26 +708,37 @@ def mb05md(a, delta, balanc='N'): :info == n+2: Matrix A is defective, possibly due to rounding errors. """ - hidden = ' (hidden by the wrapper)' - arg_list = ['balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden, - 'y', 'ldy'+hidden, 'valr', 'vali', - 'iwork'+hidden, 'dwork'+hidden, 'ldwork'+hidden, - 'info'+hidden] + hidden = " (hidden by the wrapper)" + arg_list = [ + "balanc", + "n", + "delta", + "a", + "lda" + hidden, + "v", + "ldv" + hidden, + "y", + "ldy" + hidden, + "valr", + "vali", + "iwork" + hidden, + "dwork" + hidden, + "ldwork" + hidden, + "info" + hidden, + ] n = min(a.shape) - (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md(balanc=balanc, - n=n, - delta=delta, - a=a) + (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md( + balanc=balanc, n=n, delta=delta, a=a + ) raise_if_slycot_error(INFO, arg_list, mb05md.__doc__, locals()) if not all(VALi == 0): - w = VALr + 1J*VALi + w = VALr + 1j * VALi else: w = VALr return (Ar, Vr, Yr, w) - def mb05nd(a, delta, tol=1e-7): """F, H = mb05nd(n, a, delta, tol=1e-7) @@ -699,17 +781,27 @@ def mb05nd(a, delta, tol=1e-7): representable number near the overflow threshold of the machine (see LAPACK Library Routine DLAMCH). """ - hidden = ' (hidden by the wrapper)' - arg_list = ['n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden, - 'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden, - 'dwork'+hidden, 'ldwork'+hidden] + hidden = " (hidden by the wrapper)" + arg_list = [ + "n", + "delta", + "a", + "lda" + hidden, + "ex", + "ldex" + hidden, + "exint", + "ldexin" + hidden, + "tol", + "iwork" + hidden, + "dwork" + hidden, + "ldwork" + hidden, + ] n = min(a.shape) out = _wrapper.mb05nd(n=n, delta=delta, a=a, tol=tol) raise_if_slycot_error(out[-1], arg_list, mb05nd.__doc__, locals()) return out[:-1] - def mc01td(dico, dp, p): """dp, stable, nz = mc01td(dico, dp, p) @@ -768,9 +860,9 @@ def mc01td(dico, dp, p): ``P(DB+1-j) = 0.0`` on entry for ``j = 0, 1,..., k-1`` and ``P(DB+1-k) <> 0.0``. """ - hidden = ' (hidden by the wrapper)' - arg_list = ['dico', 'dp', 'P', 'stable', 'nz', 'DWORK' + hidden, - 'IWARN', 'INFO'] + hidden = " (hidden by the wrapper)" + arg_list = ["dico", "dp", "P", "stable", + "nz", "DWORK" + hidden, "IWARN", "INFO"] (dp_out, stable_log, nz, iwarn, info) = _wrapper.mc01td(dico, dp, p) raise_if_slycot_error([iwarn, info], arg_list, mc01td.__doc__, locals()) ftrue, ffalse = _wrapper.ftruefalse() diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 7dc53ff9..254987a0 100644 --- 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) :: k + integer intent(in),required,check(n>=0) :: n + integer intent(in),required,check(n>=0) :: nrhs + double precision intent(in,out),dimension(ldt,*) :: t + integer, intent(hide),optional,check(shape(t, 0) == ldt),depend(t) :: ldt=shape(t, 0) + double precision intent(in,out),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 index cdb84433..adb2907e 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -15,49 +15,111 @@ from .test_exceptions import assert_docstring_parse +def test_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 = [ + [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 = math.mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) + print(result) + np.testing.assert_almost_equal(result, X, decimal=4) + + def test_mb03rd(): - """ Test for Schur form reduction. + """Test for Schur form reduction. RvP, 31 Jul 2019""" - test1_A = np.array([ - [ 1., -1., 1., 2., 3., 1., 2., 3.], - [ 1., 1., 3., 4., 2., 3., 4., 2.], - [ 0., 0., 1., -1., 1., 5., 4., 1.], - [ 0., 0., 0., 1., -1., 3., 1., 2.], - [ 0., 0., 0., 1., 1., 2., 3., -1.], - [ 0., 0., 0., 0., 0., 1., 5., 1.], - [ 0., 0., 0., 0., 0., 0., 0.99999999, -0.99999999 ], - [ 0., 0., 0., 0., 0., 0., 0.99999999, 0.99999999 ] - ]) + test1_A = np.array( + [ + [1.0, -1.0, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0], + [1.0, 1.0, 3.0, 4.0, 2.0, 3.0, 4.0, 2.0], + [0.0, 0.0, 1.0, -1.0, 1.0, 5.0, 4.0, 1.0], + [0.0, 0.0, 0.0, 1.0, -1.0, 3.0, 1.0, 2.0], + [0.0, 0.0, 0.0, 1.0, 1.0, 2.0, 3.0, -1.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 5.0, 1.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99999999, -0.99999999], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99999999, 0.99999999], + ] + ) test1_n = test1_A.shape[0] - test1_Ar = np.array([ - [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000 ], - [ 1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ], - ]) - - test1_Xr = np.array([ - [ 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957 ], - [ 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755 ], - [ 0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148 ], - [ 0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534 ], - [ 0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000 ] - ]) - - test1_W = np.array([1+1j, 1-1j, - 1+1j, 1-1j, - 0.99999+0.99999j, 0.99999-0.99999j, - 1., 1.]) + test1_Ar = np.array( + [ + [1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000], + [1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000], + [0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000], + [0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000], + [0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000], + [0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000], + [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850], + [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000], + ] + ) + + test1_Xr = np.array( + [ + [1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957], + [0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755], + [0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148], + [0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534], + [0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801], + [0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267], + [0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000], + [0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000], + ] + ) + + test1_W = np.array( + [ + 1 + 1j, + 1 - 1j, + 1 + 1j, + 1 - 1j, + 0.99999 + 0.99999j, + 0.99999 - 0.99999j, + 1.0, + 1.0, + ] + ) test1_pmax = 1e3 test1_tol = 0.01 From 05fd0b75963b543864f56e73355ca038701dca6d Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Wed, 20 Sep 2023 12:09:27 +0200 Subject: [PATCH 02/12] fix formatting --- slycot/__init__.py | 84 +++------ slycot/math.py | 389 ++++++++++++++++------------------------ slycot/tests/test_mb.py | 93 ++++------ 3 files changed, 227 insertions(+), 339 deletions(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index 4c3b52b8..af0ed13a 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -6,10 +6,10 @@ if __SLYCOT_SETUP__: import sys as _sys - - _sys.stderr.write("Running from Slycot source directory.\n") + _sys.stderr.write('Running from Slycot source directory.\n') del _sys else: + # import slycot.examples # The Slycot library is organised by 11-chapters. Each chapter can be identified by a single letter. @@ -26,26 +26,15 @@ # T : Transformation Routines (included) # U : Utility Routines + # Analysis routines (17/60 wrapped) - from .analysis import ( - ab01nd, - ab04md, - ab05md, - ab05nd, - ab07nd, - ab08nd, - ab08nz, - ab09ad, - ab09ax, - ab09bd, - ab09md, - ab09nd, - ab13bd, - ab13dd, - ab13ed, - ab13fd, - ab13md, - ) + from .analysis import (ab01nd, + ab04md, + ab05md, ab05nd, + ab07nd, + ab08nd, ab08nz, + ab09ad, ab09ax, ab09bd, ab09md, ab09nd, + ab13bd, ab13dd, ab13ed, ab13fd, ab13md) # Benchmark routines (0/6 wrapped) @@ -58,51 +47,36 @@ # Identification routines (0/15 wrapped) # Mathematical routines (8/281 wrapped) - from .math import mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd, mc01td + 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, - sb03md, - sb03md57, - sb03od, - sb04md, - sb04qd, - sb10ad, - sb10dd, - sb10fd, - sb10hd, - sb10yd, - sg02ad, - sg03ad, - sg03bd, - ) + from .synthesis import (sb01bd, + 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, - tb04ad, - tb05ad, - tc01od, - tc04ad, - td04ad, - tf01md, - tf01rd, - ) + from .transform import (tb01id, tb01pd, + tb03ad, + tb04ad, + tb05ad, + tc01od, tc04ad, + td04ad, + tf01md, tf01rd) # Utility routines (0/7 wrapped) + from .version import __version__ def test(): import pytest - - pytest.main(["--pyargs", "slycot"]) + pytest.main(['--pyargs', 'slycot']) diff --git a/slycot/math.py b/slycot/math.py index f7121926..1c026406 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -17,10 +17,10 @@ # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, # MA 02110-1301, USA. -import numpy as np from . import _wrapper from .exceptions import raise_if_slycot_error +import numpy as np def mb02ed(n: int, k: int, T: np.ndarray, B: np.ndarray, typet: str, nrhs: int): hidden = " (hidden by the wrapper)" @@ -44,8 +44,7 @@ def mb02ed(n: int, k: int, T: np.ndarray, B: np.ndarray, typet: str, nrhs: int): return b - -def mb03rd(n, A, X=None, jobx="U", sort="N", pmax=1.0, tol=0.0): +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]) To reduce a matrix `A` in real Schur form to a block-diagonal form @@ -245,44 +244,30 @@ def mb03rd(n, A, X=None, jobx="U", sort="N", pmax=1.0, tol=0.0): SIAM J. Numer. Anal., 20, pp. 599-610, 1983. """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "jobx", - "sort", - "n", - "pmax", - "A", - "lda" + hidden, - "X", - "ldx" + hidden, - "nblcks", - "blsize", - "wr", - "wi", - "tol", - "dwork" + hidden, - "info", - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['jobx', 'sort', 'n', 'pmax', + 'A', 'lda' + hidden, 'X', 'ldx' + hidden, + 'nblcks', 'blsize', 'wr', 'wi', 'tol', + 'dwork' + hidden, 'info'] if X is None: X = np.eye(n) Ar, Xr, nblcks, blsize, wr, wi, info = _wrapper.mb03rd( - jobx, sort, n, pmax, A, X, tol - ) + jobx, sort, n, pmax, A, X, tol) raise_if_slycot_error(info, arg_list) - if jobx == "N": + if jobx == 'N': Xr = None else: Xr = Xr[:n, :n] Ar = Ar[:n, :n] - W = wr + 1j * wi + W = wr + 1J*wi return Ar, Xr, blsize[:nblcks], W def mb03vd(n, ilo, ihi, A): - """HQ, Tau = mb03vd(n, ilo, ihi, A) + """ HQ, Tau = mb03vd(n, ilo, ihi, A) To reduce a product of p real general matrices A = A_1*A_2*...*A_p to upper Hessenberg form, H = H_1*H_2*...*H_p, where H_1 is @@ -390,20 +375,11 @@ def mb03vd(n, ilo, ihi, A): (1,2) in A_p is also unchanged for this example.) """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "n", - "p" + hidden, - "ilo", - "ihi", - "a", - "lda1" + hidden, - "lda2" + hidden, - "tau", - "ldtau" + hidden, - "dwork" + hidden, - "info", - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'p' + hidden, + 'ilo', 'ihi', 'a', + 'lda1' + hidden, 'lda2' + hidden, 'tau', + 'ldtau' + hidden, 'dwork' + hidden, 'info'] HQ, Tau, info = _wrapper.mb03vd(n, ilo, ihi, A) @@ -412,7 +388,7 @@ def mb03vd(n, ilo, ihi, A): def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): - """Q = mb03vy(n, ilo, ihi, A, Tau, [ldwork]) + """ Q = mb03vy(n, ilo, ihi, A, Tau, [ldwork]) To generate the real orthogonal matrices Q_1, Q_2, ..., Q_p, which are defined as the product of ihi-ilo elementary reflectors @@ -456,21 +432,11 @@ def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "n", - "p" + hidden, - "ilo", - "ihi", - "a", - "lda1" + hidden, - "lda2" + hidden, - "tau", - "ldtau" + hidden, - "dwork" + hidden, - "ldwork", - "info" + hidden, - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'p' + hidden, + 'ilo', 'ihi', 'a', + 'lda1' + hidden, 'lda2' + hidden, 'tau', + 'ldtau' + hidden, 'dwork' + hidden, 'ldwork', 'info' + hidden] if not ldwork: ldwork = max(1, 2 * n) @@ -483,155 +449,139 @@ def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): def mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork=None): - """T, Z, Wr = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, [ldwork]) - - To compute the Schur decomposition and the eigenvalues of a - product of matrices, H = H_1*H_2*...*H_p, with H_1 an upper - Hessenberg matrix and H_2, ..., H_p upper triangular matrices, - without evaluating the product. Specifically, the matrices Z_i - are computed, such that - - :: - - Z_1' * H_1 * Z_2 = T_1, - Z_2' * H_2 * Z_3 = T_2, - ... - Z_p' * H_p * Z_1 = T_p, - - where T_1 is in real Schur form, and T_2, ..., T_p are upper - triangular. - - The routine works primarily with the Hessenberg and triangular - submatrices in rows and columns ilo to ihi, but optionally applies - the transformations to all the rows and columns of the matrices - H_i, i = 1,...,p. The transformations can be optionally - accumulated. - - Parameters - ---------- - job : {'E', 'S'} - Indicates whether the user wishes to compute the full - Schur form or the eigenvalues only, as follows: - = 'E': Compute the eigenvalues only; - = 'S': Compute the factors T_1, ..., T_p of the full - Schur form, T = T_1*T_2*...*T_p. - compz : {'N', 'I', 'V'} - Indicates whether or not the user wishes to accumulate - the matrices Z_1, ..., Z_p, as follows: - = 'N': The matrices Z_1, ..., Z_p are not required; - = 'I': Z_i is initialized to the unit matrix and the - orthogonal transformation matrix Z_i is returned, - i = 1, ..., p; - = 'V': Z_i must contain an orthogonal matrix Q_i on - entry, and the product Q_i*Z_i is returned, - i = 1, ..., p. - n : int - The order of the matrix H. n >= 0 - ilo, ihi : int - It is assumed that all matrices H_j, j = 2, ..., p, are - already upper triangular in rows and columns [:ilo-1] and - [ihi:n], and H_1 is upper quasi-triangular in rows and - columns [:ilo-1] and [ihi:n], with H_1[ilo-1,ilo-2] = 0 - (unless ilo = 1), and H_1[ihi,ihi-1] = 0 (unless ihi = n). - The routine works primarily with the Hessenberg submatrix - in rows and columns ilo to ihi, but applies the - transformations to all the rows and columns of the - matrices H_i, i = 1,...,p, if JOB = 'S'. - 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. - iloz, ihiz : int - Specify the rows of Z to which the transformations must be - applied if compz = 'I' or compz = 'V'. - 1 <= iloz <= ilo; ihi <= ihiz <= n. - H : array_like - H[:n,:n,0] must contain the upper Hessenberg matrix H_1 and - H[:n,:n,j-1] for j > 1 must contain the upper triangular matrix - H_j, j = 2, ..., p. - Q : array_like - If compz = 'V', Q[:n,:n,:p] must contain the current matrix Q of - transformations accumulated by SLICOT Library routine - MB03VY. - If compz = 'I', Q is ignored - ldwork : int, optional - The length of the cache array. The default value is - ihi-ilo+p-1 - - Returns - ------- - T : ndarray - 3D array with the same shape as H. - If JOB = 'S', T[:n,:n,0] is upper quasi-triangular in rows - and columns [ilo-1:ihi], with any 2-by-2 diagonal blocks - corresponding to a pair of complex conjugated eigenvalues, and - T[:n,:n,j-1] for j > 1 contains the resulting upper - triangular matrix T_j. - If job = 'E', T is None - Z : ndarray - 3D array with the same shape as Q. - If compz = 'V', or compz = 'I', the leading - N-by-N-by-P part of this array contains the transformation - matrices which produced the Schur form; the - transformations are applied only to the submatrices - Z[iloz-1:ihiz,ilo-1:ihi,j-1], j = 1, ..., p. - If compz = 'N', Z is None - W : (n,) complex ndarray - The computed eigenvalues ilo to ihi. If two eigenvalues - are computed as a complex conjugate pair, they are stored - in consecutive elements of W say the i-th and - (i+1)th, with imag(W][i]) > 0 and imag(W[i+1]) < 0. - If JOB = 'S', the eigenvalues are stored in the same order - as on the diagonal of the Schur form returned in H. - - Warns - ----- - SlycotResultWarning - :info > 0: - failed to compute all the eigenvalues {ilo} to {ihi} - in a total of 30*({ihi}-{ilo}+1) iterations - the elements Wr[{info}:{ihi}] contains those - eigenvalues which have been successfully computed. + """ T, Z, Wr = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, [ldwork]) + + To compute the Schur decomposition and the eigenvalues of a + product of matrices, H = H_1*H_2*...*H_p, with H_1 an upper + Hessenberg matrix and H_2, ..., H_p upper triangular matrices, + without evaluating the product. Specifically, the matrices Z_i + are computed, such that + + :: + + Z_1' * H_1 * Z_2 = T_1, + Z_2' * H_2 * Z_3 = T_2, + ... + Z_p' * H_p * Z_1 = T_p, + + where T_1 is in real Schur form, and T_2, ..., T_p are upper + triangular. + + The routine works primarily with the Hessenberg and triangular + submatrices in rows and columns ilo to ihi, but optionally applies + the transformations to all the rows and columns of the matrices + H_i, i = 1,...,p. The transformations can be optionally + accumulated. + + Parameters + ---------- + job : {'E', 'S'} + Indicates whether the user wishes to compute the full + Schur form or the eigenvalues only, as follows: + = 'E': Compute the eigenvalues only; + = 'S': Compute the factors T_1, ..., T_p of the full + Schur form, T = T_1*T_2*...*T_p. + compz : {'N', 'I', 'V'} + Indicates whether or not the user wishes to accumulate + the matrices Z_1, ..., Z_p, as follows: + = 'N': The matrices Z_1, ..., Z_p are not required; + = 'I': Z_i is initialized to the unit matrix and the + orthogonal transformation matrix Z_i is returned, + i = 1, ..., p; + = 'V': Z_i must contain an orthogonal matrix Q_i on + entry, and the product Q_i*Z_i is returned, + i = 1, ..., p. + n : int + The order of the matrix H. n >= 0 + ilo, ihi : int + It is assumed that all matrices H_j, j = 2, ..., p, are + already upper triangular in rows and columns [:ilo-1] and + [ihi:n], and H_1 is upper quasi-triangular in rows and + columns [:ilo-1] and [ihi:n], with H_1[ilo-1,ilo-2] = 0 + (unless ilo = 1), and H_1[ihi,ihi-1] = 0 (unless ihi = n). + The routine works primarily with the Hessenberg submatrix + in rows and columns ilo to ihi, but applies the + transformations to all the rows and columns of the + matrices H_i, i = 1,...,p, if JOB = 'S'. + 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. + iloz, ihiz : int + Specify the rows of Z to which the transformations must be + applied if compz = 'I' or compz = 'V'. + 1 <= iloz <= ilo; ihi <= ihiz <= n. + H : array_like + H[:n,:n,0] must contain the upper Hessenberg matrix H_1 and + H[:n,:n,j-1] for j > 1 must contain the upper triangular matrix + H_j, j = 2, ..., p. + Q : array_like + If compz = 'V', Q[:n,:n,:p] must contain the current matrix Q of + transformations accumulated by SLICOT Library routine + MB03VY. + If compz = 'I', Q is ignored + ldwork : int, optional + The length of the cache array. The default value is + ihi-ilo+p-1 + + Returns + ------- + T : ndarray + 3D array with the same shape as H. + If JOB = 'S', T[:n,:n,0] is upper quasi-triangular in rows + and columns [ilo-1:ihi], with any 2-by-2 diagonal blocks + corresponding to a pair of complex conjugated eigenvalues, and + T[:n,:n,j-1] for j > 1 contains the resulting upper + triangular matrix T_j. + If job = 'E', T is None + Z : ndarray + 3D array with the same shape as Q. + If compz = 'V', or compz = 'I', the leading + N-by-N-by-P part of this array contains the transformation + matrices which produced the Schur form; the + transformations are applied only to the submatrices + Z[iloz-1:ihiz,ilo-1:ihi,j-1], j = 1, ..., p. + If compz = 'N', Z is None + W : (n,) complex ndarray + The computed eigenvalues ilo to ihi. If two eigenvalues + are computed as a complex conjugate pair, they are stored + in consecutive elements of W say the i-th and + (i+1)th, with imag(W][i]) > 0 and imag(W[i+1]) < 0. + If JOB = 'S', the eigenvalues are stored in the same order + as on the diagonal of the Schur form returned in H. + + Warns + ----- + SlycotResultWarning + :info > 0: + failed to compute all the eigenvalues {ilo} to {ihi} + in a total of 30*({ihi}-{ilo}+1) iterations + the elements Wr[{info}:{ihi}] contains those + eigenvalues which have been successfully computed. """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "job", - "compz", - "n", - "p" + hidden, - "ilo", - "ihi", - "iloz", - "ihiz", - "h", - "ldh1" + hidden, - "ldh2" + hidden, - "z", - "ldz1" + hidden, - "ldz2" + hidden, - "wr", - "wi", - "dwork" + hidden, - "ldwork", - "info" + hidden, - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['job', 'compz', 'n', 'p' + hidden, + 'ilo', 'ihi', 'iloz', 'ihiz', + 'h', 'ldh1' + hidden, 'ldh2' + hidden, + 'z', 'ldz1' + hidden, 'ldz2' + hidden, + 'wr', 'wi', + 'dwork' + hidden, 'ldwork', 'info' + hidden] if not ldwork: p = H.shape[2] ldwork = max(1, ihi - ilo + p - 1) T, Z, Wr, Wi, info = _wrapper.mb03wd( - job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork - ) + job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork) raise_if_slycot_error(info, arg_list, mb03rd.__doc__, locals()) - if job == "E": + if job == 'E': T = None - if compz == "N": + if compz == 'N': Z = None - W = Wr + Wi * 1j + W = Wr + Wi*1J return (T, Z, W) -def mb05md(a, delta, balanc="N"): +def mb05md(a, delta, balanc='N'): """Ar, Vr, Yr, w = mb05md(a, delta, balanc='N') Matrix exponential for a real non-defective matrix @@ -708,37 +658,26 @@ def mb05md(a, delta, balanc="N"): :info == n+2: Matrix A is defective, possibly due to rounding errors. """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "balanc", - "n", - "delta", - "a", - "lda" + hidden, - "v", - "ldv" + hidden, - "y", - "ldy" + hidden, - "valr", - "vali", - "iwork" + hidden, - "dwork" + hidden, - "ldwork" + hidden, - "info" + hidden, - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden, + 'y', 'ldy'+hidden, 'valr', 'vali', + 'iwork'+hidden, 'dwork'+hidden, 'ldwork'+hidden, + 'info'+hidden] n = min(a.shape) - (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md( - balanc=balanc, n=n, delta=delta, a=a - ) + (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md(balanc=balanc, + n=n, + delta=delta, + a=a) raise_if_slycot_error(INFO, arg_list, mb05md.__doc__, locals()) if not all(VALi == 0): - w = VALr + 1j * VALi + w = VALr + 1J*VALi else: w = VALr return (Ar, Vr, Yr, w) + def mb05nd(a, delta, tol=1e-7): """F, H = mb05nd(n, a, delta, tol=1e-7) @@ -781,27 +720,17 @@ def mb05nd(a, delta, tol=1e-7): representable number near the overflow threshold of the machine (see LAPACK Library Routine DLAMCH). """ - hidden = " (hidden by the wrapper)" - arg_list = [ - "n", - "delta", - "a", - "lda" + hidden, - "ex", - "ldex" + hidden, - "exint", - "ldexin" + hidden, - "tol", - "iwork" + hidden, - "dwork" + hidden, - "ldwork" + hidden, - ] + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden, + 'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden, + 'dwork'+hidden, 'ldwork'+hidden] n = min(a.shape) out = _wrapper.mb05nd(n=n, delta=delta, a=a, tol=tol) raise_if_slycot_error(out[-1], arg_list, mb05nd.__doc__, locals()) return out[:-1] + def mc01td(dico, dp, p): """dp, stable, nz = mc01td(dico, dp, p) @@ -860,9 +789,9 @@ def mc01td(dico, dp, p): ``P(DB+1-j) = 0.0`` on entry for ``j = 0, 1,..., k-1`` and ``P(DB+1-k) <> 0.0``. """ - hidden = " (hidden by the wrapper)" - arg_list = ["dico", "dp", "P", "stable", - "nz", "DWORK" + hidden, "IWARN", "INFO"] + hidden = ' (hidden by the wrapper)' + arg_list = ['dico', 'dp', 'P', 'stable', 'nz', 'DWORK' + hidden, + 'IWARN', 'INFO'] (dp_out, stable_log, nz, iwarn, info) = _wrapper.mc01td(dico, dp, p) raise_if_slycot_error([iwarn, info], arg_list, mc01td.__doc__, locals()) ftrue, ffalse = _wrapper.ftruefalse() diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index adb2907e..65bc18d5 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -9,7 +9,7 @@ from numpy.testing import assert_allclose from scipy.linalg import schur -from slycot import math, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd +from slycot import math, mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning from .test_exceptions import assert_docstring_parse @@ -58,68 +58,53 @@ def test_mb02ed(): [0.1653, 0.3307], ] - result = math.mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) + result = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) print(result) np.testing.assert_almost_equal(result, X, decimal=4) - def test_mb03rd(): - """Test for Schur form reduction. + """ Test for Schur form reduction. RvP, 31 Jul 2019""" - test1_A = np.array( - [ - [1.0, -1.0, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0], - [1.0, 1.0, 3.0, 4.0, 2.0, 3.0, 4.0, 2.0], - [0.0, 0.0, 1.0, -1.0, 1.0, 5.0, 4.0, 1.0], - [0.0, 0.0, 0.0, 1.0, -1.0, 3.0, 1.0, 2.0], - [0.0, 0.0, 0.0, 1.0, 1.0, 2.0, 3.0, -1.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 5.0, 1.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99999999, -0.99999999], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99999999, 0.99999999], - ] - ) + test1_A = np.array([ + [ 1., -1., 1., 2., 3., 1., 2., 3.], + [ 1., 1., 3., 4., 2., 3., 4., 2.], + [ 0., 0., 1., -1., 1., 5., 4., 1.], + [ 0., 0., 0., 1., -1., 3., 1., 2.], + [ 0., 0., 0., 1., 1., 2., 3., -1.], + [ 0., 0., 0., 0., 0., 1., 5., 1.], + [ 0., 0., 0., 0., 0., 0., 0.99999999, -0.99999999 ], + [ 0., 0., 0., 0., 0., 0., 0.99999999, 0.99999999 ] + ]) test1_n = test1_A.shape[0] - test1_Ar = np.array( - [ - [1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000], - [1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000], - [0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000], - [0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000], - [0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000], - [0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000], - [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850], - [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000], - ] - ) - - test1_Xr = np.array( - [ - [1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957], - [0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755], - [0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148], - [0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534], - [0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801], - [0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267], - [0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000], - [0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000], - ] - ) - - test1_W = np.array( - [ - 1 + 1j, - 1 - 1j, - 1 + 1j, - 1 - 1j, - 0.99999 + 0.99999j, - 0.99999 - 0.99999j, - 1.0, - 1.0, - ] - ) + test1_Ar = np.array([ + [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000 ], + [ 1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ], + ]) + + test1_Xr = np.array([ + [ 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957 ], + [ 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755 ], + [ 0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148 ], + [ 0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534 ], + [ 0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000 ] + ]) + + test1_W = np.array([1+1j, 1-1j, + 1+1j, 1-1j, + 0.99999+0.99999j, 0.99999-0.99999j, + 1., 1.]) test1_pmax = 1e3 test1_tol = 0.01 From d37340792703b91e4919c38aecbef82c4fba0f1c Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Thu, 21 Sep 2023 14:06:06 +0000 Subject: [PATCH 03/12] added numpydoc style docstring --- slycot/math.py | 78 +++++++++++++++++++++++++++++++++++++++-- slycot/tests/test_mb.py | 2 +- 2 files changed, 76 insertions(+), 4 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index 1c026406..2669462c 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -22,7 +22,78 @@ import numpy as np -def mb02ed(n: int, k: int, T: np.ndarray, B: np.ndarray, typet: str, nrhs: int): + +def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): + """ B = 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 : 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. + 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 + On exit, 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 + ----- + SlycotResultWarning + :info < 0: + if INFO = -i, the i-th argument had an illegal value; + :info = 1: + the reduction algorithm failed. The Toeplitz matrix + associated with T is not (numerically) positive definite. + 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", @@ -38,11 +109,12 @@ def mb02ed(n: int, k: int, T: np.ndarray, B: np.ndarray, typet: str, nrhs: int): "info", ] - t, b, info = _wrapper.mb02ed(typet="C", k=k, n=n, nrhs=nrhs, t=T, b=B) + T, X, info = _wrapper.mb02ed(typet="C", k=k, n=n, nrhs=nrhs, t=T, b=B) raise_if_slycot_error(info, arg_list) - return b + 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/tests/test_mb.py b/slycot/tests/test_mb.py index 65bc18d5..39ba7c78 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -58,7 +58,7 @@ def test_mb02ed(): [0.1653, 0.3307], ] - result = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) + result,_ = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) print(result) np.testing.assert_almost_equal(result, X, decimal=4) From 15188ef79ddab48f7100483c56f9d09d0a0417c2 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Thu, 21 Sep 2023 16:14:14 +0200 Subject: [PATCH 04/12] fixed typos in docstring --- slycot/math.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index 2669462c..c2a37b42 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -40,27 +40,27 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): Note: the notation x / y means that x corresponds to typet = 'R' and y corresponds to typet = 'C'. T : ndarray - The leading K-by-N*K or N*K-by-K part of this array must contain the first + 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 + The leading nrhs-by-n*k or 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. + The number of blocks in T. n >= 0. k : int - The number of rows/columns in T, equal to the blocksize. K >= 0. + The number of rows/columns in T, equal to the blocksize. k >= 0. nrhs : int - The number of right-hand sides. NRHS >= 0. + The number of right-hand sides. nrhs >= 0. Returns ------- X : ndarray - Leading NRHS-by-N*K / N*K-by-NRHS part of + Leading nrhs-by-n*k / n*k-by-nrhs part of this array contains the solution matrix X. T: ndarray - On exit, 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). + On exit, 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 ----- From bad3c74d1116802c52af8d862367a67c33a927a9 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 18:33:50 +0000 Subject: [PATCH 05/12] Added new test, added checking of input variables --- slycot/math.py | 29 +++++--- slycot/src/math.pyf | 4 +- slycot/tests/test_mb.py | 149 ++++++++++++++++++++++++++++++++++++---- 3 files changed, 156 insertions(+), 26 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index c2a37b42..5887653a 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -18,13 +18,13 @@ # MA 02110-1301, USA. from . import _wrapper -from .exceptions import raise_if_slycot_error +from .exceptions import SlycotParameterError, raise_if_slycot_error import numpy as np def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): - """ B = mb02ed(typet, T, B, n, k, nrhs) + """ 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. @@ -58,18 +58,20 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): Leading nrhs-by-n*k / n*k-by-nrhs part of this array contains the solution matrix X. T: ndarray - On exit, if no error is thrown and nrhs > 0, then the leading + 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 ----- - SlycotResultWarning - :info < 0: - if INFO = -i, the i-th argument had an illegal value; + SlycotArithmeticError :info = 1: - the reduction algorithm failed. The Toeplitz matrix - associated with T is not (numerically) positive definite. + the reduction algorithm failed. The Toeplitz matrix associated + with T is not numerically positive definite. + SlycotParameterError + :info < 0: + if INFO = -i, the i-th argument had an illegal value. + Notes ----- The algorithm uses Householder transformations, modified hyperbolic rotations, @@ -109,9 +111,16 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): "info", ] - T, X, info = _wrapper.mb02ed(typet="C", k=k, n=n, nrhs=nrhs, t=T, b=B) + 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) - raise_if_slycot_error(info, arg_list) + 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 diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 254987a0..61a9d675 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -17,9 +17,9 @@ subroutine mb02ed(typet,k,n,nrhs,t,ldt,b,ldb,dwork,ldwork,info) ! in MB02ED.f 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),dimension(ldt,*) :: t + 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),dimension(ldb,*) :: b + 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) diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index 39ba7c78..ef28c68a 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -10,12 +10,13 @@ from scipy.linalg import schur from slycot import math, mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd -from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning +from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning, SlycotParameterError from .test_exceptions import assert_docstring_parse -def test_mb02ed(): +def test_mb02ed_ex(): + """Test MB02ED using the example given in the MB02ED SLICOT Documentation""" n = 3 k = 3 nrhs = 2 @@ -46,22 +47,142 @@ def test_mb02ed(): [1.0000, 2.0000], ] ) - X = [ - [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], - ] + 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) - print(result) 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(): + """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(): """ Test for Schur form reduction. From 3a97769f235ccea36600cc78e714927816d745f5 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 19:16:21 +0000 Subject: [PATCH 06/12] remove checks in wrapper --- slycot/math.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index 5887653a..f4c35f86 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -111,13 +111,6 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): "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) - 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()) From da77ca9db0cfb47efaf1c8d90e2c0b65b5018795 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 19:19:09 +0000 Subject: [PATCH 07/12] changed Warns to Raises in docstring --- slycot/math.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index f4c35f86..a999b507 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -62,8 +62,8 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): 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 - ----- + Raises + ------ SlycotArithmeticError :info = 1: the reduction algorithm failed. The Toeplitz matrix associated From 7da2c4cb0b5dccc64f0da875b6ac106b5fba9e4a Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 19:36:52 +0000 Subject: [PATCH 08/12] remove unused import --- slycot/math.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/slycot/math.py b/slycot/math.py index a999b507..a5574d77 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -18,7 +18,7 @@ # MA 02110-1301, USA. from . import _wrapper -from .exceptions import SlycotParameterError, raise_if_slycot_error +from .exceptions import raise_if_slycot_error import numpy as np From 56609fce448dce81799c950f392c10ec79bbc05d Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 20:18:35 +0000 Subject: [PATCH 09/12] Update docstring, fix checks in signature, add tests --- slycot/math.py | 10 ++++++++-- slycot/src/math.pyf | 4 ++-- slycot/tests/test_mb.py | 20 ++++++++++++-------- 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index a5574d77..cde7c620 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -69,8 +69,14 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): the reduction algorithm failed. The Toeplitz matrix associated with T is not numerically positive definite. SlycotParameterError - :info < 0: - if INFO = -i, the i-th argument had an illegal value. + :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 ----- diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 61a9d675..ed6b9771 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -14,9 +14,9 @@ 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(k>=0) :: k integer intent(in),required,check(n>=0) :: n - integer intent(in),required,check(n>=0) :: nrhs + integer intent(in),required,check(nrhs>=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 diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index ef28c68a..45038d55 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -110,17 +110,21 @@ def test_mb02ed_parameter_errors(): ] ) - # Test for negative number of columns + # 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: - result,_ = mb02ed(T=T, B=B, n=n, k=-1, typet=TYPET, nrhs=nrhs) + mb02ed(T=T, B=B, n=n, k=-1, typet=TYPET, nrhs=nrhs) assert cm.value.info == -2 - # Test for negative number of blocks + #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) + 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 + #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) + mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=-1) assert cm.value.info == -4 @@ -174,8 +178,8 @@ def test_mb02ed_matrix_error(): ) with pytest.raises(SlycotArithmeticError, - match = r"the reduction algorithm failed." - r"The Toeplitz matrix associated with T" + match = "\nthe 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 From bf0d1339e273d773bea6f2d1315c1b4f6c7a95fa Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 20:24:30 +0000 Subject: [PATCH 10/12] remove checks in signature --- slycot/src/math.pyf | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index ed6b9771..fe78cdde 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -14,9 +14,9 @@ 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,check(k>=0) :: k - integer intent(in),required,check(n>=0) :: n - integer intent(in),required,check(nrhs>=0) :: nrhs + 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 From eb937648e45d0ca4a3cc0ea0d2e21438a904f1bc Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 22:51:58 +0200 Subject: [PATCH 11/12] Apply suggestions from code review Co-authored-by: Ben Greiner --- slycot/math.py | 8 ++++---- slycot/tests/test_mb.py | 5 +---- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index cde7c620..23775fc9 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -23,7 +23,7 @@ import numpy as np -def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): +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 @@ -39,10 +39,10 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): 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 : ndarray + T : array_like 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 + B : array_like The leading nrhs-by-n*k or n*k-by-nrhs part of this array must contain the right-hand side matrix B. n : int @@ -66,7 +66,7 @@ def mb02ed(typet: str,T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): ------ SlycotArithmeticError :info = 1: - the reduction algorithm failed. The Toeplitz matrix associated + The reduction algorithm failed. The Toeplitz matrix associated with T is not numerically positive definite. SlycotParameterError :info = -1: diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index 45038d55..5dc5dcd5 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -128,7 +128,6 @@ def test_mb02ed_parameter_errors(): assert cm.value.info == -4 - def test_mb02ed_matrix_error(): """Test for a negative definite input matrix in MB02ED""" n = 3 @@ -178,15 +177,13 @@ def test_mb02ed_matrix_error(): ) with pytest.raises(SlycotArithmeticError, - match = "\nthe reduction algorithm failed. " + 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. From 05797af95dd377e1a13d898c8255e00c81e244a9 Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sat, 23 Sep 2023 21:00:25 +0000 Subject: [PATCH 12/12] update docstring --- slycot/math.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/slycot/math.py b/slycot/math.py index 23775fc9..d7447047 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -40,10 +40,10 @@ def mb02ed(typet: str, T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int) 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 or n*k-by-k part of this array must contain the first + 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 or n*k-by-nrhs part of this array must contain the + 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.