diff --git a/README.rst b/README.rst index f4b11619..75c583b2 100644 --- a/README.rst +++ b/README.rst @@ -1,10 +1,10 @@ Slycot ============= -.. image:: https://travis-ci.org/jgoppert/Slycot.svg - :target: https://travis-ci.org/jgoppert/Slycot -.. image:: https://coveralls.io/repos/jgoppert/Slycot/badge.png - :target: https://coveralls.io/r/jgoppert/Slycot +.. image:: https://travis-ci.org/python-control/Slycot.svg?branch=master + :target: https://travis-ci.org/python-control/Slycot +.. image:: https://coveralls.io/repos/python-control/Slycot/badge.png + :target: https://coveralls.io/r/python-control/Slycot Python wrapper for selected SLICOT routines, notably including solvers for Riccati, Lyapunov and Sylvester equations. diff --git a/slycot/__init__.py b/slycot/__init__.py index 6595c40d..96b54b31 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -25,8 +25,8 @@ # Mathematical routines (3/81 wrapped) from .math import mc01td, mb05md, mb05nd - # Synthesis routines (11/50 wrapped) - from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md + # Synthesis routines (12/50 wrapped) + from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od from .synthesis import sb04md,sb04qd,sb10ad,sb10hd,sg03ad from .synthesis import sg02ad diff --git a/slycot/analysis.py b/slycot/analysis.py index 86a818de..44391aa9 100644 --- a/slycot/analysis.py +++ b/slycot/analysis.py @@ -876,7 +876,7 @@ def ab09bd(dico,job,equil,n,m,p,A,B,C,D,nr=None,tol1=0,tol2=0,ldwork=None): e = ArithmeticError('The computation of Hankel singular values failed') e.info = out[-1] raise e - nr,A,B,C,D,hsv = out[:-2] + Nr,A,B,C,D,hsv = out[:-2] return nr, A[:Nr,:Nr], B[:Nr,:], C[:,:Nr],D[:,:], hsv # to be replaced by python wrappers diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf index e89ec8bc..f86773fb 100644 --- a/slycot/src/analysis.pyf +++ b/slycot/src/analysis.pyf @@ -231,4 +231,30 @@ subroutine ab09ax(dico,job,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,hsv,t,ldt,ti,ldti,t integer optional :: ldwork = max(1,n*(max(n,max(m,p))+5)+n*(n+1)/2) integer intent(out) :: iwarn integer intent(out) :: info -end subroutine ab09ax \ No newline at end of file +end subroutine ab09ax +subroutine ab09bd(dico,job,equil,ordsel,n,m,p,nr,a,lda,b,ldb,c,ldc,d,ldd,hsv,tol1,tol2,iwork,dwork,ldwork,iwarn,info) !in :balred:AB09BD.f + character intent(in) :: dico + character intent(in) :: job + character intent(in) :: equil + character intent(in) :: ordsel + integer check(n>=0) :: n + integer check(m>=0) :: m + integer check(p>=0) :: p + integer intent(in,out) :: nr + 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) + double precision intent(out),dimension(n),depend(n) :: hsv + double precision :: tol1 =0.0 + double precision :: tol2 =0.0 + integer intent(hide,cache),dimension(max(m,p)) :: iwork + double precision intent(hide,cache),dimension(ldwork) :: dwork + integer optional :: ldwork = max(1,n*(2*n+max(n,max(m,p))+5)+n*(n+1)/2) + integer intent(out) :: iwarn + integer intent(out) :: info +end subroutine ab09bd diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf index fb035932..db78bcce 100644 --- a/slycot/src/synthesis.pyf +++ b/slycot/src/synthesis.pyf @@ -347,6 +347,34 @@ subroutine sb03md(dico,job,fact,trana,n,c,ldc,a,lda,u,ldu,scale,sep,ferr,wr,wi,i integer optional,check(ldwork>=max(2*n*n,3*n)),depend(n) :: ldwork=max(2*n*n,3*n) integer intent(out) :: info end subroutine sb03md +<<<<<<< HEAD +<<<<<<< HEAD +subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info) ! in SB03OD.f +======= +subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info) +>>>>>>> f1b51c5... Fixed sb03od bug in synthesis.pyf. +======= +subroutine sb03od(dico,fact,trans,n,m,a,lda,q,ldq,b,ldb,scale,wr,wi,dwork,ldwork,info) ! in SB03OD.f +>>>>>>> 47a3169... Added wrapper to SLICOT routine SB03OD. + fortranname sb03od + character :: dico + character :: fact='N' + character :: trans='N' + integer check(n>0) :: n + integer check(m>0) :: m + double precision dimension(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=shape(a,0) + double precision dimension(n,n),depend(n) :: q + integer intent(hide),depend(q) :: ldq=shape(q,0) + double precision intent(in,out,copy),dimension(n,n),depend(n) :: b + integer intent(hide),depend(b) :: ldb=shape(b,0) + double precision intent(out) :: scale + double precision intent(out),dimension(n),depend(n) :: wr + double precision intent(out),dimension(n),depend(n) :: wi + double precision intent(hide,cache),dimension(ldwork) :: dwork + integer optional,check(ldwork>=max(1,4*n + min(m,n))),depend(n,m) :: ldwork=max(1,4*n + min(m,n)) + integer intent(out) :: info +end subroutine sb03od subroutine sb04md(n,m,a,lda,b,ldb,c,ldc,z,ldz,iwork,dwork,ldwork,info) ! in SB04MD.f integer check(n>0) :: n integer check(m>0) :: m diff --git a/slycot/synthesis.py b/slycot/synthesis.py index da8ee69d..f088e67a 100644 --- a/slycot/synthesis.py +++ b/slycot/synthesis.py @@ -881,6 +881,252 @@ def sb03md(n,C,A,U,dico,job='X',fact='N',trana='N',ldwork=None): w.imag = wi[0:n] return X,scale,sep,ferr,w +def sb03od(n,m,A,Q,B,dico,fact='N',trans='N',ldwork=None): + """ U,scale,w = sb03od(dico,n,m,A,Q,B,[fact,trans,ldwork]) + + To solve for X = op(U)'*op(U) either the stable non-negative + definite continuous-time Lyapunov equation + 2 + op(A)'*X + X*op(A) = -scale *op(B)'*op(B) (1) + + or the convergent non-negative definite discrete-time Lyapunov + equation + 2 + op(A)'*X*op(A) - X = -scale *op(B)'*op(B) (2) + + where op(K) = K or K' (i.e., the transpose of the matrix K), A is + an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper + triangular matrix containing the Cholesky factor of the solution + matrix X, X = op(U)'*op(U), and scale is an output scale factor, + set less than or equal to 1 to avoid overflow in X. If matrix B + has full rank then the solution matrix X will be positive-definite + and hence the Cholesky factor U will be nonsingular, but if B is + rank deficient then X may be only positive semi-definite and U + will be singular. + + In the case of equation (1) the matrix A must be stable (that + is, all the eigenvalues of A must have negative real parts), + and for equation (2) the matrix A must be convergent (that is, + all the eigenvalues of A must lie inside the unit circle). + + Required arguments + ------------------ + + n : input int + The order of the matrix A and the number of columns in + matrix op(B). n >= 0. + m : input int + The number of rows in matrix op(B). m >= 0. + A : input rank-2 array('d'), shape (n,n) + On entry, the leading n-by-n part of this array must + contain the matrix A. If fact = 'F', then A contains + an upper quasi-triangular matrix S in Schur canonical + form; the elements below the upper Hessenberg part of the + array A are not referenced. + On exit, the leading n-by-n upper Hessenberg part of this + array contains the upper quasi-triangular matrix S in + Schur canonical form from the Shur factorization of A. + The contents of array A is not modified if fact = 'F'. + Q : input rank-2 array('d'), shape (n,n) + On entry, if fact = 'F', then the leading n-by-n part of + this array must contain the orthogonal matrix Q of the + Schur factorization of A. + Otherwise, Q need not be set on entry. + On exit, the leading n-by-n part of this array contains + the orthogonal matrix Q of the Schur factorization of A. + The contents of array Q is not modified if fact = 'F'. + B : input rank-2 array('d'), shape (m,n) + On entry, if trans = 'N', the leading m-by-n part of this + array must contain the coefficient matrix B of the + equation. + On entry, if trans = 'T', the leading N-by-m part of this + array must contain the coefficient matrix B of the + equation. + On exit, the leading n-by-n part of this array contains + the upper triangular Cholesky factor U of the solution + matrix X of the problem, X = op(U)'*op(U). + If m = 0 and n > 0, then U is set to zero. + dico : input string(len=1) + Specifies the type of Lyapunov equation to be solved as + follows: + = 'C': Equation (1), continuous-time case; + = 'D': Equation (2), discrete-time case. + + Optional arguments + ------------------ + + fact := 'N' input string(len=1) + Specifies whether or not the real Schur factorization + of the matrix A is supplied on entry, as follows: + = 'F': On entry, A and Q contain the factors from the + real Schur factorization of the matrix A; + = 'N': The Schur factorization of A will be computed + and the factors will be stored in A and Q. + trans := 'N' input string(len=1) + Specifies the form of op(K) to be used, as follows: + = 'N': op(K) = K (No transpose); + = 'T': op(K) = K**T (Transpose). + ldwork := None input int + The length of the array DWORK. + If m > 0, ldwork >= max(1,4*n + min(m,n)); + If m = 0, ldwork >= 1. + For optimum performance ldwork should sometimes be larger. + + Return objects + ______________ + + U : rank-2 array('d'), shape (n,n) + The leading n-by-n part of this array contains + the upper triangular Cholesky factor U of the solution + matrix X of the problem, X = op(U)'*op(U). + scale : float + The scale factor, scale, set less than or equal to 1 to + prevent the solution overflowing. + w : rank-1 array('c'), shape (n) + If fact = 'N', this array contains the eigenvalues of A. + + Raises + ______ + + ValueError : e + e.info contains information about the exact type of exception + < 0: if info = -i, the i-th argument had an illegal value; + = 1: if the Lyapunov equation is (nearly) singular + (warning indicator); + if DICO = 'C' this means that while the matrix A + (or the factor S) has computed eigenvalues with + negative real parts, it is only just stable in the + sense that small perturbations in A can make one or + more of the eigenvalues have a non-negative real + part; + if DICO = 'D' this means that while the matrix A + (or the factor S) has computed eigenvalues inside + the unit circle, it is nevertheless only just + convergent, in the sense that small perturbations + in A can make one or more of the eigenvalues lie + outside the unit circle; + perturbed values were used to solve the equation; + = 2: if FACT = 'N' and DICO = 'C', but the matrix A is + not stable (that is, one or more of the eigenvalues + of A has a non-negative real part), or DICO = 'D', + but the matrix A is not convergent (that is, one or + more of the eigenvalues of A lies outside the unit + circle); however, A will still have been factored + and the eigenvalues of A returned in WR and WI. + = 3: if FACT = 'F' and DICO = 'C', but the Schur factor S + supplied in the array A is not stable (that is, one + or more of the eigenvalues of S has a non-negative + real part), or DICO = 'D', but the Schur factor S + supplied in the array A is not convergent (that is, + one or more of the eigenvalues of S lies outside the + unit circle); + = 4: if FACT = 'F' and the Schur factor S supplied in + the array A has two or more consecutive non-zero + elements on the first sub-diagonal, so that there is + a block larger than 2-by-2 on the diagonal; + = 5: if FACT = 'F' and the Schur factor S supplied in + the array A has a 2-by-2 diagonal block with real + eigenvalues instead of a complex conjugate pair; + = 6: if FACT = 'N' and the LAPACK Library routine DGEES + has failed to converge. This failure is not likely + to occur. The matrix B will be unaltered but A will + be destroyed. + """ + hidden = ' (hidden by the wrapper)' + arg_list = ['dico','fact', 'trans', 'n', 'm', 'a', 'lda'+hidden, 'q', + 'ldq'+hidden, 'b', 'ldb'+hidden, 'scale', 'wr'+hidden, + 'wi'+hidden, 'dwork'+hidden, 'ldwork', 'info'+hidden] + if ldwork is None: + if m > 0: + ldwork = max(1,4*n + min(m,n)) + elif m == 0: + ldwork = 1 + if dico != 'C' and dico != 'D': + raise ValueError('dico must be either D or C') + out = _wrapper.sb03od(dico,n,m,A,Q,B,fact=fact,trans=trans,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 + if out[-1] == 1: + if dico == 'D': + error_text = """this means that while the matrix A + (or the factor S) has computed eigenvalues inside + the unit circle, it is nevertheless only just + convergent, in the sense that small perturbations + in A can make one or more of the eigenvalues lie + outside the unit circle; + perturbed values were used to solve the equation;""" + else: + error_text = """this means that while the matrix A + (or the factor S) has computed eigenvalues with + negative real parts, it is only just stable in the + sense that small perturbations in A can make one or + more of the eigenvalues have a non-negative real + part;""" + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 2: + if dico == 'D': + error_text = """the matrix A is not convergent (that is, one or + more of the eigenvalues of A lies outside the unit + circle); however, A will still have been factored + and the eigenvalues of A returned in WR and WI.""" + else: + error_text = """the matrix A is + not stable (that is, one or more of the eigenvalues + of A has a non-negative real part).""" + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 3: + if dico == 'D': + error_text = """the Schur factor S + supplied in the array A is not convergent (that is, + one or more of the eigenvalues of S lies outside the + unit circle).""" + else: + error_text = """the Schur factor S + supplied in the array A is not stable (that is, one + or more of the eigenvalues of S has a non-negative + real part).""" + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 4: + if fact == 'F': + error_text = """the Schur factor S supplied in + the array A has two or more consecutive non-zero + elements on the first sub-diagonal, so that there is + a block larger than 2-by-2 on the diagonal.""" + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 5: + if fact == 'F': + error_text = """the Schur factor S supplied in + the array A has a 2-by-2 diagonal block with real + eigenvalues instead of a complex conjugate pair.""" + e = ValueError(error_text) + e.info = out[-1] + raise e + if out[-1] == 6: + if fact == 'N': + error_text = """the LAPACK Library routine DGEES + has failed to converge. This failure is not likely + to occur. The matrix B will be unaltered but A will + be destroyed.""" + e = ValueError(error_text) + e.info = out[-1] + raise e + U,scale,wr,wi = out[:-1] + w = _np.zeros(n,'complex64') + w.real = wr[0:n] + w.imag = wi[0:n] + return U,scale,w + def sb04md(n,m,A,B,C,ldwork=None): """X = sb04md(n,m,A,B,C[,ldwork])