From 3598ac9e06a60151f75321da95c79311aad78b5c Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Thu, 20 Jul 2023 22:29:00 +0200 Subject: [PATCH 1/8] Add ab04md wrapper --- slycot/__init__.py | 3 ++- slycot/analysis.py | 17 +++++++++++++++++ slycot/src/analysis.pyf | 20 ++++++++++++++++++++ 3 files changed, 39 insertions(+), 1 deletion(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index 69527783..32b4a5db 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -27,8 +27,9 @@ # U : Utility Routines - # Analysis routines (16/60 wrapped) + # Analysis routines (17/60 wrapped) from .analysis import (ab01nd, + ab04md, ab05md, ab05nd, ab07nd, ab08nd, ab08nz, diff --git a/slycot/analysis.py b/slycot/analysis.py index 32050236..d2edb67e 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -139,6 +139,23 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None): Z = None return Ac, Bc, ncont, indcon, nblk, Z, tau +def ab04md(type_bn, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): + """ a,b,c,d = ab04md() + """ + + hidden = ' (hidden by the wrapper)' + arg_list = ['type_bn', 'n', 'm', 'p', 'alpha', 'beta', + 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden, + 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden] + + if ldwork is None: + ldwork = max(n, 3*m) + + out = _wrapper.ab04md(type_bn, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork) + info=out[-1] + raise_if_slycot_error(info, arg_list) + + return out[:-1] def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'): """ n,a,b,c,d = ab05md(n1,m1,p1,n2,p2,a1,b1,c1,d1,a2,b2,c2,d2,[uplo]) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 633ff6ec..28629e38 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -19,6 +19,26 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo integer :: ldwork = max(n,3*m) integer intent(out) :: info end subroutine ab01nd +subroutine ab04md(type_bn,n,m,p,alpha,beta,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in AB04MD.f + character :: type_bn + integer check(n>=0) :: n + integer check(m>=0) :: m + integer check(p>=0) :: p + double precision intent(in) :: alpha + double precision intent(in) :: beta + double precision intent(in,out,copy), dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda = shape(a,0) + double precision intent(in,out,copy), dimension(n,m),depend(n,m) :: b + integer intent(hide),depend(b) :: ldb = shape(b,0) + double precision intent(in,out,copy), dimension(p,n),depend(n,p) :: c + integer intent(hide),depend(c) :: ldc = shape(c,0) + double precision intent(in,out,copy), dimension(p,m),depend(m,p) :: d + integer intent(hide),depend(d) :: ldd = shape(d,0) + integer intent(hide,cache),dimension(n),depend(n) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork + integer :: ldwork = max(1,n) + integer intent(out) :: info +end subroutine ab04md subroutine ab05md(uplo,over,n1,m1,p1,n2,p2,a1,lda1,b1,ldb1,c1,ldc1,d1,ldd1,a2,lda2,b2,ldb2,c2,ldc2,d2,ldd2,n,a,lda,b,ldb,c,ldc,d,ldd,dwork,ldwork,info) ! in AB05MD.f character :: uplo = 'U' character intent(hide) :: over = 'N' ! not sure how the overlap works From 52cc0f8fae101ee03dd2013bb83e43cc51527a4c Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Sat, 29 Jul 2023 18:22:47 +0200 Subject: [PATCH 2/8] Unittest added for ab04md. --- slycot/tests/CMakeLists.txt | 1 + slycot/tests/test_ab04md.py | 71 +++++++++++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+) create mode 100644 slycot/tests/test_ab04md.py diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 2aafab3f..80d751da 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -2,6 +2,7 @@ set(PYSOURCE __init__.py test_ab01.py + test_ab04md.py test_ab08n.py test_ag08bd.py test_examples.py diff --git a/slycot/tests/test_ab04md.py b/slycot/tests/test_ab04md.py new file mode 100644 index 00000000..39b156da --- /dev/null +++ b/slycot/tests/test_ab04md.py @@ -0,0 +1,71 @@ +# =================================================== +# ab08n* tests + +import unittest +from slycot import analysis +import numpy as np + +from scipy.linalg import eig +from numpy.testing import assert_equal, assert_allclose + + +class test_ab04md(unittest.TestCase): + """Test ab04md. + + Example data taken from + https://www.slicot.org/objects/software/shared/doc/AB04MD.html. + """ + + Ac = np.array([[1.0, 0.5], + [0.5, 1.0]]) + Bc = np.array([[0.0, -1.0], + [1.0, 0.0]]) + Cc = np.array([[-1.0, 0.0], + [0.0, 1.0]]) + Dc = np.array([[1.0, 0.0], + [0.0, -1.0]]) + + Ad = np.array([[-1.0, -4.0], + [-4.0, -1.0]]) + Bd = np.array([[2.8284, 0.0], + [0.0, -2.8284]]) + Cd = np.array([[0.0, 2.8284], + [-2.8284, 0.0]]) + Dd = np.array([[-1.0, 0.0], + [0.0, -3.0]]) + + def test_ab04md_cont_disc_cont(self): + """Test transformation from continuous - to discrete - to continuous time. + """ + + n, m = self.Bc.shape + p = self.Cc.shape[0] + + Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md('C',n,m,p,self.Ac,self.Bc,self.Cc,self.Dc) + + Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md('D',n,m,p,Ad_t,Bd_t,Cd_t,Dd_t) + + assert_allclose(self.Ac, Ac_t) + assert_allclose(self.Bc, Bc_t) + assert_allclose(self.Cc, Cc_t) + assert_allclose(self.Dc, Dc_t) + + def test_ab04md_disc_cont_disc(self): + """Test transformation from discrete - to continuous - to discrete time. + """ + + n, m = self.Bc.shape + p = self.Cc.shape[0] + + Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md('D',n,m,p,self.Ad,self.Bd,self.Cd,self.Dd) + + Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md('C',n,m,p,Ac_t,Bc_t,Cc_t,Dc_t) + + assert_allclose(self.Ad, Ad_t) + assert_allclose(self.Bd, Bd_t) + assert_allclose(self.Cd, Cd_t) + assert_allclose(self.Dd, Dd_t) + + +if __name__ == "__main__": + unittest.main() From 416b352ec96fe5b6221fe6c4c99d2a78479424a4 Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Sat, 29 Jul 2023 19:38:21 +0200 Subject: [PATCH 3/8] Add docstring, Change parameter name from type_bn to type_t --- slycot/analysis.py | 66 +++++++++++++++++++++++++++++++++++++---- slycot/src/analysis.pyf | 4 +-- 2 files changed, 63 insertions(+), 7 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index d2edb67e..7da8eae8 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -139,19 +139,75 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None): Z = None return Ac, Bc, ncont, indcon, nblk, Z, tau -def ab04md(type_bn, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): - """ a,b,c,d = ab04md() +def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): + """ At,Bt,Ct,Dt = ab04md(type_bn, n, m, p, A, B, C, D, [alpha, beta,ldwork]) + + Parameters + ---------- + type_t : {'D','C'} + Indicates the type of the original system and the + transformation to be performed as follows: + = 'D': discrete-time -> continuous-time; + = 'C': continuous-time -> discrete-time. + n : input int + The order of the matrix A, the number of rows of matrix B and + the number of columns of matrix C. It represents the dimension of + the state vector. n > 0. + m : input int + The number of columns of matrix B. It represents the dimension of + the input vector. m > 0. + p : input int + The number of rows of matrix C. It represents the dimension of + the output vector. p > 0. + A : input rank-2 array('d') with bounds (n,n) + The leading n-by-n part of this array must contain the system state + matrix A. + B : input rank-2 array('d') with bounds (n,m) + The leading n-by-m part of this array must contain the system input + matrix B. + C : input rank-2 array('d') with bounds (p,n) + The leading p-by-n part of this array must contain the system output + matrix C. + D : input rank-2 array('d') with bounds (p,m) + The leading p-by-m part of this array must contain the system direct + transmission matrix D. + alpha : double + Parameter specifying the bilinear transformation. + Recommended values for stable systems: alpha = 1, alpha != 0, + beta : double + Parameter specifying the bilinear transformation. + Recommended values for stable systems: beta = 1, beta != 0, + ldwork : int + The length of the cache array. + ldwork >= max(1, n) + Returns + ------- + At : output rank-2 array('d') with bounds (n,n) + The state matrix At of the transformed system. + Bt : output rank-2 array('d') with bounds (n,m) + The input matrix Bt of the transformed system. + Ct : output rank-2 array('d') with bounds (p,n) + The output matrix Ct of the transformed system. + Dt : output rank-2 array('d') with bounds (p,m) + The transmission matrix Dt of the transformed system. + Raises + ------ + SlycotArithmeticError + :info == 1: + If the matrix (ALPHA*I + A) is exactly singular + :info == 2: + If the matrix (BETA*I - A) is exactly singular. """ hidden = ' (hidden by the wrapper)' - arg_list = ['type_bn', 'n', 'm', 'p', 'alpha', 'beta', + arg_list = ['type_t', 'n', 'm', 'p', 'alpha', 'beta', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden, 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden] if ldwork is None: - ldwork = max(n, 3*m) + ldwork = max(1, n) - out = _wrapper.ab04md(type_bn, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork) + out = _wrapper.ab04md(type_t, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork) info=out[-1] raise_if_slycot_error(info, arg_list) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 28629e38..b2bb5a12 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -19,8 +19,8 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo integer :: ldwork = max(n,3*m) integer intent(out) :: info end subroutine ab01nd -subroutine ab04md(type_bn,n,m,p,alpha,beta,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in AB04MD.f - character :: type_bn +subroutine ab04md(type_t,n,m,p,alpha,beta,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in AB04MD.f + character :: type_t integer check(n>=0) :: n integer check(m>=0) :: m integer check(p>=0) :: p From 6643167d1d034f8c3a8ac8e86190400d037dc40a Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Mon, 7 Aug 2023 21:52:58 +0200 Subject: [PATCH 4/8] Update slycot/analysis.py Co-authored-by: Ben Greiner --- slycot/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 7da8eae8..41162852 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -149,7 +149,7 @@ def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): transformation to be performed as follows: = 'D': discrete-time -> continuous-time; = 'C': continuous-time -> discrete-time. - n : input int + n : int The order of the matrix A, the number of rows of matrix B and the number of columns of matrix C. It represents the dimension of the state vector. n > 0. From fe5c56aaf5b011d765ab77410543021496d4c70e Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Mon, 7 Aug 2023 22:11:15 +0200 Subject: [PATCH 5/8] Update docstring (remove input mark from int types) --- slycot/analysis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 41162852..ea1227d2 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -153,10 +153,10 @@ def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): The order of the matrix A, the number of rows of matrix B and the number of columns of matrix C. It represents the dimension of the state vector. n > 0. - m : input int + m : int The number of columns of matrix B. It represents the dimension of the input vector. m > 0. - p : input int + p : int The number of rows of matrix C. It represents the dimension of the output vector. p > 0. A : input rank-2 array('d') with bounds (n,n) From 8ad02364625a54de094901f9b5a309d8f04aeebe Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Mon, 7 Aug 2023 22:33:37 +0200 Subject: [PATCH 6/8] Improve docstring --- slycot/analysis.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index ea1227d2..4fbe7678 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -159,16 +159,16 @@ def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): p : int The number of rows of matrix C. It represents the dimension of the output vector. p > 0. - A : input rank-2 array('d') with bounds (n,n) + A : (n,n) ndarray The leading n-by-n part of this array must contain the system state matrix A. - B : input rank-2 array('d') with bounds (n,m) + B : (n,m) ndarray The leading n-by-m part of this array must contain the system input matrix B. - C : input rank-2 array('d') with bounds (p,n) + C : (p,n) ndarray The leading p-by-n part of this array must contain the system output matrix C. - D : input rank-2 array('d') with bounds (p,m) + D : (p,m) ndarray The leading p-by-m part of this array must contain the system direct transmission matrix D. alpha : double @@ -182,13 +182,13 @@ def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): ldwork >= max(1, n) Returns ------- - At : output rank-2 array('d') with bounds (n,n) + At : (n,n) ndarray The state matrix At of the transformed system. - Bt : output rank-2 array('d') with bounds (n,m) + Bt : (n,m) ndarray The input matrix Bt of the transformed system. - Ct : output rank-2 array('d') with bounds (p,n) + Ct : (p,n) ndarray The output matrix Ct of the transformed system. - Dt : output rank-2 array('d') with bounds (p,m) + Dt : (p,m) ndarray The transmission matrix Dt of the transformed system. Raises ------ From 2289d93bbf8eb579327da9e6087d7753b7888932 Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Wed, 23 Aug 2023 21:55:47 +0200 Subject: [PATCH 7/8] Update docstring, minor fixes --- slycot/analysis.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 4fbe7678..f84ff43f 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -140,7 +140,7 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None): return Ac, Bc, ncont, indcon, nblk, Z, tau def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): - """ At,Bt,Ct,Dt = ab04md(type_bn, n, m, p, A, B, C, D, [alpha, beta,ldwork]) + """ At,Bt,Ct,Dt = ab04md(type_t, n, m, p, A, B, C, D, [alpha, beta,ldwork]) Parameters ---------- @@ -159,27 +159,29 @@ def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None): p : int The number of rows of matrix C. It represents the dimension of the output vector. p > 0. - A : (n,n) ndarray + A : (n,n) array_like The leading n-by-n part of this array must contain the system state matrix A. - B : (n,m) ndarray + B : (n,m) array_like The leading n-by-m part of this array must contain the system input matrix B. - C : (p,n) ndarray + C : (p,n) array_like The leading p-by-n part of this array must contain the system output matrix C. - D : (p,m) ndarray + D : (p,m) array_like The leading p-by-m part of this array must contain the system direct transmission matrix D. - alpha : double + alpha : double, optional Parameter specifying the bilinear transformation. Recommended values for stable systems: alpha = 1, alpha != 0, - beta : double + Default is 1.0. + beta : double, optional Parameter specifying the bilinear transformation. Recommended values for stable systems: beta = 1, beta != 0, - ldwork : int + Default is 1.0. + ldwork : int, optional The length of the cache array. - ldwork >= max(1, n) + ldwork >= max(1, n), default is max(1, n) Returns ------- At : (n,n) ndarray From 118bf588f69168f6afc99d4fdd01bf40bc7241ab Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 24 Aug 2023 10:13:04 +0200 Subject: [PATCH 8/8] cleanup test file --- slycot/tests/test_ab04md.py | 47 +++++++++++++++++-------------------- 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/slycot/tests/test_ab04md.py b/slycot/tests/test_ab04md.py index 39b156da..7ce59a0e 100644 --- a/slycot/tests/test_ab04md.py +++ b/slycot/tests/test_ab04md.py @@ -1,38 +1,33 @@ -# =================================================== -# ab08n* tests - -import unittest from slycot import analysis import numpy as np -from scipy.linalg import eig -from numpy.testing import assert_equal, assert_allclose +from numpy.testing import assert_allclose -class test_ab04md(unittest.TestCase): +class test_ab04md: """Test ab04md. - + Example data taken from - https://www.slicot.org/objects/software/shared/doc/AB04MD.html. + https://www.slicot.org/objects/software/shared/doc/AB04MD.html """ Ac = np.array([[1.0, 0.5], - [0.5, 1.0]]) + [0.5, 1.0]]) Bc = np.array([[0.0, -1.0], - [1.0, 0.0]]) + [1.0, 0.0]]) Cc = np.array([[-1.0, 0.0], - [0.0, 1.0]]) + [0.0, 1.0]]) Dc = np.array([[1.0, 0.0], - [0.0, -1.0]]) + [0.0, -1.0]]) Ad = np.array([[-1.0, -4.0], - [-4.0, -1.0]]) + [-4.0, -1.0]]) Bd = np.array([[2.8284, 0.0], - [0.0, -2.8284]]) + [0.0, -2.8284]]) Cd = np.array([[0.0, 2.8284], - [-2.8284, 0.0]]) + [-2.8284, 0.0]]) Dd = np.array([[-1.0, 0.0], - [0.0, -3.0]]) + [0.0, -3.0]]) def test_ab04md_cont_disc_cont(self): """Test transformation from continuous - to discrete - to continuous time. @@ -41,15 +36,17 @@ def test_ab04md_cont_disc_cont(self): n, m = self.Bc.shape p = self.Cc.shape[0] - Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md('C',n,m,p,self.Ac,self.Bc,self.Cc,self.Dc) + Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md( + 'C', n, m, p, self.Ac, self.Bc, self.Cc, self.Dc) - Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md('D',n,m,p,Ad_t,Bd_t,Cd_t,Dd_t) + Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md( + 'D', n, m, p, Ad_t, Bd_t, Cd_t, Dd_t) assert_allclose(self.Ac, Ac_t) assert_allclose(self.Bc, Bc_t) assert_allclose(self.Cc, Cc_t) assert_allclose(self.Dc, Dc_t) - + def test_ab04md_disc_cont_disc(self): """Test transformation from discrete - to continuous - to discrete time. """ @@ -57,15 +54,13 @@ def test_ab04md_disc_cont_disc(self): n, m = self.Bc.shape p = self.Cc.shape[0] - Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md('D',n,m,p,self.Ad,self.Bd,self.Cd,self.Dd) + Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md( + 'D', n, m, p, self.Ad, self.Bd, self.Cd, self.Dd) - Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md('C',n,m,p,Ac_t,Bc_t,Cc_t,Dc_t) + Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md( + 'C', n, m, p, Ac_t, Bc_t, Cc_t, Dc_t) assert_allclose(self.Ad, Ad_t) assert_allclose(self.Bd, Bd_t) assert_allclose(self.Cd, Cd_t) assert_allclose(self.Dd, Dd_t) - - -if __name__ == "__main__": - unittest.main()