diff --git a/slycot/analysis.py b/slycot/analysis.py index daf05c94..eb18f35e 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -1571,4 +1571,117 @@ def ab13fd(n, A, tol = 0.0): else: raise RuntimeError("unknown error code %r" % out[-1]) +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. + m : input int + The number of columns of matrix B. m >= 0. + p : input int + The number of rows of matrix C. p >= 0. + A : rank-2 array('d') with bounds (l,n) + The leading l-by-n part of this array must + contain the state dynamics matrix A of the system. + E : rank-2 array('d') with bounds (l,n) + The leading l-by-n part of this array must + contain the descriptor matrix E of the system. + B : rank-2 array('d') with bounds (l,m) + The leading l-by-m part of this array must + contain the input/state matrix B of the system. + C : 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 : 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 system + matrix as follows: + = 'S': Perform balancing (scaling); + = 'N': Do not perform balancing. + tol := 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. + If the user sets TOL <= 0, then default tolerances are + used instead, as follows: TOLDEF = L*N*EPS in TG01FD + (to determine the rank of E) and TOLDEF = (L+P)*(N+M)*EPS + in the rest, where EPS is the machine precision + (see LAPACK Library routine DLAMCH). TOL < 1. + ldwork : input int + The length of the cache array. + ldwork >= max( 4*(l,n), ldw ), if equil = 'S', + ldwork >= ldw, if equil = 'N', where + ldw = max(l+p,m+n)*(m+n) + max(1,5*max(l+p,m+n)). + For optimum performance ldwork should be larger. + Return objects: + Af : rank-2 array('d') + the leading NFZ-by-NFZ part of this array + contains the matrix Af of the reduced pencil. + Ef : rank-2 array('d') + the leading NFZ-by-NFZ part of this array + contains the matrix Ef of the reduced pencil. + nrank : output int + The normal rank of the system pencil. + niz : output int + The number of infinite zeros. + infz : rank-1 array('i') + 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 in the Smith form, where i = 1,2,...,DINFZ. + kronr : rank-1 array('i') + The leading NKROR elements of this array contain the + right Kronecker (column) indices. + infe : rank-1 array('i') + The leading NINFE elements of infe contain the + multiplicities of infinite eigenvalues. + """ + 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': + ldwork = max(4*(l+n), ldw) + else: #equil == 'N' + ldwork = ldw + + [Af,Ef,nfz,nrank,niz,dinfz,nkror,ninfe,nkrol,infz,kronr,infe,kronl,info]= _wrapper.ag08bd(equil,l,n,m,p,A,E,B,C,D,tol,ldwork) + + if info < 0: + error_text = "The following argument had an illegal value: "+arg_list[-info-1] + e = ValueError(error_text) + e.info = info + raise e + if info != 0: + e = ArithmeticError('ag08bd failed') + e.info = info + raise e + + return Af[:nfz,:nfz],Ef[:nfz,:nfz],nrank,niz,infz[:dinfz],kronr[:nkror],infe[:ninfe],kronl[:nkrol] + # to be replaced by python wrappers diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index 25c8edb2..58a0656b 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -386,3 +386,36 @@ subroutine ab13fd(n,a,lda,beta,omega,tol,dwork,ldwork,cwork,lcwork,info) ! in AB integer intent(hide),depend(n) :: lcwork = max(1,n*(n+3)) integer intent(out) :: info end subroutine ab13fd +subroutine ag08bd(equil,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,d,ldd,nfz,nrank,niz,dinfz,nkror,ninfe,nkrol,infz,kronr,infe,kronl,tol,iwork,dwork,ldwork,info) ! in AG08BD.f + character intent(in) :: equil + integer intent(in),required,check(l>=0) :: l + integer intent(in),required,check(n>=0) :: n + integer intent(in),required,check(m>=0) :: m + integer intent(in),required,check(p>=0) :: p + double precision intent(in,out,copy),dimension(l,n),depend(l,n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(l,n),depend(l,m) :: e + integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1) + double precision intent(in,copy),dimension(l,m),depend(l,m) :: b + integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1) + double precision intent(in,copy),dimension(p,n),depend(p,n) :: c + integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1) + double precision intent(in,copy),dimension(p,m),depend(p,m) :: d + integer intent(hide),depend(d) :: ldd=MAX(shape(d,0),1) + integer intent(out) :: nfz + integer intent(out) :: nrank + integer intent(out) :: niz + integer intent(out) :: dinfz + integer intent(out) :: nkror + integer intent(out) :: ninfe + integer intent(out) :: nkrol + integer intent(out),dimension(n+1),depend(n) :: infz + integer intent(out),dimension(n+m+1),depend(n,m) :: kronr + integer intent(out),dimension(1+MIN(l+p,n+m)),depend(l,p,n,m) :: infe + integer intent(out),dimension(l+p+1),depend(l,p) :: kronl + double precision intent(in) :: tol + integer intent(cache,hide),dimension(ldwork),depend(ldwork) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork + integer required intent(in) :: ldwork + integer intent(out) :: info +end subroutine ag08bd diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf index 2fdb1308..56e53721 100644 --- a/slycot/src/synthesis.pyf +++ b/slycot/src/synthesis.pyf @@ -525,6 +525,25 @@ subroutine sb10hd(n,m,np,ncon,nmeas,a,lda,b,ldb,c,ldc,d,ldd,ak,ldak,bk,ldbk,ck,l logical intent(hide), dimension(2*n), depend(n) :: bwork integer intent(out) :: info end subroutine sb10hd +subroutine sb10jd(n,m,np,a,lda,b,ldb,c,ldc,d,ldd,e,lde,nsys,dwork,ldwork,info) ! in SB10JD.f + integer intent(in),required :: n + integer intent(in),required :: m + integer intent(in),required :: np + double precision intent(in,out,copy),dimension(lda,n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(ldb,m) :: b + integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1) + double precision intent(in,out,copy),dimension(ldc,n) :: c + integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1) + double precision intent(in,out,copy),dimension(ldd,m) :: d + integer intent(hide),depend(d) :: ldd=MAX(shape(d,0),1) + double precision intent(in,copy),dimension(lde,n) :: e + integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1) + integer intent(out) :: nsys + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork +integer required intent(in) :: ldwork + integer intent(out) :: info +end subroutine sb10jd subroutine sg03ad(dico,job,fact,trans,uplo,n,a,lda,e,lde,q,ldq,z,ldz,x,ldx,scale,sep,ferr,alphar,alphai,beta,iwork,dwork,ldwork,info) ! in SG03AD.f character :: dico character :: job diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf index 819d6247..48c6d44f 100644 --- a/slycot/src/transform.pyf +++ b/slycot/src/transform.pyf @@ -421,4 +421,90 @@ subroutine tb01pd(job,equil,n,m,p,a,lda,b,ldb,c,ldc,nr,tol,iwork,dwork,ldwork,in integer optional,depend(n,m,p) :: ldwork=max(1,n+max(n,max(3*m,3*p))) integer intent(out) :: info end subroutine tb01pd - +subroutine tg01fd_nn(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ldz,ranke,rnka22,tol,iwork,dwork,ldwork,info) ! in TG01FD.f + fortranname tg01fd + character intent(hide) :: compq = 'N' + character intent(hide) :: compz = 'N' + character intent(in),required :: joba + integer intent(in),required,check(l>=0) :: l + integer intent(in),required,check(n>=0) :: n + integer intent(in),required,check(m>=0) :: m + integer intent(in),required,check(p>=0) :: p + double precision intent(in,out,copy),dimension(l,n),depend(l,n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(l,n),depend(l,m) :: e + integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1) + double precision intent(in,out,copy),dimension(l,m),depend(l,m) :: b + integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1) + double precision intent(in,out,copy),dimension(p,n),depend(p,n) :: c + integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1) + double precision intent(hide),dimension(0,0) :: q + integer intent(hide),depend(q) :: ldq=1 + double precision intent(hide),dimension(0,0) :: z + integer intent(hide),depend(z) :: ldz=1 + integer intent(out) :: ranke + integer intent(out) :: rnka22 + double precision intent(in) :: tol + integer intent(cache,hide),dimension(ldwork),depend(ldwork) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork +integer required intent(in) :: ldwork + integer intent(out) :: info +end subroutine tg01fd_nn +subroutine tg01fd_ii(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ldz,ranke,rnka22,tol,iwork,dwork,ldwork,info) ! in TG01FD.f + fortranname tg01fd + character intent(hide) :: compq = 'I' + character intent(hide) :: compz = 'I' + character intent(in),required :: joba + integer intent(in),required,check(l>=0) :: l + integer intent(in),required,check(n>=0) :: n + integer intent(in),required,check(m>=0) :: m + integer intent(in),required,check(p>=0) :: p + double precision intent(in,out,copy),dimension(l,n),depend(l,n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(l,n),depend(l,m) :: e + integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1) + double precision intent(in,out,copy),dimension(l,m),depend(l,m) :: b + integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1) + double precision intent(in,out,copy),dimension(p,n),depend(p,n) :: c + integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1) + double precision intent(out),dimension(l,l) :: q + integer intent(hide),depend(q) :: ldq=MAX(shape(q,0),1) + double precision intent(out),dimension(n,n) :: z + integer intent(hide),depend(z) :: ldz=MAX(shape(z,0),1) + integer intent(out) :: ranke + integer intent(out) :: rnka22 + double precision intent(in) :: tol + integer intent(cache,hide),dimension(ldwork),depend(ldwork) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork +integer required intent(in) :: ldwork + integer intent(out) :: info +end subroutine tg01fd_ii +subroutine tg01fd_uu(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ldz,ranke,rnka22,tol,iwork,dwork,ldwork,info) ! in TG01FD.f + fortranname tg01fd + character intent(hide) :: compq = 'U' + character intent(hide) :: compz = 'U' + character intent(in),required :: joba + integer intent(in),required,check(l>=0) :: l + integer intent(in),required,check(n>=0) :: n + integer intent(in),required,check(m>=0) :: m + integer intent(in),required,check(p>=0) :: p + double precision intent(in,out,copy),dimension(l,n),depend(l,n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(l,n),depend(l,m) :: e + integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1) + double precision intent(in,out,copy),dimension(l,m),depend(l,m) :: b + integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1) + double precision intent(in,out,copy),dimension(p,n),depend(p,n) :: c + integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1) + double precision intent(in,out,copy),dimension(l,l) :: q + integer intent(hide),depend(q) :: ldq=MAX(shape(q,0),1) + double precision intent(in,out,copy),dimension(n,n) :: z + integer intent(hide),depend(z) :: ldz=MAX(shape(z,0),1) + integer intent(out) :: ranke + integer intent(out) :: rnka22 + double precision intent(in) :: tol + integer intent(cache,hide),dimension(ldwork),depend(ldwork) :: iwork + double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork +integer required intent(in) :: ldwork + integer intent(out) :: info +end subroutine tg01fd_uu diff --git a/slycot/synthesis.py b/slycot/synthesis.py index 75daa347..2c0e442f 100644 --- a/slycot/synthesis.py +++ b/slycot/synthesis.py @@ -1838,6 +1838,85 @@ def sb10hd(n,m,np,ncon,nmeas,A,B,C,D,tol=0.0,ldwork=None): raise e return out[:-1] + +def sb10jd(n,m,np,A,B,C,D,E,ldwork=None): + """ A,B,C,D = sb10jd(n,m,np,A,B,C,D,E,[ldwork]) + + To convert the descriptor state-space system + + E*dx/dt = A*x + B*u + y = C*x + D*u + + into regular state-space form + + dx/dt = Ad*x + Bd*u + y = Cd*x + Dd*u . + + Required arguments: + n : input int + The order of the descriptor system. n >= 0. + m : input int + The column size of the matrix B. m >= 0. + np : input int + The row size of the matrix C. np >= 0. + A : rank-2 array('d') with bounds (n,n) + The leading n-by-n part of this array must + contain the state matrix A of the descriptor system. + B : rank-2 array('d') with bounds (l,m) + The leading n-by-m part of this array must + contain the input matrix B of the descriptor system. + C : rank-2 array('d') with bounds (np,n) + The leading np-by-n part of this array must + contain the output matrix C of the descriptor system. + D : rank-2 array('d') with bounds (np,m) + The leading np-by-m part of this array must + contain the matrix D of the descriptor system. + E : rank-2 array('d') with bounds (l,n) + The leading n-by-n part of this array must + contain the matrix E of the descriptor system. + Optional arguments: + ldwork : input int + The length of the cache array. + ldwork >= max( 1, 2*n*n + 2*n + n*MAX( 5, n + m + np ) ). + For good performance, ldwork must generally be larger. + Return objects: + A : rank-2 array('d') with bounds (nsys,nsys) + The leading nsys-by-nsys part of this array + contains the state matrix Ad of the converted system. + B : rank-2 array('d') with bounds (nsys,m) + The leading NSYS-by-M part of this array + contains the input matrix Bd of the converted system. + C : rank-2 array('d') with bounds (np,nsys) + The leading NP-by-NSYS part of this array + contains the output matrix Cd of the converted system. + D : rank-2 array('d') with bounds (np,m) + The leading NP-by-M part of this array contains + the matrix Dd of the converted system. + """ + + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'm', 'np', 'A', 'lda'+hidden, 'B', 'ldb'+hidden, 'C', 'ldc'+hidden, 'D', 'ldd'+hidden, 'E', 'lde'+hidden, 'nsys', 'dwork'+hidden, 'ldwork', 'info'] + + if ldwork is None: + ldwork = max(1, 2 * n * n + 2 * n + n * max(5, n + m + np)) + + A,B,C,D,nsys,info = _wrapper.sb10jd(n,m,np,A,B,C,D,E,ldwork) + + if info < 0: + error_text = "The following argument had an illegal value: "+arg_list[-info-1] + e = ValueError(error_text) + e.info = info + raise e + elif info == 1: + e = ArithmeticError("The sb10jd algorithm did not converge") + e.info = 1 + raise e + elif info != 0: + e = ArithmeticError('sb10jd failed') + e.info = info + raise e + + return A[:nsys,:nsys],B[:nsys,:m],C[:np, :nsys],D[:np, :m] def sg03ad(dico,job,fact,trans,uplo,N,A,E,Q,Z,X,ldwork=None): """ A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = diff --git a/slycot/tests/test_ag08bd.py b/slycot/tests/test_ag08bd.py new file mode 100644 index 00000000..f1101765 --- /dev/null +++ b/slycot/tests/test_ag08bd.py @@ -0,0 +1,118 @@ +# =================================================== +# ag08bd tests + +import unittest +from slycot import analysis +import numpy as np + +from numpy.testing import assert_raises, assert_almost_equal, assert_equal + +# test1 input parameters + +test1_l = 9 +test1_n = 9 +test1_m = 3 +test1_p = 3 +test1_tol = 1.0e-7 +test1_equil = 'N' + +test1_A = np.eye(9, dtype=int) + +test1_E = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0]]) + +test1_B = np.array([[-1, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, -1, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, -1], + [ 0, 0, 0], + [ 0, 0, 0]]) + +test1_C = np.array([[ 0, 1, 1, 0, 3, 4, 0, 0, 2], + [ 0, 1, 0, 0, 4, 0, 0, 2, 0], + [ 0, 0, 1, 0, -1, 4, 0, -2, 2]]) + +test1_D = np.array([[ 1, 2, -2], + [ 0, -1, -2], + [ 0, 0, 0]]) + + +class test_tg01fd(unittest.TestCase): + """ test1 to 4: Verify ag08bd with input parameters according to example in documentation """ + + def test1_ag08bd(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. + + Af,Ef,nrank,niz,infz,kronr,infe,kronl = analysis.ag08bd(l=test1_l,n=test1_n,m=0,p=0,A=test1_A,E=test1_E,B=np.zeros((test1_l,1)),C=np.zeros((1,test1_n)),D=np.zeros((1,1)),equil=test1_equil, tol=test1_tol) + + assert_equal(Af, np.zeros((0,0))) + assert_equal(Ef, np.zeros((0,0))) + assert_equal(nrank, 9) + assert_equal(niz, 6) + assert_equal(infz, [0,3]) + assert_equal(kronr, []) + assert_equal(infe, [3,3,3]) + assert_equal(kronl, []) + + def test2_ag08bd(self): + #test [A-lambda*E;C] + #B,D must have correct dimensions as before + + Af,Ef,nrank,niz,infz,kronr,infe,kronl = analysis.ag08bd(l=test1_l,n=test1_n,m=0,p=test1_p,A=test1_A,E=test1_E,B=np.zeros((test1_l,1)),C=test1_C,D=np.zeros((test1_p,1)),equil=test1_equil, tol=test1_tol) + + assert_equal(Af, np.zeros((0,0))) + assert_equal(Ef, np.zeros((0,0))) + assert_equal(nrank, 9) + assert_equal(niz, 4) + assert_equal(infz, [0,2]) + assert_equal(kronr, []) + assert_equal(infe, [1,3,3]) + assert_equal(kronl, [0,1,1]) + + def test3_ag08bd(self): + #test [A-lambda*E,B] + #C,D must have correct dimensions as before + + Af,Ef,nrank,niz,infz,kronr,infe,kronl = analysis.ag08bd(l=test1_l,n=test1_n,m=test1_m,p=0,A=test1_A,E=test1_E,B=test1_B,C=np.zeros((1,test1_n)),D=np.zeros((1,test1_m)),equil=test1_equil, tol=test1_tol) + + assert_equal(Af, np.zeros((0,0))) + assert_equal(Ef, np.zeros((0,0))) + assert_equal(nrank, 9) + assert_equal(niz, 0) + assert_equal(infz, []) + assert_equal(kronr, [2,2,2]) + assert_equal(infe, [1,1,1]) + assert_equal(kronl, []) + + def test4_ag08bd(self): + #test [A-lambda*E,B;C,D] + + Af,Ef,nrank,niz,infz,kronr,infe,kronl = analysis.ag08bd(l=test1_l,n=test1_n,m=test1_m,p=test1_p,A=test1_A,E=test1_E,B=test1_B,C=test1_C,D=test1_D,equil=test1_equil, tol=test1_tol) + + assert_almost_equal(Af, [[0.77045021]]) + assert_almost_equal(Ef, [[0.77045021]]) + assert_equal(nrank, 11) + assert_equal(niz, 2) + assert_equal(infz, [0,1]) + assert_equal(kronr, [2]) + assert_equal(infe, [1,1,1,1,3]) + assert_equal(kronl, [1]) + + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/tests/test_sb10jd.py b/slycot/tests/test_sb10jd.py new file mode 100644 index 00000000..c252dcee --- /dev/null +++ b/slycot/tests/test_sb10jd.py @@ -0,0 +1,72 @@ +# =================================================== +# sb10jd tests + +import unittest +from slycot import synthesis +import numpy as np +from numpy.testing import assert_raises, assert_almost_equal, assert_equal + +# test1 input parameters + +test1_n = 6 +test1_m = 1 +test1_np = 6 + +test1_A = np.array([[ 0, 0, 0, -1, 1, 0], + [ 0, 32, 0, 0, -1, 1], + [ 0, 0, 1, 0, 0, 0], + [ 0, 0, 0, 1, 0, 0], + [-1, 1, 0, 0, 0, 0], + [ 0, -1, 1, 0, 0, 0]]) + + +test1_E = np.array([[ 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, -10, 0, 10], + [ 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0]]) + +test1_B = np.array([[-7.1], + [ 0. ], + [ 0. ], + [ 0. ], + [ 0. ], + [ 0. ]]) + +test1_C = np.eye(6) + +test1_D = np.zeros((7,1)) + +# test1 expected results + +test1_Aexp = np.array([[-0.00312500]]) +test1_Bexp = np.array([[ 0.05899985]]) +test1_Cexp = np.array([[-1.17518847e-02], + [-1.17518847e-02], + [-1.17518847e-02], + [ 0.00000000e+00], + [ 0.00000000e+00], + [ 3.76060309e-01]]) +test1_Dexp = np.array([[ 2.21875000e-01], + [ 2.21875000e-01], + [ 2.21875000e-01], + [ 0.00000000e+00], + [ 7.10000000e+00], + [ 0.00000000e+00]]) + +class test_sb10jd(unittest.TestCase): + def test1_sb10jd(self): + """ verify the output of sb10jd for a descriptor system """ + A,B,C,D = synthesis.sb10jd(test1_n,test1_m,test1_np,test1_A,test1_B,test1_C,test1_D,test1_E) + assert_almost_equal(A, test1_Aexp) + assert_almost_equal(B, test1_Bexp) + assert_almost_equal(C, test1_Cexp) + assert_almost_equal(D, test1_Dexp) + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/tests/test_tg01fd.py b/slycot/tests/test_tg01fd.py new file mode 100644 index 00000000..bc8fe2cb --- /dev/null +++ b/slycot/tests/test_tg01fd.py @@ -0,0 +1,115 @@ +# =================================================== +# tg01fd tests + +import unittest +from slycot import transform +import numpy as np + +from numpy.testing import assert_raises, assert_almost_equal, assert_equal + +# test1 input parameters +test1_l = 4 +test1_n = 4 +test1_m = 2 +test1_p = 2 +test1_tol = 0.0 +test1_A = np.array([[-1, 0, 0, 3], + [ 0, 0, 1, 2], + [ 1, 1, 0, 4], + [ 0, 0, 0, 0]]) + +test1_E = np.array([[1, 2, 0, 0], + [0, 1, 0, 1], + [3, 9, 6, 3], + [0, 0, 2, 0]]) + +test1_B = np.array([[1, 0], + [0, 0], + [0, 1], + [1, 1]]) + +test1_C = np.array([[-1, 0, 1, 0], + [ 0, 1, -1, 1]]) + +#test1 expected output +test1_Aexp = np.array([[ 2.02781052, 0.10783277, 3.90616686, -2.15710472], + [-0.09804588, 0.25437761, 1.60529591, -0.12692683], + [ 0.27131089, 0.77603837, -0.36920735, -0.48533567], + [ 0.06900656, -0.56694671, -2.19740106, 0.3086067 ]]) + +test1_Eexp = np.array([[10.15874008, 5.82296975, 1.30205562, 0. ], + [ 0. , -2.468405 , -0.18960188, 0. ], + [ 0. , 0. , 1.03378058, 0. ], + [ 0. , 0. , 0. , 0. ]]) + +test1_Bexp = np.array([[-0.21566555, -0.97049496], + [ 0.30148458, 0.95156071], + [ 0.75952691, 0.09906873], + [ 1.13389342, 0.37796447]]) + +test1_Cexp = np.array([[ 3.65148372e-01, -1.00000000e+00, -4.47213595e-01, -8.16496581e-01], + [-1.09544512e+00, 1.00000000e+00, -8.94427191e-01, 2.22044605e-16]]) + +test1_Qexp = np.array([[-0.21566555, -0.50875523, 0.61092382, 0.56694671], + [-0.10783277, -0.25437761, -0.77603837, 0.56694671], + [-0.97049496, 0.1413209 , -0.04953436, -0.18898224], + [ 0. , 0.81023981, 0.14860309, 0.56694671]]) +test1_Zexp = np.array([[-3.65148372e-01, -1.35772740e-16, 4.47213595e-01, 8.16496581e-01], + [-9.12870929e-01, 0.00000000e+00, 0.00000000e+00, -4.08248290e-01], + [ 6.19714937e-17, -1.00000000e+00, 0.00000000e+00, -1.38572473e-16], + [-1.82574186e-01, -6.78863700e-17, -8.94427191e-01, 4.08248290e-01]]) + +test1_ranke_exp = 3 +test1_rnka22_exp = 1 + +class test_tg01fd(unittest.TestCase): + + def test1_tg01fd(self): + """ test1: Verify from tg01fd with input parameters according to test in documentation """ + A,E,B,C,ranke,rnka22,Q,Z = transform.tg01fd(l=test1_l,n=test1_n,m=test1_m,p=test1_p,A=test1_A,E=test1_E,B=test1_B,C=test1_C,compq='I',compz='I',joba='T',tol=test1_tol) + assert_almost_equal(A, test1_Aexp) + assert_almost_equal(E, test1_Eexp) + assert_almost_equal(B, test1_Bexp) + assert_almost_equal(C, test1_Cexp) + assert_almost_equal(Q, test1_Qexp) + assert_almost_equal(Z, test1_Zexp) + assert_equal(test1_ranke_exp, ranke) + assert_equal(test1_rnka22_exp, rnka22) + + def test2_tg01fd(self): + """ verify that Q and Z output with compq and compz set to 'U' equals the dot product of Q and Z input and Q and Z output with compq and compz set to 'I' """ + + l = 30 + n = 30 + m = 70 + p = 44 + + np.random.seed(0) + + Ain = np.random.rand(l, n) + Ein = np.random.rand(l, n) + Bin = np.random.rand(n, m) + Cin = np.random.rand(p, n) + Qin = np.random.randn(l,l) + Zin = np.random.randn(n,n) + + A_1,E_1,B_1,C_1,ranke_1,rnka22_1,Q_1,Z_1= transform.tg01fd(l=l,n=n,m=m,p=p,A=Ain,E=Ein,B=Bin,C=Cin,compq='I', compz='I', joba='T', tol=0.0) + + A_2,E_2,B_2,C_2,ranke_2,rnka22_2,Q_2,Z_2= transform.tg01fd(l=l,n=n,m=m,p=p,A=Ain,E=Ein,B=Bin,C=Cin,Q=Qin,Z=Zin,compq='U', compz='U', joba='T', tol=0.0) + + assert_equal(A_1, A_2) + assert_equal(E_1, E_2) + assert_equal(B_1, B_2) + assert_equal(C_1, C_2) + assert_equal(ranke_1, ranke_2) + assert_equal(rnka22_1, rnka22_2) + + assert_almost_equal(np.dot(Qin, Q_1), Q_2) + assert_almost_equal(np.dot(Zin, Z_1), Z_2) + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/transform.py b/slycot/transform.py index 485f2f3d..f3731893 100644 --- a/slycot/transform.py +++ b/slycot/transform.py @@ -1050,5 +1050,199 @@ def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None): raise e return out[:-1] +def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ldwork=None): + """ A,E,B,C,ranke,rnka22,Q,Z = tg01fd(l,n,m,p,A,E,B,C,[Q,Z,compq,compz,joba,tol,ldwork]) + + To compute for the descriptor system (A-lambda E,B,C) + the orthogonal transformation matrices Q and Z such that the + transformed system (Q'*A*Z-lambda Q'*E*Z, Q'*B, C*Z) is + in a SVD-like coordinate form with + + ( A11 A12 ) ( Er 0 ) + Q'*A*Z = ( ) , Q'*E*Z = ( ) , + ( A21 A22 ) ( 0 0 ) + + where Er is an upper triangular invertible matrix. + Optionally, the A22 matrix can be further reduced to the form + + ( Ar X ) + A22 = ( ) , + ( 0 0 ) + + with Ar an upper triangular invertible matrix, and X either a full + or a zero matrix. + The left and/or right orthogonal transformations performed + to reduce E and A22 can be optionally accumulated. + + 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. + m : input int + The number of columns of matrix B. m >= 0. + p : input int + The number of rows of matrix C. p >= 0. + A : rank-2 array('d') with bounds (l,n) + The leading l-by-n part of this array must + contain the state dynamics matrix A. + E : rank-2 array('d') with bounds (l,n) + The leading l-by-n part of this array must + contain the descriptor matrix E. + B : rank-2 array('d') with bounds (l,m) + The leading L-by-M part of this array must + contain the input/state matrix B. + C : rank-2 array('d') with bounds (p,n) + The leading P-by-N part of this array must + contain the state/output matrix C. + Optional arguments: + Q : rank-2 array('d') with bounds (l,l) + If COMPQ = 'N': Q is not referenced. + If COMPQ = 'I': Q need not be set. + If COMPQ = 'U': The leading l-by-l part of this + array must contain an orthogonal matrix + Q1. + Z : rank-2 array('d') with bounds (n,n) + If COMPZ = 'N': Z is not referenced. + If COMPZ = 'I': Z need not be set. + If COMPZ = 'U': The leading n-by-n part of this + array must contain an orthogonal matrix + Z1. + compq := 'N' input string(len=1) + = 'N': do not compute Q. + = 'I': Q is initialized to the unit matrix, and the + orthogonal matrix Q is returned. + = 'U': Q must contain an orthogonal matrix Q1 on entry, + and the product Q1*Q is returned. + compz := 'N' input string(len=1) + = 'N': do not compute Z. + = 'I': Z is initialized to the unit matrix, and the + orthogonal matrix Z is returned. + = 'U': Z must contain an orthogonal matrix Z1 on entry, + and the product Z1*Z is returned. + joba := 'N' input string(len=1) + = 'N': do not reduce A22. + = 'R': reduce A22 to a SVD-like upper triangular form. + = 'T': reduce A22 to an upper trapezoidal form. + tol := 0 input float + The tolerance to be used in determining the rank of E + and of A22. If the user sets TOL > 0, then the given + value of TOL is used as a lower bound for the + reciprocal condition numbers of leading submatrices + of R or R22 in the QR decompositions E * P = Q * R of E + or A22 * P22 = Q22 * R22 of A22. + A submatrix whose estimated condition number is less than + 1/TOL is considered to be of full rank. If the user sets + TOL <= 0, then an implicitly computed, default tolerance, + defined by TOLDEF = L*N*EPS, is used instead, where + EPS is the machine precision (see LAPACK Library routine + DLAMCH). TOL < 1. + ldwork : input int + The length of the cache array. + ldwork >= MAX( 1, n+p, MIN(l,n)+MAX(3*n-1,m,l) ). + For optimal performance, ldwork should be larger. + Return objects: + A : rank-2 array('d') with bounds (l,n) + On entry, the leading L-by-N part of this array must + contain the state dynamics matrix A. + On exit, the leading L-by-N part of this array contains + the transformed matrix Q'*A*Z. If JOBA = 'T', this matrix + is in the form + + ( A11 * * ) + Q'*A*Z = ( * Ar X ) , + ( * 0 0 ) + + where A11 is a RANKE-by-RANKE matrix and Ar is a + RNKA22-by-RNKA22 invertible upper triangular matrix. + If JOBA = 'R' then A has the above form with X = 0. + E : rank-2 array('d') with bounds (l,n) + The leading L-by-N part of this array contains + the transformed matrix Q'*E*Z. + + ( Er 0 ) + Q'*E*Z = ( ) , + ( 0 0 ) + + where Er is a RANKE-by-RANKE upper triangular invertible + matrix. + B : rank-2 array('d') with bounds (l,m) + The leading L-by-M part of this array contains + the transformed matrix Q'*B. + C : rank-2 array('d') with bounds (p,n) + The leading P-by-N part of this array contains + the transformed matrix C*Z. + Q : rank-2 array('d') with bounds (l,l) + If COMPQ = 'N': Q is not referenced. + If COMPQ = 'I': The leading L-by-L part of this + array contains the orthogonal matrix Q, + where Q' is the product of Householder + transformations which are applied to A, + E, and B on the left. + If COMPQ = 'U': The leading L-by-L part of this + array contains the orthogonal matrix + Q1*Q. + Z : rank-2 array('d') with bounds (n,n) + If COMPZ = 'N': Z is not referenced. + If COMPZ = 'I': The leading N-by-N part of this + array contains the orthogonal matrix Z, + which is the product of Householder + transformations applied to A, E, and C + on the right. + If COMPZ = 'U': The leading N-by-N part of this + array contains the orthogonal matrix + Z1*Z. + ranke : output int + The estimated rank of matrix E, and thus also the order + of the invertible upper triangular submatrix Er. + rnka22 : output int + If JOBA = 'R' or 'T', then RNKA22 is the estimated rank of + matrix A22, and thus also the order of the invertible + upper triangular submatrix Ar. + If JOBA = 'N', then RNKA22 is not referenced. + """ + + hidden = ' (hidden by the wrapper)' + arg_list = ['compq', 'compz', 'joba', 'l', 'n', 'm', 'p', 'A', 'lda'+hidden, 'E','lde'+hidden,'B','ldb'+hidden,'C','ldc'+hidden,'Q','ldq'+hidden,'Z','ldz'+hidden,'ranke','rnka22','tol','iwork'+hidden, 'dwork'+hidden, 'ldwork', 'info'] + + + if compq != 'N' and compq != 'I' and compq != 'U': + raise ValueError('Parameter compq had an illegal value') + + if compz != 'N' and compz != 'I' and compz != 'U': + raise ValueError('Parameter compz had an illegal value') + + if joba != 'N' and joba != 'R' and joba != 'T': + raise ValueError('Parameter joba had an illegal value') + + if ldwork is None: + ldwork = max(1, n+p, min(l,n) + max(3*n-1, m, l)) + + + if compq == 'N' and compz == 'N': + A,E,B,C,ranke,rnka22,info = _wrapper.tg01fd_nn(joba,l,n,m,p,A,E,B,C,tol,ldwork) + Q = None + Z = None + elif compq == 'I' and compz == 'I': + A,E,B,C,Q,Z,ranke,rnka22,info = _wrapper.tg01fd_ii(joba,l,n,m,p,A,E,B,C,tol,ldwork) + elif compq == 'U' and compz == 'U': + A,E,B,C,Q,Z,ranke,rnka22,info = _wrapper.tg01fd_uu(joba,l,n,m,p,A,E,B,C,Q,Z,tol,ldwork) + else: + raise ValueError("The combination of compq and compz in not implemented") + + if info < 0: + error_text = "The following argument had an illegal value: "+arg_list[-info-1] + e = ValueError(error_text) + e.info = info + raise e + if info != 0: + e = ArithmeticError('tg01fd failed') + e.info = info + raise e + + if joba == 'N': + rnka22 = None + + return A,E,B,C,ranke,rnka22,Q,Z # to be replaced by python wrappers