From 9485d94aaae2221dd82ee0bd7c66e28d610e390f Mon Sep 17 00:00:00 2001 From: Alexander Bodenseher Date: Sun, 24 Sep 2023 16:06:32 +0200 Subject: [PATCH] Implement MB02ED (#214) --- slycot/__init__.py | 10 +-- slycot/math.py | 101 +++++++++++++++++++++++ slycot/src/math.pyf | 14 ++++ slycot/tests/test_mb.py | 173 +++++++++++++++++++++++++++++++++++++++- 4 files changed, 291 insertions(+), 7 deletions(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index 9d8704a8..af0ed13a 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -35,7 +35,7 @@ ab08nd, ab08nz, ab09ad, ab09ax, ab09bd, ab09md, ab09nd, ab13bd, ab13dd, ab13ed, ab13fd, ab13md) - + # Benchmark routines (0/6 wrapped) # Adaptive control routines (0/0 wrapped) @@ -46,8 +46,8 @@ # 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) @@ -55,13 +55,13 @@ # 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, diff --git a/slycot/math.py b/slycot/math.py index ae9aab18..d7447047 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -23,6 +23,107 @@ import numpy as np +def mb02ed(typet: str, T: np.ndarray, B: np.ndarray, n: int, k: int, nrhs: int): + """ X, T = mb02ed(typet, T, B, n, k, nrhs) + + Solve a system of linear equations T*X = B or X*T = B with a positive + definite block Toeplitz matrix T. + + Parameters + ---------- + typet: str + Specifies the type of T: + - 'R': T contains the first block row of an s.p.d. block Toeplitz matrix, + and the system X*T = B is solved. + - 'C': T contains the first block column of an s.p.d. block Toeplitz matrix, + and the system T*X = B is solved. + Note: the notation x / y means that x corresponds to + typet = 'R' and y corresponds to typet = 'C'. + T : array_like + The leading k-by-n*k / n*k-by-k part of this array must contain the first + block row/column of an s.p.d. block Toeplitz matrix. + B : array_like + The leading nrhs-by-n*k / n*k-by-nrhs part of this array must contain the + right-hand side matrix B. + n : int + The number of blocks in T. n >= 0. + k : int + The number of rows/columns in T, equal to the blocksize. k >= 0. + nrhs : int + The number of right-hand sides. nrhs >= 0. + + Returns + ------- + X : ndarray + Leading nrhs-by-n*k / n*k-by-nrhs part of + this array contains the solution matrix X. + T: ndarray + If no error is thrown and nrhs > 0, then the leading + k-by-n*k / n*k-by-k part of this array contains the last + row / column of the Cholesky factor of inv(T). + + Raises + ------ + SlycotArithmeticError + :info = 1: + The reduction algorithm failed. The Toeplitz matrix associated + with T is not numerically positive definite. + SlycotParameterError + :info = -1: + typet must be either "R" or "C" + :info = -2: + k must be >= 0 + :info = -3: + n must be >= 0 + :info = -4: + nrhs must be >= 0 + + Notes + ----- + The algorithm uses Householder transformations, modified hyperbolic rotations, + and block Gaussian eliminations in the Schur algorithm [1], [2]. + + References + ---------- + [1] Kailath, T. and Sayed, A. + Fast Reliable Algorithms for Matrices with Structure. + SIAM Publications, Philadelphia, 1999. + + [2] Kressner, D. and Van Dooren, P. + Factorizations and linear system solvers for matrices with Toeplitz structure. + SLICOT Working Note 2000-2, 2000. + + Numerical Aspects + ----------------- + The implemented method is numerically equivalent to forming the Cholesky factor R and the + inverse Cholesky factor of T using the generalized Schur algorithm and solving the systems + of equations R*X = L*B or X*R = B*L by a blocked backward substitution algorithm. + The algorithm requires O(K * N^2 + K * N * NRHS) floating-point operations. + + """ + + hidden = " (hidden by the wrapper)" + arg_list = [ + "typet", + "k", + "n", + "nrhs", + "t", + "ldt" + hidden, + "b", + "ldb" + hidden, + "ldwork" + hidden, + "dwork" + hidden, + "info", + ] + + T, X, info = _wrapper.mb02ed(typet=typet, k=k, n=n, nrhs=nrhs, t=T, b=B) + + raise_if_slycot_error(info, arg_list, docstring=mb02ed.__doc__, checkvars=locals()) + + return X, T + + def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0): """Ar, Xr, blsize, W = mb03rd(n, A, [X, jobx, sort, pmax, tol]) diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 7dc53ff9..fe78cdde 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),required :: k + integer intent(in),required :: n + integer intent(in),required :: nrhs + double precision intent(in,out,copy),dimension(ldt,*) :: t + integer, intent(hide),optional,check(shape(t, 0) == ldt),depend(t) :: ldt=shape(t, 0) + double precision intent(in,out,copy),dimension(ldb,*) :: b + integer, intent(hide),optional,check(shape(b, 0) == ldb),depend(b) :: ldb=shape(b, 0) + double precision intent(cache,hide),dimension(ldwork) :: dwork + integer optional,check(ldwork>=n*k*k+(n+2)*k), depend(n,k) :: ldwork=max(1,n*k*k+(n+2)*k) + integer intent(out):: info +end subroutine mb02ed + subroutine mb03rd(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f character intent(in) :: jobx character intent(in),required :: sort diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index cdb84433..5dc5dcd5 100644 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -9,12 +9,181 @@ from numpy.testing import assert_allclose from scipy.linalg import schur -from slycot import math, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd -from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning +from slycot import math, mb02ed, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd +from slycot.exceptions import SlycotArithmeticError, SlycotResultWarning, SlycotParameterError from .test_exceptions import assert_docstring_parse +def test_mb02ed_ex(): + """Test MB02ED using the example given in the MB02ED SLICOT Documentation""" + n = 3 + k = 3 + nrhs = 2 + TYPET = "C" + T = np.array( + [ + [3.0000, 1.0000, 0.2000], + [1.0000, 4.0000, 0.4000], + [0.2000, 0.4000, 5.0000], + [0.1000, 0.1000, 0.2000], + [0.2000, 0.0400, 0.0300], + [0.0500, 0.2000, 0.1000], + [0.1000, 0.0300, 0.1000], + [0.0400, 0.0200, 0.2000], + [0.0100, 0.0300, 0.0200], + ] + ) + B = np.array( + [ + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + ] + ) + X = np.array( + [ + [0.2408, 0.4816], + [0.1558, 0.3116], + [0.1534, 0.3068], + [0.2302, 0.4603], + [0.1467, 0.2934], + [0.1537, 0.3075], + [0.2349, 0.4698], + [0.1498, 0.2995], + [0.1653, 0.3307], + ] + ) + + result,_ = mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) + np.testing.assert_almost_equal(result, X, decimal=4) + +def test_mb02ed_parameter_errors(): + """Test for errors in the input parameters of MB02ED""" + n = 3 + k = 3 + nrhs = 2 + TYPET = "C" + T = np.array( + [ + [3.0000, 1.0000, 0.2000], + [1.0000, 4.0000, 0.4000], + [0.2000, 0.4000, 5.0000], + [0.1000, 0.1000, 0.2000], + [0.2000, 0.0400, 0.0300], + [0.0500, 0.2000, 0.1000], + [0.1000, 0.0300, 0.1000], + [0.0400, 0.0200, 0.2000], + [0.0100, 0.0300, 0.0200], + ] + ) + B = np.array( + [ + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + ] + ) + X = np.array( + [ + [0.2408, 0.4816], + [0.1558, 0.3116], + [0.1534, 0.3068], + [0.2302, 0.4603], + [0.1467, 0.2934], + [0.1537, 0.3075], + [0.2349, 0.4698], + [0.1498, 0.2995], + [0.1653, 0.3307], + ] + ) + + # Test for wrong parameter typet + with pytest.raises(expected_exception=SlycotParameterError, match='typet must be either "R" or "C"') as cm: + mb02ed(T=T, B=B, n=n, k=k, typet='U', nrhs=nrhs) + assert cm.value.info == -1 + #Test for negative number of columns + with pytest.raises(expected_exception=SlycotParameterError, match="k must be >= 0") as cm: + mb02ed(T=T, B=B, n=n, k=-1, typet=TYPET, nrhs=nrhs) + assert cm.value.info == -2 + #Test for negative number of blocks + with pytest.raises(expected_exception=SlycotParameterError, match="n must be >= 0") as cm: + mb02ed(T=T, B=B, n=-1, k=k, typet=TYPET, nrhs=nrhs) + assert cm.value.info == -3 + #Test for negative number of right hand sides + with pytest.raises(expected_exception=SlycotParameterError, match="nrhs must be >= 0") as cm: + mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=-1) + assert cm.value.info == -4 + + +def test_mb02ed_matrix_error(): + """Test for a negative definite input matrix in MB02ED""" + n = 3 + k = 3 + nrhs = 2 + TYPET = "C" + T = np.array( + [ + [3.0000, 1.0000, 0.2000], + [1.0000, 4.0000, 0.4000], + [0.2000, 0.4000, 5.0000], + [0.1000, 0.1000, 0.2000], + [0.2000, 0.0400, 0.0300], + [0.0500, 0.2000, 0.1000], + [0.1000, 0.0300, 0.1000], + [0.0400, 0.0200, 0.2000], + [0.0100, 0.0300, 0.0200], + ] + ) + # Create a negative definite matrix + T = -1 * T + B = np.array( + [ + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + [1.0000, 2.0000], + ] + ) + X = np.array( + [ + [0.2408, 0.4816], + [0.1558, 0.3116], + [0.1534, 0.3068], + [0.2302, 0.4603], + [0.1467, 0.2934], + [0.1537, 0.3075], + [0.2349, 0.4698], + [0.1498, 0.2995], + [0.1653, 0.3307], + ] + ) + + with pytest.raises(SlycotArithmeticError, + match = "The reduction algorithm failed. " + "The Toeplitz matrix associated\nwith T " + r"is not numerically positive definite.") as cm: + mb02ed(T=T, B=B, n=n, k=k, typet=TYPET, nrhs=nrhs) + assert cm.value.info == 1 + + def test_mb03rd(): """ Test for Schur form reduction.