From 2f78306efc326d2a9263b2a742c5dbd02a28d11f Mon Sep 17 00:00:00 2001 From: lytex Date: Wed, 1 Apr 2020 18:58:43 +0200 Subject: [PATCH 01/10] add ab08nz function to support regular pencil for complex state-space systems --- slycot/__init__.py | 3 +- slycot/analysis.py | 88 +++++++++++++++++++++++++++++++++++++++++ slycot/src/analysis.pyf | 31 +++++++++++++++ 3 files changed, 121 insertions(+), 1 deletion(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index 11b13d12..aae4d81e 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -14,7 +14,8 @@ # Analysis routines (14/40 wrapped) from .analysis import ab01nd,ab05md,ab05nd,ab07nd,ab08nd, \ - ab09ad,ab09ax,ab09bd,ab09md,ab09nd,ab13bd,ab13dd,ab13ed,ab13fd + ab09ad,ab09ax,ab09bd,ab09md,ab09nd,ab13bd,ab13dd,ab13ed,ab13fd, \ + ab08nz # Data analysis routines (0/7 wrapped) diff --git a/slycot/analysis.py b/slycot/analysis.py index 1b496d69..edfd8bad 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -475,6 +475,94 @@ def ab08nd(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): raise e return out[:-1] +def ab08nz(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): + """ nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = ab08nz(n,m,p,A,B,C,D,[equil,tol,ldwork]) + + To construct for a linear multivariable system described by a state-space + model (A,B,C,D) a regular pencil (Af - lambda*Bf ) which has the invariant + zeros of the system as generalized eigenvalues. + The routine also computes the orders of the infinite zeros and the + right and left Kronecker indices of the system (A,B,C,D). + + Required arguments: + n : input int + The number of state variables. n >= 0. + m : input int + The number of system inputs. m >= 0. + p : input int + The number of system outputs. p >= 0. + A : input rank-2 array('d') with bounds (n,n) + The leading n-by-n part of this array must contain the state + dynamics matrix A of the system. + B : input rank-2 array('d') with bounds (n,m) + The leading n-by-m part of this array must contain the input/state + matrix B of the system. + C : input rank-2 array('d') with bounds (p,n) + The leading p-by-n part of this array must contain the state/output + matrix C of the system. + D : input rank-2 array('d') with bounds (p,m) + The leading p-by-m part of this array must contain the direct + transmission matrix D of the system. + Optional arguments: + equil := 'N' input string(len=1) + Specifies whether the user wishes to balance the compound matrix + as follows: + = 'S': Perform balancing (scaling); + = 'N': Do not perform balancing. + tol := 0.0 input float + A tolerance used in rank decisions to determine the effective rank, + which is defined as the order of the largest leading (or trailing) + triangular submatrix in the QR (or RQ) factorization with column + (or row) pivoting whose estimated condition number is less than 1/tol. + ldwork := None input int + The length of the cache array. The default value is n + 3*max(m,p), + for better performance should be larger. + Return objects: + nu : int + The number of (finite) invariant zeros. + rank : int + The normal rank of the transfer function matrix. + dinfz : int + The maximum degree of infinite elementary divisors. + nkror : int + The number of right Kronecker indices. + nkrol : int + The number of left Kronecker indices. + infz : rank-1 array('i') with bounds (n) + The leading dinfz elements of infz contain information on the + infinite elementary divisors as follows: the system has infz(i) + infinite elementary divisors of degree i, where i = 1,2,...,dinfz. + kronr : rank-1 array('i') with bounds (max(n,m)+1) + the leading nkror elements of this array contain the right kronecker + (column) indices. + kronl : rank-1 array('i') with bounds (max(n,p)+1) + the leading nkrol elements of this array contain the left kronecker + (row) indices. + Af : rank-2 array('d') with bounds (max(1,n+m),n+min(p,m)) + the leading nu-by-nu part of this array contains the coefficient + matrix Af of the reduced pencil. the remainder of the leading + (n+m)-by-(n+min(p,m)) part is used as internal workspace. + Bf : rank-2 array('d') with bounds (max(1,n+p),n+m) + The leading nu-by-nu part of this array contains the coefficient + matrix Bf of the reduced pencil. the remainder of the leading + (n+p)-by-(n+m) part is used as internal workspace. + """ + hidden = ' (hidden by the wrapper)' + arg_list = ['equil', 'n', 'm', 'p', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, + 'C', 'LDC'+hidden, 'D', 'LDD'+hidden, 'nu', 'rank', 'dinfz', 'nkror', + 'nkrol', 'infz', 'kronr', 'kronl', 'Af', 'LDAF'+hidden, 'Bf', + 'LDBF'+hidden, 'tol', 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', + 'INFO'+hidden] + if ldwork is None: + ldwork = n+3*max(m,p) #only an upper bound + out = _wrapper.ab08nz(n,m,p,A,B,C,D,equil=equil,tol=tol,ldwork=ldwork) + if out[-1] < 0: + error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] + e = ValueError(error_text) + e.info = out[-1] + raise e + return out[:-1] + def ab09ad(dico,job,equil,n,m,p,A,B,C,nr=None,tol=0,ldwork=None): """ nr,Ar,Br,Cr,hsv = ab09ad(dico,job,equil,n,m,p,A,B,C,[nr,tol,ldwork]) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 58a0656b..fc412fde 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -184,6 +184,37 @@ subroutine ab08nd(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkr integer optional :: ldwork = n + 3*max(m,p) integer intent(out) :: info end subroutine ab08nd +subroutine ab08nz(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkrol,infz,kronr,kronl,af,ldaf,bf,ldbf,tol,iwork,dwork,ldwork,info) ! in :new:AB08NZ.f + character :: equil='N' + integer check(n>=0) :: n + integer check(m>=0) :: m + integer check(p>=0) :: p + complex*16 dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=shape(a,0) + complex*16 dimension(n,m),depend(n,m) :: b + integer intent(hide),depend(b) :: ldb=shape(b,0) + complex*16 dimension(p,n),depend(n,p) :: c + integer intent(hide),depend(c) :: ldc=shape(c,0) + complex*16 dimension(p,m),depend(m,p) :: d + integer intent(hide),depend(d) :: ldd=shape(d,0) + integer intent(out) :: nu + integer intent(out) :: rank_bn + integer intent(out) :: dinfz + integer intent(out) :: nkror + integer intent(out) :: nkrol + integer intent(out),dimension(n),depend(n) :: infz + integer intent(out),dimension(max(n,m)+1),depend([n,m]) :: kronr + integer intent(out),dimension(max(n,p)+1),depend([n,p]) :: kronl + complex*16 intent(out),dimension(max(1,n+m),n+min(p,m)) :: af + integer intent(hide),depend(af) :: ldaf=shape(af,0) + complex*16 intent(out),dimension(max(1,n+p),n+m) :: bf + integer intent(hide),depend(bf) :: ldbf=shape(bf,0) + double precision :: tol=0.0 + integer intent(hide,cache),dimension(max(m,p)) :: iwork + double precision intent(hide,cache),dimension(ldwork) :: dwork + integer optional :: ldwork = n + 3*max(m,p) + integer intent(out) :: info +end subroutine ab08nz subroutine ab09ad(dico,job,equil,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,hsv,tol,iwork,dwork,ldwork,iwarn,info) !in :balred:AB09AD.f character intent(in) :: dico character intent(in) :: job From 42ee06fc5b6c216f061c57c44e79e06127334440 Mon Sep 17 00:00:00 2001 From: lytex Date: Wed, 1 Apr 2020 21:42:12 +0200 Subject: [PATCH 02/10] add ab08n* tests (ab08nd executes successfully, ab08nz SEGFAULTs) --- slycot/tests/test_ab08n.py | 58 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 slycot/tests/test_ab08n.py diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py new file mode 100644 index 00000000..1823f1cc --- /dev/null +++ b/slycot/tests/test_ab08n.py @@ -0,0 +1,58 @@ +# =================================================== +# ag08bd tests + +import unittest +from slycot import analysis +import numpy as np + +from numpy.testing import assert_raises, assert_almost_equal, assert_equal + +# test input parameters + +test_A = np.array([[1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 3, 0, 0, 0], + [0, 0, 0,-4, 0, 0], + [0, 0, 0, 0,-1, 0], + [0, 0, 0, 0, 1, 3]]) + +test_B = np.array([[0 , -1], + [-1, 0], + [ 1, -1], + [ 0, 0], + [ 0, 1], + [-1, -1]]) + +test_C = np.array([[1, 0, 0, 1, 0, 0], + [0, 1, 0, 1, 0, 1], + [0, 0, 1, 0, 0, 1]]) + +test_D = np.zeros((3, 2)) + +test_A = test_A.astype(np.complex128) +test_B = test_B.astype(np.complex128) +test_C = test_C.astype(np.complex128) +test_D = test_D.astype(np.complex128) + + +class test_ab08n(unittest.TestCase): + """ test1 to 4: Verify ag08bd with input parameters according to example in documentation """ + + def test_ab08nd(self): + #test [A-lambda*E] + #B,C,D must have correct dimensions according to l,n,m and p, but cannot have zero length in any dimenstion. Then the wrapper will complain. The length is then set to one. + + nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = analysis.ab08nd(6,2,3,test_A,test_B,test_C,test_D) + + def test_ab08nz(self): + #test [A-lambda*E] + #B,C,D must have correct dimensions according to l,n,m and p, but cannot have zero length in any dimenstion. Then the wrapper will complain. The length is then set to one. + nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = analysis.ab08nz(6,2,3,test_A,test_B,test_C,test_D) + + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() From 203dc752e40f6d8cbb2c65803e4361aa095a7103 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 03:03:15 +0200 Subject: [PATCH 03/10] clean __init__ --- slycot/__init__.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/slycot/__init__.py b/slycot/__init__.py index aae4d81e..715a6214 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -12,10 +12,11 @@ # import slycot.examples - # Analysis routines (14/40 wrapped) - from .analysis import ab01nd,ab05md,ab05nd,ab07nd,ab08nd, \ - ab09ad,ab09ax,ab09bd,ab09md,ab09nd,ab13bd,ab13dd,ab13ed,ab13fd, \ - ab08nz + # Analysis routines (15/40 wrapped) + from .analysis import ab01nd, ab05md, ab05nd, ab07nd, ab08nd, ab08nz + from .analysis import ab09ad, ab09ax, ab09bd, ab09md, ab09nd + from .analysis import ab13bd, ab13dd, ab13ed, ab13fd + # Data analysis routines (0/7 wrapped) From 8f06d66047f22afd955d6afff09006643628f54c Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 03:03:39 +0200 Subject: [PATCH 04/10] rework signatures for ab08nd and ab08nz --- slycot/src/analysis.pyf | 65 +++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index fc412fde..2904949c 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -155,17 +155,17 @@ subroutine ab07nd(n,m,a,lda,b,ldb,c,ldc,d,ldd,rcond,iwork,dwork,ldwork,info) ! i end subroutine ab07nd subroutine ab08nd(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkrol,infz,kronr,kronl,af,ldaf,bf,ldbf,tol,iwork,dwork,ldwork,info) ! in :new:AB08ND.f character :: equil='N' - integer check(n>=0) :: n - integer check(m>=0) :: m - integer check(p>=0) :: p - double precision dimension(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=shape(a,0) - double precision dimension(n,m),depend(n,m) :: b - integer intent(hide),depend(b) :: ldb=shape(b,0) - double precision dimension(p,n),depend(n,p) :: c - integer intent(hide),depend(c) :: ldc=shape(c,0) - double precision dimension(p,m),depend(m,p) :: d - integer intent(hide),depend(d) :: ldd=shape(d,0) + integer intent(in),check(n>=0),required :: n + integer intent(in),check(m>=0),required :: m + integer intent(in),check(p>=0),required :: p + double precision intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a + integer intent(hide),check(lda>=max(1,n)) :: lda=shape(a,0) + double precision intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b + integer intent(hide),check(ldb>=max(1,n)) :: ldb=shape(b,0) + double precision intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c + integer intent(hide),check(ldc>=max(1,p)) :: ldc=shape(c,0) + double precision intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d + integer intent(hide),check(ldd>=max(1,p)) :: ldd=shape(d,0) integer intent(out) :: nu integer intent(out) :: rank_bn integer intent(out) :: dinfz @@ -184,35 +184,36 @@ subroutine ab08nd(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkr integer optional :: ldwork = n + 3*max(m,p) integer intent(out) :: info end subroutine ab08nd -subroutine ab08nz(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkrol,infz,kronr,kronl,af,ldaf,bf,ldbf,tol,iwork,dwork,ldwork,info) ! in :new:AB08NZ.f +subroutine ab08nz(equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nu,rank_bn,dinfz,nkror,nkrol,infz,kronr,kronl,af,ldaf,bf,ldbf,tol,iwork,dwork,zwork,lzwork,info) ! in AB08NZ.f character :: equil='N' - integer check(n>=0) :: n - integer check(m>=0) :: m - integer check(p>=0) :: p - complex*16 dimension(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=shape(a,0) - complex*16 dimension(n,m),depend(n,m) :: b - integer intent(hide),depend(b) :: ldb=shape(b,0) - complex*16 dimension(p,n),depend(n,p) :: c - integer intent(hide),depend(c) :: ldc=shape(c,0) - complex*16 dimension(p,m),depend(m,p) :: d - integer intent(hide),depend(d) :: ldd=shape(d,0) + integer intent(in),check(n>=0),required :: n + integer intent(in),check(m>=0),required :: m + integer intent(in),check(p>=0),required :: p + complex*16 intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a + integer intent(hide),check(lda>=max(1,n)) :: lda=shape(a,0) + complex*16 intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b + integer intent(hide),check(ldb>=max(1,n)) :: ldb=shape(b,0) + complex*16 intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c + integer intent(hide),check(ldc>=max(1,p)) :: ldc=shape(c,0) + complex*16 intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d + integer intent(hide),check(ldd>=max(1,p)) :: ldd=shape(d,0) integer intent(out) :: nu integer intent(out) :: rank_bn integer intent(out) :: dinfz integer intent(out) :: nkror integer intent(out) :: nkrol - integer intent(out),dimension(n),depend(n) :: infz - integer intent(out),dimension(max(n,m)+1),depend([n,m]) :: kronr - integer intent(out),dimension(max(n,p)+1),depend([n,p]) :: kronl + integer intent(out),dimension(n) :: infz + integer intent(out),dimension(max(n,m)+1) :: kronr + integer intent(out),dimension(max(n,p)+1) :: kronl complex*16 intent(out),dimension(max(1,n+m),n+min(p,m)) :: af - integer intent(hide),depend(af) :: ldaf=shape(af,0) + integer intent(hide),check(ldaf>=max(1,n+m)) :: ldaf=shape(af,0) complex*16 intent(out),dimension(max(1,n+p),n+m) :: bf - integer intent(hide),depend(bf) :: ldbf=shape(bf,0) - double precision :: tol=0.0 - integer intent(hide,cache),dimension(max(m,p)) :: iwork - double precision intent(hide,cache),dimension(ldwork) :: dwork - integer optional :: ldwork = n + 3*max(m,p) + integer intent(hide),check(ldbf>=max(1,n+p)) :: ldbf=shape(bf,0) + double precision intent(in) :: tol = 0.0 + integer intent(hide),cache,dimension(max(m,p)) :: iwork + double precision intent(hide),cache,dimension(max(n,2*max(p,m))) :: dwork + complex*16 intent(out),cache,dimension(lzwork) :: zwork + integer intent(in),check(lzwork>=max(min(p,m) + max(3*m-1,n), max(min(p,n) + max(3*p-1,max(n+p,n+m)), min(m,n) + max(3*m-1,n+m)))) :: lzwork = max(min(p,m) + max(3*m-1,n), max(min(p,n) + max(3*p-1,max(n+p,n+m)), min(m,n) + max(3*m-1,n+m))) integer intent(out) :: info end subroutine ab08nz subroutine ab09ad(dico,job,equil,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,hsv,tol,iwork,dwork,ldwork,iwarn,info) !in :balred:AB09AD.f From e2471dc43cd1f8b0dc2c0d45a5f5e9dc5ec9820b Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 03:04:21 +0200 Subject: [PATCH 05/10] rework analysis.ab08nz() function --- slycot/analysis.py | 66 ++++++++++++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 26 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index edfd8bad..96acc378 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -475,8 +475,8 @@ def ab08nd(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): raise e return out[:-1] -def ab08nz(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): - """ nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = ab08nz(n,m,p,A,B,C,D,[equil,tol,ldwork]) +def ab08nz(n, m, p, A, B, C, D, equil='N', tol=0., lzwork=None): + """ nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = ab08nz(n,m,p,A,B,C,D,[equil,tol,lzwork]) To construct for a linear multivariable system described by a state-space model (A,B,C,D) a regular pencil (Af - lambda*Bf ) which has the invariant @@ -514,9 +514,12 @@ def ab08nz(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): which is defined as the order of the largest leading (or trailing) triangular submatrix in the QR (or RQ) factorization with column (or row) pivoting whose estimated condition number is less than 1/tol. - ldwork := None input int - The length of the cache array. The default value is n + 3*max(m,p), - for better performance should be larger. + lzwork := None input int + The length of the cache array zwork. The default value is calculated + to MAX( 1, MIN(P,M) + MAX(3*M-1,N), + MIN(P,N) + MAX(3*P-1,N+P,N+M), + MIN(M,N) + MAX(3*M-1,N+M) ) + for optimum performance should be larger. Return objects: nu : int The number of (finite) invariant zeros. @@ -546,22 +549,33 @@ def ab08nz(n,m,p,A,B,C,D,equil='N',tol=0,ldwork=None): The leading nu-by-nu part of this array contains the coefficient matrix Bf of the reduced pencil. the remainder of the leading (n+p)-by-(n+m) part is used as internal workspace. + lzwork_opt : int + The optimal value of lzwork. """ hidden = ' (hidden by the wrapper)' - arg_list = ['equil', 'n', 'm', 'p', 'A', 'LDA'+hidden, 'B', 'LDB'+hidden, - 'C', 'LDC'+hidden, 'D', 'LDD'+hidden, 'nu', 'rank', 'dinfz', 'nkror', - 'nkrol', 'infz', 'kronr', 'kronl', 'Af', 'LDAF'+hidden, 'Bf', - 'LDBF'+hidden, 'tol', 'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', - 'INFO'+hidden] - if ldwork is None: - ldwork = n+3*max(m,p) #only an upper bound - out = _wrapper.ab08nz(n,m,p,A,B,C,D,equil=equil,tol=tol,ldwork=ldwork) - if out[-1] < 0: - error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] + arg_list = ['equil', 'n', 'm', 'p', + 'a', 'lda' + hidden, 'b', 'ldb' + hidden, + 'c', 'ldc' + hidden, 'd', 'ldd' + hidden, + 'nu', 'rank', 'dinfz', 'nkror', 'nkrol', 'infz', 'kronr', + 'kronl', 'af', 'ldaf' + hidden, 'bf', 'ldbf' + hidden, + 'tol', 'iwork' + hidden, 'dwork' + hidden, 'zwork', + 'lzwork', 'info'] + if lzwork is None: + lzwork = max(min(p, m) + max(3*m-1, n), + min(p, n) + max(3*p-1, n+p, n+m), + min(m, n) + max(3*m-1, n+m)) + out = _wrapper.ab08nz(n, m, p, A, B, C, D, + equil=equil, tol=tol, lzwork=lzwork) + nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf, zwork, info \ + = out + if info < 0: + error_text = "The following argument had an illegal value: " + \ + arg_list[info-1] e = ValueError(error_text) - e.info = out[-1] + e.info = info raise e - return out[:-1] + return (nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf, + int(zwork[0])) def ab09ad(dico,job,equil,n,m,p,A,B,C,nr=None,tol=0,ldwork=None): """ nr,Ar,Br,Cr,hsv = ab09ad(dico,job,equil,n,m,p,A,B,C,[nr,tol,ldwork]) @@ -817,7 +831,7 @@ def ab09ax(dico,job,n,m,p,A,B,C,nr=None,tol=0.0,ldwork=None): def ab09bd(dico,job,equil,n,m,p,A,B,C,D,nr=None,tol1=0,tol2=0,ldwork=None): """ nr,Ar,Br,Cr,Dr,hsv = ab09bd(dico,job,equil,n,m,p,A,B,C,D,[nr,tol1,tol2,ldwork]) - + To compute a reduced order model (Ar,Br,Cr,Dr) for a stable original state-space representation (A,B,C,D) by using either the square-root or the balancing-free square-root Singular @@ -1000,7 +1014,7 @@ def ab09md(dico,job,equil,n,m,p,A,B,C,alpha=None,nr=None,tol=0,ldwork=None): The number of system outputs. p >= 0. A : input rank-2 array('d'), dimension (n,n) On entry, the leading N-by-N part of this array must - contain the state dynamics matrix A. + contain the state dynamics matrix A. B : input rank-2 array('d'), dimension (n,m) On entry, the leading N-by-M part of this array must contain the original input/state matrix B. @@ -1144,7 +1158,7 @@ def ab09md(dico,job,equil,n,m,p,A,B,C,alpha=None,nr=None,tol=0,ldwork=None): def ab09nd(dico,job,equil,n,m,p,A,B,C,D,alpha=None,nr=None,tol1=0,tol2=0,ldwork=None): """ nr,Ar,Br,Cr,Dr,ns,hsv = ab09nd(dico,job,equil,n,m,p,A,B,C,D,[alpha,nr,tol1,tol2,ldwork]) - + To compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using either the square-root or the balancing-free square-root Singular @@ -1662,23 +1676,23 @@ def ag08bd(l,n,m,p,A,E,B,C,D,equil='N',tol=0.0,ldwork=None): """ Af,Ef,nrank,niz,infz,kronr,infe,kronl = ag08bd(l,n,m,p,A,E,B,C,D,[equil,tol,ldwork]) To extract from the system pencil - + ( A-lambda*E B ) S(lambda) = ( ) ( C D ) - + a regular pencil Af-lambda*Ef which has the finite Smith zeros of S(lambda) as generalized eigenvalues. The routine also computes the orders of the infinite Smith zeros and determines the singular and infinite Kronecker structure of system pencil, i.e., the right and left Kronecker indices, and the multiplicities of infinite eigenvalues. - + Required arguments: l : input int The number of rows of matrices A, B, and E. l >= 0. n : input int - The number of columns of matrices A, E, and C. n >= 0. + The number of columns of matrices A, E, and C. n >= 0. m : input int The number of columns of matrix B. m >= 0. p : input int @@ -1746,10 +1760,10 @@ def ag08bd(l,n,m,p,A,E,B,C,D,equil='N',tol=0.0,ldwork=None): """ hidden = ' (hidden by the wrapper)' arg_list = ['equil', 'l', 'n', 'm', 'p', 'A', 'lda'+hidden, 'E', 'lde'+hidden, 'B', 'ldb'+hidden, 'C', 'ldc'+hidden, 'D', 'ldd'+hidden, 'nfz', 'nrank', 'niz', 'dinfz', 'nkror', 'ninfe', 'nkrol', 'infz', 'kronr', 'infe', 'kronl', 'tol', 'iwork'+hidden, 'dwork'+hidden, 'ldwork', 'info'] - + if equil != 'S' and equil != 'N': raise ValueError('Parameter equil had an illegal value') - + if ldwork is None: ldw = max(l+p,m+n)*(m+n) + max(1,5*max(l+p,m+n)) if equil == 'S': From da58b7d2f098a47605dba4695f1a4cd09a862c48 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 03:04:43 +0200 Subject: [PATCH 06/10] extend tests for ab08nX --- slycot/tests/test_ab08n.py | 102 ++++++++++++++++++++++--------------- 1 file changed, 62 insertions(+), 40 deletions(-) diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py index 1823f1cc..a5773e4e 100644 --- a/slycot/tests/test_ab08n.py +++ b/slycot/tests/test_ab08n.py @@ -5,53 +5,75 @@ from slycot import analysis import numpy as np -from numpy.testing import assert_raises, assert_almost_equal, assert_equal +from scipy.linalg import eig +from numpy.testing import assert_equal, assert_allclose -# test input parameters -test_A = np.array([[1, 0, 0, 0, 0, 0], - [0, 1, 0, 0, 0, 0], - [0, 0, 3, 0, 0, 0], - [0, 0, 0,-4, 0, 0], - [0, 0, 0, 0,-1, 0], - [0, 0, 0, 0, 1, 3]]) - -test_B = np.array([[0 , -1], - [-1, 0], - [ 1, -1], - [ 0, 0], - [ 0, 1], - [-1, -1]]) - -test_C = np.array([[1, 0, 0, 1, 0, 0], - [0, 1, 0, 1, 0, 1], - [0, 0, 1, 0, 0, 1]]) - -test_D = np.zeros((3, 2)) - -test_A = test_A.astype(np.complex128) -test_B = test_B.astype(np.complex128) -test_C = test_C.astype(np.complex128) -test_D = test_D.astype(np.complex128) - - class test_ab08n(unittest.TestCase): - """ test1 to 4: Verify ag08bd with input parameters according to example in documentation """ + """ ag08nX with input parameters according to example in documentation """ + + A = np.diag([1., 1., 3., -4., -1., 3.]) + + B = np.array([[ 0., -1.], + [-1., 0.], + [ 1., -1.], + [ 0., 0.], + [ 0., 1.], + [-1., -1.]]) + + C = np.array([[1., 0., 0., 1., 0., 0.], + [0., 1., 0., 1., 0., 1.], + [0., 0., 1., 0., 0., 1.]]) + + D = np.zeros((3, 2)) + + def normalize(self, w): + wi = np.flip(np.argsort(np.abs(w))) + wn = w[wi]/w[wi[0]] + return wn + + def ab08nX(self, ab08fun, A, B, C, D): + n = 6 + m = 2 + p = 3 + # Check the observability and compute the ordered set of + # the observability indices (call the routine with M = 0). + out = ab08fun(n, 0, p, A, B, C, D) + nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf = out[:10] + + assert_equal(kronl[:nkrol], np.array([1, 2, 2])) + assert_equal(n-nu, 5) + assert_allclose(Af[:nu, :nu], np.array([[-1.]])) + # Check the controllability and compute the ordered set of + # the controllability indices (call the routine with P = 0) + out = ab08fun(n, m, 0, A, B, C, D) + nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf = out[:10] + assert_equal(kronr[:nkror], np.array([2, 3])) + assert_equal(n-nu, 5) + assert_allclose(Af[:nu, :nu], np.array([[-4.]])) + # Compute the structural invariants of the given system. + out = ab08fun(n, m, p, A, B, C, D) + nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf = out[:10] + assert_equal(nu, 2) + # Compute the invariant zeros of the given system. + w = eig(Af[:nu, :nu], Bf[:nu, :nu], left=False, right=False) + w_ref = np.array([-2., 1.]) + assert_allclose(self.normalize(w), self.normalize(w_ref)) + # the examples value of infinite zeros does not match the code + # compare output formats to given strings + # assert_equal(sum(infz[:dinfz]), 2) + # assert_equal([[infz[i], i+1] for i in range(dinfz)], [[1, 1]]) + assert_equal(nkror, 0) + assert_equal(nkrol, 1) + assert_equal(kronl[:nkrol], np.array([2])) def test_ab08nd(self): - #test [A-lambda*E] - #B,C,D must have correct dimensions according to l,n,m and p, but cannot have zero length in any dimenstion. Then the wrapper will complain. The length is then set to one. - - nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = analysis.ab08nd(6,2,3,test_A,test_B,test_C,test_D) + self.ab08nX(analysis.ab08nd, self.A, self.B, self.C, self.D) def test_ab08nz(self): - #test [A-lambda*E] - #B,C,D must have correct dimensions according to l,n,m and p, but cannot have zero length in any dimenstion. Then the wrapper will complain. The length is then set to one. - nu,rank,dinfz,nkror,nkrol,infz,kronr,kronl,Af,Bf = analysis.ab08nz(6,2,3,test_A,test_B,test_C,test_D) - - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + Ac, Bc, Cc, Dc = [M.astype(np.complex128) for M in [self.A, self.B, + self.C, self.D]] + self.ab08nX(analysis.ab08nz, Ac, Bc, Cc, Dc) if __name__ == "__main__": From 96ae4bfdfb4e950cd7791971093b1ff9c8ae2b64 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 03:13:05 +0200 Subject: [PATCH 07/10] fix docstring --- slycot/tests/test_ab08n.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/slycot/tests/test_ab08n.py b/slycot/tests/test_ab08n.py index a5773e4e..bdeb0b88 100644 --- a/slycot/tests/test_ab08n.py +++ b/slycot/tests/test_ab08n.py @@ -10,7 +10,8 @@ class test_ab08n(unittest.TestCase): - """ ag08nX with input parameters according to example in documentation """ + """ Test regular pencil construction ab08nX with input parameters + according to example in documentation """ A = np.diag([1., 1., 3., -4., -1., 3.]) @@ -68,9 +69,11 @@ def ab08nX(self, ab08fun, A, B, C, D): assert_equal(kronl[:nkrol], np.array([2])) def test_ab08nd(self): + "Test Construct regular pencil for real matrices" self.ab08nX(analysis.ab08nd, self.A, self.B, self.C, self.D) def test_ab08nz(self): + "Test Construct regular pencil for (pseudo) complex matrices" Ac, Bc, Cc, Dc = [M.astype(np.complex128) for M in [self.A, self.B, self.C, self.D]] self.ab08nX(analysis.ab08nz, Ac, Bc, Cc, Dc) From 96d98908861247b4b6c5bd0f7da86406bda161a4 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 3 Apr 2020 13:24:11 +0200 Subject: [PATCH 08/10] remove complex warning in ab08nz --- slycot/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 96acc378..31bf0238 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -575,7 +575,7 @@ def ab08nz(n, m, p, A, B, C, D, equil='N', tol=0., lzwork=None): e.info = info raise e return (nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf, - int(zwork[0])) + int(zwork[0].real)) def ab09ad(dico,job,equil,n,m,p,A,B,C,nr=None,tol=0,ldwork=None): """ nr,Ar,Br,Cr,hsv = ab09ad(dico,job,equil,n,m,p,A,B,C,[nr,tol,ldwork]) From 021a5c47c5e8dd753a4bc99b61e47b6f904ac5f7 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Fri, 10 Apr 2020 13:15:31 +0200 Subject: [PATCH 09/10] add the test file to CMakeLists.txt --- slycot/tests/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 8c7f26d1..0b021358 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -2,6 +2,7 @@ set(PYSOURCE __init__.py test.py + test_ab08n.py test_ag08bd.py test_mb.py test_mc.py From d397dbbffc7ce033141db99b2a4e141dfbf55897 Mon Sep 17 00:00:00 2001 From: bnavigator Date: Sat, 11 Apr 2020 16:21:40 +0200 Subject: [PATCH 10/10] readd tolerance calculation docstring, reformat lzwork docstring --- slycot/analysis.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/slycot/analysis.py b/slycot/analysis.py index 31bf0238..76c63ba8 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -514,12 +514,21 @@ def ab08nz(n, m, p, A, B, C, D, equil='N', tol=0., lzwork=None): which is defined as the order of the largest leading (or trailing) triangular submatrix in the QR (or RQ) factorization with column (or row) pivoting whose estimated condition number is less than 1/tol. + If tol is set to less than SQRT((N+P)*(N+M))*EPS + then the tolerance is taken as SQRT((N+P)*(N+M))*EPS, + where EPS is the machine precision (see LAPACK Library + Routine DLAMCH). lzwork := None input int - The length of the cache array zwork. The default value is calculated - to MAX( 1, MIN(P,M) + MAX(3*M-1,N), - MIN(P,N) + MAX(3*P-1,N+P,N+M), - MIN(M,N) + MAX(3*M-1,N+M) ) - for optimum performance should be larger. + The length of the internal cache array ZWORK. The default value is + calculated to + MAX( 1, + MIN(P,M) + MAX(3*M-1,N), + MIN(P,N) + MAX(3*P-1,N+P,N+M), + MIN(M,N) + MAX(3*M-1,N+M) ) + For optimum performance lzwork should be larger. + If lzwork = -1, then a workspace query is assumed; + the routine only calculates the optimal size of the + ZWORK array, and returns this value in lzwork_opt Return objects: nu : int The number of (finite) invariant zeros. @@ -577,6 +586,7 @@ def ab08nz(n, m, p, A, B, C, D, equil='N', tol=0., lzwork=None): return (nu, rank, dinfz, nkror, nkrol, infz, kronr, kronl, Af, Bf, int(zwork[0].real)) + def ab09ad(dico,job,equil,n,m,p,A,B,C,nr=None,tol=0,ldwork=None): """ nr,Ar,Br,Cr,hsv = ab09ad(dico,job,equil,n,m,p,A,B,C,[nr,tol,ldwork])