Skip to content

Complex pyf #1

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
lytex merged 5 commits into from
Apr 3, 2020
Merged

Complex pyf #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f963190
clean __init__
bnavigator Apr 3, 2020
66b40b8
rework signatures for ab08nd and ab08nz
bnavigator Apr 3, 2020
7581eae
rework analysis.ab08nz() function
bnavigator Apr 3, 2020
6adae10
extend tests for ab08nX
bnavigator Apr 3, 2020
2fba484
fix docstring
bnavigator Apr 3, 2020
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
9 changes: 5 additions & 4 deletions slycot/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
66 changes: 40 additions & 26 deletions slycot/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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':
Expand Down
65 changes: 33 additions & 32 deletions slycot/src/analysis.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
105 changes: 65 additions & 40 deletions slycot/tests/test_ab08n.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,53 +5,78 @@
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 """
""" Test regular pencil construction ab08nX 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)
"Test Construct regular pencil for real matrices"
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)
"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)


if __name__ == "__main__":
Expand Down