From 048592db253037d760cc390f5b9416fd2283f244 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 19 May 2020 02:29:46 +0200 Subject: [PATCH 1/6] add test for ab01nd --- slycot/tests/CMakeLists.txt | 1 + slycot/tests/test_ab01.py | 66 +++++++++++++++++++++++++++++++++++++ slycot/tests/test_ab08n.py | 4 +-- 3 files changed, 69 insertions(+), 2 deletions(-) create mode 100755 slycot/tests/test_ab01.py diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 7a5d126e..2aafab3f 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -1,6 +1,7 @@ set(PYSOURCE __init__.py + test_ab01.py test_ab08n.py test_ag08bd.py test_examples.py diff --git a/slycot/tests/test_ab01.py b/slycot/tests/test_ab01.py new file mode 100755 index 00000000..ec544b13 --- /dev/null +++ b/slycot/tests/test_ab01.py @@ -0,0 +1,66 @@ +""" +Test ab01 wrappers + +@author: bnavigator +""" + +from numpy import array +from numpy.testing import assert_allclose, assert_equal + +from scipy.linalg.lapack import dorgqr + +from slycot.analysis import ab01nd + + +def test_ab01nd(): + """SLICOT doc example + + http://slicot.org/objects/software/shared/doc/AB01ND.html""" + + # Example program data + n = 3 + m = 2 + tol = 0.0 + + A = array([[-1., 0., 0.], + [-2., -2., -2.], + [-1., 0., -3.]]) + B = array([[1., 0.], + [0, 2.], + [0., 1.]]) + + for jobz in ['N', 'I', 'F']: + Ac, Bc, ncont, indcon, nblk, Z, tau = ab01nd(n, m, A, B, + jobz=jobz, tol=tol) + + # The transformed state dynamics matrix of a controllable realization + assert_allclose(Ac[:ncont, :ncont], array([[-3.0000, 2.2361], + [ 0.0000, -1.0000]]), + atol=0.0001) + + # and the dimensions of its diagonal blocks are + assert_equal(nblk[:indcon], array([2])) + + # The transformed input/state matrix B of a controllable realization + assert_allclose(Bc[:ncont, :],array([[ 0.0000, -2.2361], + [ 1.0000, 0.0000]]), + atol=0.0001) + + # The controllability index of the transformed system representation + assert indcon == 1 + + if jobz == 'N': + assert Z is None + continue + elif jobz == 'I': + Z_ = Z + elif jobz == 'F': + Z_, _, info = dorgqr(Z, tau) + assert info == 0 + + # The similarity transformation matrix Z + assert_allclose(Z_, array([[ 0.0000, 1.0000, 0.0000], + [-0.8944, 0.0000, -0.4472], + [-0.4472, 0.0000, 0.8944]]), + atol=0.0001) + diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py index bdeb0b88..81d1000f 100644 --- a/slycot/tests/test_ab08n.py +++ b/slycot/tests/test_ab08n.py @@ -1,5 +1,5 @@ # =================================================== -# ag08bd tests +# ab08nX tests import unittest from slycot import analysis @@ -9,7 +9,7 @@ from numpy.testing import assert_equal, assert_allclose -class test_ab08n(unittest.TestCase): +class test_ab08nX(unittest.TestCase): """ Test regular pencil construction ab08nX with input parameters according to example in documentation """ From 3708ea19fe7fa6d4bb830093852f86dea91c2fac Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 19 May 2020 02:31:21 +0200 Subject: [PATCH 2/6] fix ab01nd: no tuple assignment. +numpydoc --- slycot/analysis.py | 149 ++++++++++++++++++++++++--------------------- 1 file changed, 79 insertions(+), 70 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 4ae19b27..f38fe793 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -21,19 +21,23 @@ from .exceptions import raise_if_slycot_error, SlycotParameterError -def ab01nd(n,m,A,B,jobz='N',tol=0,ldwork=None): +def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None): """ Ac,Bc,ncont,indcon,nblk,Z,tau = ab01nd_i(n,m,A,B,[jobz,tol,ldwork]) To find a controllable realization for the linear time-invariant multi-input system + :: + dX/dt = A * X + B * U, where A and B are N-by-N and N-by-M matrices, respectively, which are reduced by this routine to orthogonal canonical form using (and optionally accumulating) orthogonal similarity - transformations. Specifically, the pair (A, B) is reduced to - the pair (Ac, Bc), Ac = Z' * A * Z, Bc = Z' * B, given by + transformations. Specifically, the pair ``(A, B)`` is reduced to + the pair ``(Ac, Bc), Ac = Z' * A * Z, Bc = Z' * B``, given by + + :: [ Acont * ] [ Bcont ] Ac = [ ], Bc = [ ], @@ -49,69 +53,73 @@ def ab01nd(n,m,A,B,jobz='N',tol=0,ldwork=None): [ . . . . . ] [ . ] [ 0 0 . . . Ap,p-1 App ] [ 0 ] - where the blocks B1, A21, ..., Ap,p-1 have full row ranks and - p is the controllability index of the pair. The size of the - block Auncont is equal to the dimension of the uncontrollable - subspace of the pair (A, B). + where the blocks ``B1, A21, ..., Ap,p-1`` have full row ranks and + `p` is the controllability index of the pair. The size of the + block `Auncont` is equal to the dimension of the uncontrollable + subspace of the pair ``(A, B)``. - Required arguments: - n : input int - The order of the original state-space representation, i.e. - the order of the matrix A. N > 0. - m : input int - The number of system inputs, or of columns of B. M > 0. - A : rank-2 array('d') with bounds (n,n) - The leading n-by-n part of this array must contain the original - state dynamics matrix A. - B : rank-2 array('d') with bounds (n,m) - The leading n-by-m part of this array must contain the input - matrix B. - Optional arguments: - jobz := 'N' input string(len=1) - Indicates whether the user wishes to accumulate in a matrix Z - the orthogonal similarity transformations for reducing the system, - as follows: - = 'N': Do not form Z and do not store the orthogonal transformations; - = 'F': Do not form Z, but store the orthogonal transformations in - the factored form; - = 'I': Z is initialized to the unit matrix and the orthogonal - transformation matrix Z is returned. - tol := 0 input float - The tolerance to be used in rank determination when transforming - (A, B). If tol <= 0 a default value is used. - ldwork := max(n,3*m) input int - The length of the cache array. ldwork >= max(n, 3*m). - For optimum performance it should be larger. - Return objects: - Ac : rank-2 array('d') with bounds (n,n) - The leading ncont-by-ncont part contains the upper block - Hessenberg state dynamics matrix Acont in Ac, given by Z'*A*Z, - of a controllable realization for the original system. The - elements below the first block-subdiagonal are set to zero. - Bc : rank-2 array('d') with bounds (n,m) - The leading ncont-by-m part of this array contains the transformed - input matrix Bcont in Bc, given by Z'*B, with all elements but the - first block set to zero. - ncont : int - The order of the controllable state-space representation. - indcon : int - The controllability index of the controllable part of the system - representation. - nblk : rank-1 array('i') with bounds (n) - The leading indcon elements of this array contain the the orders of - the diagonal blocks of Acont. - Z : rank-2 array('d') with bounds (n,n) - If jobz = 'I', then the leading N-by-N part of this array contains - the matrix of accumulated orthogonal similarity transformations - which reduces the given system to orthogonal canonical form. - If jobz = 'F', the elements below the diagonal, with the array tau, - represent the orthogonal transformation matrix as a product of - elementary reflectors. The transformation matrix can then be - obtained by calling the LAPACK Library routine DORGQR. - If jobz = 'N', the array Z is None. - tau : rank-1 array('d') with bounds (n) - The elements of tau contain the scalar factors of the - elementary reflectors used in the reduction of B and A.""" + Parameters + ---------- + n : int + The order of the original state-space representation, i.e. + the order of the matrix A. ``n > 0``. + m : int + The number of system inputs, or of columns of B. ``m > 0``. + A : (n,n) array_like + The original state dynamics matrix A. + B : (n,m) array_like + The input matrix B. + jobz : {'N', 'F', 'I'}, optional + Indicates whether the user wishes to accumulate in a matrix Z + the orthogonal similarity transformations for reducing the system, + as follows: + + := 'N': Do not form Z and do not store the orthogonal transformations; + (default) + := 'F': Do not form Z, but store the orthogonal transformations in + the factored form; + := 'I': Z is initialized to the unit matrix and the orthogonal + transformation matrix Z is returned. + tol : float, optional + The tolerance to be used in rank determination when transforming + ``(A, B)``. If ``tol <= 0`` a default value is used. + ldwork : int, optional + The length of the cache array. ``ldwork >= max(n, 3*m)``. + For optimum performance it should be larger. + default: ``ldwork = max(n, 3*m)`` + + Returns + ------- + Ac : (n,n) ndarray + The leading ncont-by-ncont part contains the upper block + Hessenberg state dynamics matrix Acont in Ac, given by Z'*A*Z, + of a controllable realization for the original system. The + elements below the first block-subdiagonal are set to zero. + Bc : (n,m) ndarray + The leading ncont-by-m part of this array contains the transformed + input matrix Bcont in Bc, given by ``Z'*B``, with all elements but the + first block set to zero. + ncont : int + The order of the controllable state-space representation. + indcon : int + The controllability index of the controllable part of the system + representation. + nblk : (n,) int ndarray + The leading indcon elements of this array contain the the orders of + the diagonal blocks of Acont. + Z : (n,n) ndarray + - If jobz = 'I', then the leading N-by-N part of this array contains + the matrix of accumulated orthogonal similarity transformations + which reduces the given system to orthogonal canonical form. + - If jobz = 'F', the elements below the diagonal, with the array tau, + represent the orthogonal transformation matrix as a product of + elementary reflectors. The transformation matrix can then be + obtained by calling the LAPACK Library routine DORGQR. + - If jobz = 'N', the array Z is `None`. + tau : (n, ) ndarray + The elements of tau contain the scalar factors of the + elementary reflectors used in the reduction of B and A. + """ hidden = ' (hidden by the wrapper)' arg_list = ['jobz', 'n', 'm', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, @@ -125,12 +133,13 @@ def ab01nd(n,m,A,B,jobz='N',tol=0,ldwork=None): if ldwork is None: ldwork = max(n, 3*m) - out = wrappermap[jobz](n, m, A, B, tol=tol, ldwork=ldwork) - raise_if_slycot_error(out[-1], arg_list) - # sets Z to None + Ac, Bc, ncont, indcon, nblk, Z, tau, info = \ + wrappermap[jobz](n, m, A, B, tol=tol, ldwork=ldwork) + raise_if_slycot_error(info, arg_list) + if jobz == "N": - out[5] = None - return out[:-1] + Z = None + return Ac, Bc, ncont, indcon, nblk, Z, tau def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'): From 8a99d919ac8c311108aef81b8908ef5fe8ed13c9 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 19 May 2020 02:43:08 +0200 Subject: [PATCH 3/6] merge the wrapper suite for mb01nd --- slycot/analysis.py | 8 ++----- slycot/src/analysis.pyf | 49 +++-------------------------------------- 2 files changed, 5 insertions(+), 52 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index f38fe793..6073c013 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -126,15 +126,11 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None): 'ncont', 'indcon', 'nblk', 'Z', 'LDZ'+hidden, 'tau', 'tol', 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden] - wrappermap = {"N": _wrapper.ab01nd_n, - "I": _wrapper.ab01nd_i, - "F": _wrapper.ab01nd_f} - if ldwork is None: ldwork = max(n, 3*m) - Ac, Bc, ncont, indcon, nblk, Z, tau, info = \ - wrappermap[jobz](n, m, A, B, tol=tol, ldwork=ldwork) + Ac, Bc, ncont, indcon, nblk, Z, tau, info = _wrapper.ab01nd( + jobz, n, m, A, B, tol=tol, ldwork=ldwork) raise_if_slycot_error(info, arg_list) if jobz == "N": diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 2904949c..01172cfe 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -1,7 +1,6 @@ ! -*- f90 -*- -subroutine ab01nd_i(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f - fortranname ab01nd - character intent(hide) :: jobz = 'I' +subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f + character :: jobz integer check(n>0) :: n integer check(n>0) :: m double precision intent(in,out),dimension(n,n),depend(n) :: a @@ -19,49 +18,7 @@ subroutine ab01nd_i(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,d double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork integer :: ldwork = max(n,3*m) integer intent(out) :: info -end subroutine ab01nd_i -subroutine ab01nd_f(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f - fortranname ab01nd - character intent(hide) :: jobz = 'F' - integer check(n>0) :: n - integer check(n>0) :: m - double precision intent(in,out),dimension(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=shape(a,0) - double precision intent(in,out),dimension(n,m),depend(n,m) :: b - integer intent(hide),depend(b) :: ldb=shape(b,0) - integer intent(out) :: ncont - integer intent(out) :: indcon - integer intent(out),dimension(n),depend(n) :: nblk - double precision intent(out),dimension(n,n),depend(n) :: z - integer intent(hide),depend(z) :: ldz=shape(z,0) - double precision intent(out),dimension(n),depend(n) :: tau - double precision :: tol = 0 - integer intent(hide,cache),dimension(m),depend(m) :: iwork - double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork - integer :: ldwork = max(n,3*m) - integer intent(out) :: info -end subroutine ab01nd_f -subroutine ab01nd_n(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f - fortranname ab01nd - character intent(hide) :: jobz = 'N' - integer check(n>0) :: n - integer check(n>0) :: m - double precision intent(in,out),dimension(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=shape(a,0) - double precision intent(in,out),dimension(n,m),depend(n,m) :: b - integer intent(hide),depend(b) :: ldb=shape(b,0) - integer intent(out) :: ncont - integer intent(out) :: indcon - integer intent(out),dimension(n),depend(n) :: nblk - double precision intent(out),dimension(1,1),depend(n) :: z - integer intent(hide),depend(z) :: ldz=shape(z,0) - double precision intent(out),dimension(n),depend(n) :: tau - double precision :: tol = 0 - integer intent(hide,cache),dimension(m),depend(m) :: iwork - double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork - integer :: ldwork = max(n,3*m) - integer intent(out) :: info -end subroutine ab01nd_n +end subroutine ab01nd 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 0fd455a4e3b2f78887ff0e5b1acb862c708c8347 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 19 May 2020 03:18:20 +0200 Subject: [PATCH 4/6] allocate less memory fo Z if unneeded --- slycot/src/analysis.pyf | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 01172cfe..3115d983 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -10,8 +10,8 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo integer intent(out) :: ncont integer intent(out) :: indcon integer intent(out),dimension(n),depend(n) :: nblk - double precision intent(out),dimension(n,n),depend(n) :: z - integer intent(hide),depend(z) :: ldz=shape(z,0) + double precision intent(out),dimension(ldz,n),depend(n,ldz) :: z + integer intent(hide),depend(n) :: ldz = (*jobz == 'N' ? 1 : n) double precision intent(out),dimension(n),depend(n) :: tau double precision :: tol = 0 integer intent(hide,cache),dimension(m),depend(m) :: iwork From 822f69d52e1d2d61f6fd5d924928f4985a78cf26 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 19 May 2020 05:31:33 +0200 Subject: [PATCH 5/6] fix typo in wrapper --- slycot/src/analysis.pyf | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 3115d983..3b36aeb3 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -2,7 +2,7 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwork,ldwork,info) ! in AB01ND.f character :: jobz integer check(n>0) :: n - integer check(n>0) :: m + integer check(m>0) :: m double precision intent(in,out),dimension(n,n),depend(n) :: a integer intent(hide),depend(a) :: lda=shape(a,0) double precision intent(in,out),dimension(n,m),depend(n,m) :: b @@ -11,7 +11,7 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo integer intent(out) :: indcon integer intent(out),dimension(n),depend(n) :: nblk double precision intent(out),dimension(ldz,n),depend(n,ldz) :: z - integer intent(hide),depend(n) :: ldz = (*jobz == 'N' ? 1 : n) + integer intent(hide),depend(n, jobz) :: ldz = (*jobz == 'N' ? 1 : n) double precision intent(out),dimension(n),depend(n) :: tau double precision :: tol = 0 integer intent(hide,cache),dimension(m),depend(m) :: iwork From 9862672ee3a97d221641bc8e8074ff857ab94ff0 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Tue, 19 May 2020 12:51:20 +0200 Subject: [PATCH 6/6] update header comment in test_ab08n AB08NX is a Fortran subroutine --- slycot/tests/test_ab08n.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py index 81d1000f..06fd0c0d 100644 --- a/slycot/tests/test_ab08n.py +++ b/slycot/tests/test_ab08n.py @@ -1,5 +1,5 @@ # =================================================== -# ab08nX tests +# ab08n* tests import unittest from slycot import analysis