Skip to content

add support for ag08bd, sb10jd and tg01fd #34

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
murrayrm merged 10 commits into from
Jul 29, 2018
Merged

add support for ag08bd, sb10jd and tg01fd #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c80fad0
add support for ag08bd, sb10jd and tg01fd
marcuslil Jul 17, 2018
388f8ee
correct ldwork for ag08bd
marcuslil Jul 17, 2018
3edd2c9
sb10jd, tg01fd
marcuslil Jul 17, 2018
34b3369
add tests for tg01fd
marcuslil Jul 18, 2018
4dc5451
test_tg01fd.py: make python2.7 compatibale and add a seed before crea…
marcuslil Jul 19, 2018
fb04bd6
correct documentation and spaces
marcuslil Jul 19, 2018
1b2c5b1
add ag08bd tests
marcuslil Jul 19, 2018
b804b92
add test for sb10jd
marcuslil Jul 19, 2018
40ccf91
test_sb10jd: round some low numbers
marcuslil Jul 19, 2018
dcbee80
change ld*=shape(*,0) to ld*=MAX(shape(*,0),1) for ag08bd, sb10jd and…
marcuslil Jul 22, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 113 additions & 0 deletions slycot/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
33 changes: 33 additions & 0 deletions slycot/src/analysis.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
19 changes: 19 additions & 0 deletions slycot/src/synthesis.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
88 changes: 87 additions & 1 deletion slycot/src/transform.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
79 changes: 79 additions & 0 deletions slycot/synthesis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 =
Expand Down
Loading