1import numpy as np
2import ctypes as ct
3
4__all__ = ["ab09ad","ab09ax","tb01id","tb01wd"]
5
6def c_arr(w=False):
7    flags = ['F']
8    if w:
9        flags.append('W')
10    return np.ctypeslib.ndpointer(dtype=ct.c_double, ndim=1, flags=flags)
11def c_mat(w=False):
12    flags = ['F']
13    if w:
14        flags.append('W')
15    return np.ctypeslib.ndpointer(dtype=ct.c_double, ndim=2, flags=flags)
16
17slicot = ct.cdll.LoadLibrary("libslicot.so.0")
18
19################################################################
20## AB09AD
21#
22slicot.ab09ad_.restype = None
23slicot.ab09ad_.argtypes = [
24    ct.c_char_p, #DICO,
25    ct.c_char_p, #JOB,
26    ct.c_char_p, #EQUIL,
27    ct.c_char_p, #ORDSEL,
28    ct.POINTER(ct.c_int), #N,
29    ct.POINTER(ct.c_int), #M,
30    ct.POINTER(ct.c_int), #P,
31    ct.POINTER(ct.c_int), #NR,
32    c_mat(True), #A,
33    ct.POINTER(ct.c_int), #LDA,
34    c_mat(True), #B,
35    ct.POINTER(ct.c_int), #LDB,
36    c_mat(True), #C,
37    ct.POINTER(ct.c_int), #LDC,
38    c_arr(True), #HSV,
39    ct.POINTER(ct.c_double), #TOL,
40    np.ctypeslib.ndpointer(dtype=ct.c_int, ndim=1, flags=["F","W"]), #IWORK,
41    c_arr(True), #DWORK,
42    ct.POINTER(ct.c_int), #LDWORK,
43    ct.POINTER(ct.c_int), #IWARN,
44    ct.POINTER(ct.c_int), #INFO
45    ]
46
47def ab09ad(dico, job, equil, A, B, C, nr=None, tol=0):
48    ordsel = 'A' if nr is None else 'F'
49    Ar = np.array(A, dtype=np.float64, order='F', copy=True)
50    Br = np.array(B, dtype=np.float64, order='F', copy=True)
51    Cr = np.array(C, dtype=np.float64, order='F', copy=True)
52    _nr = ct.c_int(nr or 0)
53    hsv = np.zeros(Ar.shape[1], dtype=np.float64, order='F')
54    iwork = np.zeros(Ar.shape[1], dtype=np.int32, order='F')
55    dwork = np.zeros(Ar.size * 5, dtype=np.float64, order='F')
56    iwarn = ct.c_int()
57    info = ct.c_int()
58    n = ct.c_int(Ar.shape[0])
59    m = ct.c_int(Br.shape[1])
60    p = ct.c_int(Cr.shape[0])
61    lda = ct.c_int(Ar.shape[0])
62    ldb = ct.c_int(Br.shape[0])
63    ldc = ct.c_int(Cr.shape[0])
64    _tol = ct.c_double(tol)
65    ldwork = ct.c_int(dwork.shape[0])
66    slicot.ab09ad_(
67        dico, job, equil, ordsel, ct.byref(n), ct.byref(m), ct.byref(p),
68        ct.byref(_nr), Ar, ct.byref(lda), Br, ct.byref(ldb), Cr, ct.byref(ldc),
69        hsv, ct.byref(_tol), iwork, dwork, ct.byref(ldwork), ct.byref(iwarn), ct.byref(info))
70    if iwarn.value != 0:
71        print "ab09ad:", (
72            "The selected order NR is greater than the order of a minimal\n"
73            "realization of the given system. The order has been set according\n"
74            "to the minimal realization of the system.")
75    if info.value != 0:
76        if info.value < 0:
77            raise ValueError("%dth input argument is illegal" % -info.value())
78        if info.value == 1:
79            raise ValueError("the reduction of A to the real Schur form failed")
80        if info.value == 2:
81            raise ValueError("the state matrix A is not stable (if DICO = 'C') or not convergent (if DICO = 'D')")
82        if info.value == 3:
83            raise ValueError("the computation of Hankel singular values failed.")
84        raise ValueError("unknown error value %d" % info.value)
85    nr = _nr.value
86    if nr != Ar.shape[1]:
87        Ar = Ar[:nr,:nr]
88        Br = Br[:nr]
89        Cr = Cr[:,:nr]
90    return Ar, Br, Cr, hsv
91
92################################################################
93## AB09AX
94##
95slicot.ab09ax_.restype = None
96slicot.ab09ax_.argtypes = [
97    ct.c_char_p, #DICO,
98    ct.c_char_p, #JOB,
99    ct.c_char_p, #ORDSEL,
100    ct.POINTER(ct.c_int), #N,
101    ct.POINTER(ct.c_int), #M,
102    ct.POINTER(ct.c_int), #P,
103    ct.POINTER(ct.c_int), #NR,
104    c_mat(True), #A,
105    ct.POINTER(ct.c_int), #LDA,
106    c_mat(True), #B,
107    ct.POINTER(ct.c_int), #LDB,
108    c_mat(True), #C,
109    ct.POINTER(ct.c_int), #LDC,
110    c_arr(True), #HSV,
111    c_mat(True), #T,
112    ct.POINTER(ct.c_int), #LDT,
113    c_mat(True), #TI,
114    ct.POINTER(ct.c_int), #LDTI,
115    ct.POINTER(ct.c_double), #TOL,
116    np.ctypeslib.ndpointer(dtype=ct.c_int, ndim=1, flags=["F","W"]), #IWORK,
117    c_arr(True), #DWORK,
118    ct.POINTER(ct.c_int), #LDWORK,
119    ct.POINTER(ct.c_int), #IWARN,
120    ct.POINTER(ct.c_int), #INFO
121    ]
122
123def ab09ax(dico, job, A, B, C, nr=None, tol=0):
124    ordsel = 'A' if nr is None else 'F'
125    Ar = np.array(A, dtype=np.float64, order='F', copy=True)
126    Br = np.array(B, dtype=np.float64, order='F', copy=True)
127    Cr = np.array(C, dtype=np.float64, order='F', copy=True)
128    _nr = ct.c_int(nr or 0)
129    hsv = np.zeros(Ar.shape[1], dtype=np.float64, order='F')
130    T = np.zeros_like(Ar)
131    ldt = ct.c_int(Ar.shape[0])
132    Ti = np.zeros_like(Ar)
133    ldti = ct.c_int(Ar.shape[0])
134    iwork = np.zeros(Ar.shape[1], dtype=np.int32, order='F')
135    n = Ar.shape[0]
136    m = Br.shape[1]
137    p = Cr.shape[0]
138    minsz = max(1,n*(max(n,m,p)+5) + n*(n+1)/2)
139    dwork = np.zeros(minsz * 2, dtype=np.float64, order='F')
140    iwarn = ct.c_int()
141    info = ct.c_int()
142    n = ct.c_int(Ar.shape[0])
143    m = ct.c_int(Br.shape[1])
144    p = ct.c_int(Cr.shape[0])
145    lda = ct.c_int(Ar.shape[0])
146    ldb = ct.c_int(Br.shape[0])
147    ldc = ct.c_int(Cr.shape[0])
148    _tol = ct.c_double(tol)
149    ldwork = ct.c_int(dwork.shape[0])
150    slicot.ab09ax_(
151        dico, job, ordsel, ct.byref(n), ct.byref(m), ct.byref(p),
152        ct.byref(_nr), Ar, ct.byref(lda), Br, ct.byref(ldb), Cr, ct.byref(ldc),
153        hsv, T, ct.byref(ldt), Ti, ct.byref(ldti), ct.byref(_tol),
154        iwork, dwork, ct.byref(ldwork), ct.byref(iwarn), ct.byref(info))
155    if iwarn.value != 0:
156        print "ab09ad:", (
157            "The selected order NR is greater than the order of a minimal\n"
158            "realization of the given system. The order has been set according\n"
159            "to the minimal realization of the system.")
160    if info.value != 0:
161        if info.value < 0:
162            raise ValueError("%dth input argument is illegal" % -info.value())
163        if info.value == 1:
164            raise ValueError("the state matrix A is not stable (if DICO = 'C') or not convergent (if DICO = 'D')")
165        if info.value == 2:
166            raise ValueError("the computation of Hankel singular values failed.")
167        raise ValueError("unknown error value %d" % info.value)
168    nr = _nr.value
169    if nr != Ar.shape[1]:
170        Ar = Ar[:nr,:nr]
171        Br = Br[:nr]
172        Cr = Cr[:,:nr]
173        T = T[:,:nr]
174        Ti = Ti[:nr,:]
175    return T, Ti, Ar, Br, Cr, hsv
176
177################################################################
178## TB01ID
179##
180slicot.tb01id_.restype = None
181slicot.tb01id_.argtypes = [
182    ct.c_char_p, #JOB,
183    ct.POINTER(ct.c_int), #N,
184    ct.POINTER(ct.c_int), #M,
185    ct.POINTER(ct.c_int), #P,
186    ct.POINTER(ct.c_double), #MAXRED,
187    c_mat(True), #A,
188    ct.POINTER(ct.c_int), #LDA,
189    c_mat(True), #B,
190    ct.POINTER(ct.c_int), #LDB,
191    c_mat(True), #C,
192    ct.POINTER(ct.c_int), #LDC,
193    c_arr(True), #SCALE,
194    ct.POINTER(ct.c_int), #INFO
195    ]
196
197def tb01id(job, A, B, C, maxred=0):
198    Ar = np.array(A, dtype=np.float64, order='F', copy=True)
199    Br = np.array(B, dtype=np.float64, order='F', copy=True)
200    Cr = np.array(C, dtype=np.float64, order='F', copy=True)
201    n = ct.c_int(Ar.shape[0])
202    m = ct.c_int(Br.shape[1])
203    p = ct.c_int(Cr.shape[0])
204    lda = ct.c_int(Ar.shape[0])
205    ldb = ct.c_int(Br.shape[0])
206    ldc = ct.c_int(Cr.shape[0])
207    scale = np.zeros(Ar.shape[1], dtype=np.float64, order='F')
208    info = ct.c_int()
209    maxred_ = ct.c_double(maxred)
210    slicot.tb01id_(job, ct.byref(n), ct.byref(m), ct.byref(p), ct.byref(maxred_),
211                   Ar, ct.byref(lda), Br, ct.byref(ldb), Cr, ct.byref(ldc),
212                   scale, ct.byref(info))
213    if info.value != 0:
214        if info.value < 0:
215            raise ValueError("%dth input argument is illegal" % -info.value())
216        raise ValueError("unknown error value %d" % info.value)
217    return Ar, Br, Cr, maxred_.value, scale
218
219################################################################
220## TB01WD
221##
222slicot.tb01wd_.restype = None
223slicot.tb01wd_.argtypes = [
224    ct.POINTER(ct.c_int), #N,
225    ct.POINTER(ct.c_int), #M,
226    ct.POINTER(ct.c_int), #P,
227    c_mat(True), #A,
228    ct.POINTER(ct.c_int), #LDA,
229    c_mat(True), #B,
230    ct.POINTER(ct.c_int), #LDB,
231    c_mat(True), #C,
232    ct.POINTER(ct.c_int), #LDC,
233    c_mat(True), #U,
234    ct.POINTER(ct.c_int), #LDU,
235    c_arr(True), #WR,
236    c_arr(True), #WI,
237    c_arr(True), #DWORK,
238    ct.POINTER(ct.c_int), #LDWORK,
239    ct.POINTER(ct.c_int), #INFO
240    ]
241
242def tb01wd(A, B, C):
243    Ar = np.array(A, dtype=np.float64, order='F', copy=True)
244    Br = np.array(B, dtype=np.float64, order='F', copy=True)
245    Cr = np.array(C, dtype=np.float64, order='F', copy=True)
246    n = ct.c_int(Ar.shape[0])
247    m = ct.c_int(Br.shape[1])
248    p = ct.c_int(Cr.shape[0])
249    lda = ct.c_int(Ar.shape[0])
250    ldb = ct.c_int(Br.shape[0])
251    ldc = ct.c_int(Cr.shape[0])
252    U = np.zeros_like(Ar)
253    ldu = ct.c_int(Ar.shape[0])
254    WR = np.zeros(Ar.shape[0], dtype=np.float64, order='F')
255    WI = np.zeros(Ar.shape[0], dtype=np.float64, order='F')
256    dwork = np.zeros(Ar.shape[0] * 5, dtype=np.float64, order='F')
257    ldwork = ct.c_int(dwork.shape[0])
258    info = ct.c_int()
259    slicot.tb01wd_(
260        ct.byref(n), ct.byref(m), ct.byref(p),
261        Ar, ct.byref(lda), Br, ct.byref(ldb), Cr, ct.byref(ldc), U, ct.byref(ldu),
262        WR, WI, dwork, ct.byref(ldwork), ct.byref(info))
263    if info.value != 0:
264        if info.value < 0:
265            raise ValueError("%dth input argument is illegal" % -info.value())
266        raise ValueError("not all eigenvalues computed (%d not computed)" % info.value)
267    return Ar, Br, Cr, U, WR, WI
268