Skip to content

Commit fcf4396

Browse files
bnavigatorbnavigator
committed
mb03rd schur to block-diagonal transform
Based on PR python-control#73 by @repagh Moved from transform to math single wrapper for all jobx parameter valies docstring in numpydoc (python-control#100) run all jobx and sort parameter values in test check the returned complex eigenvalues in test
1 parent de9090a commit fcf4396

File tree

7 files changed

+350
-244
lines changed

7 files changed

+350
-244
lines changed

‎slycot/__init__.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,8 +24,8 @@
2424

2525
# Identification routines (0/5 wrapped)
2626

27-
# Mathematical routines (6/81 wrapped)
28-
from .math import mc01td, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
27+
# Mathematical routines (7/81 wrapped)
28+
from .math import mc01td, mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd
2929

3030
# Synthesis routines (14/50 wrapped)
3131
from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od

‎slycot/math.py‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,234 @@
2020
from . import _wrapper
2121
import warnings
2222

23+
import numpy as np
24+
25+
26+
def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0):
27+
"""Ar, Xr, blsize, W = mb03rd(n, A, [X, jobx, sort, pmax, tol])
28+
29+
To reduce a matrix `A` in real Schur form to a block-diagonal form
30+
using well-conditioned non-orthogonal similarity transformations.
31+
The condition numbers of the transformations used for reduction
32+
are roughly bounded by `pmax`*`pmax`, where `pmax` is a given value.
33+
The transformations are optionally postmultiplied in a given
34+
matrix `X`. The real Schur form is optionally ordered, so that
35+
clustered eigenvalues are grouped in the same block.
36+
37+
Parameters
38+
----------
39+
n : int
40+
The order of the matrices `A` and `X`. `n` >= 0.
41+
A : (n, n) array_like
42+
The matrix `A` to be block-diagonalized, in real Schur form.
43+
X : (n, n) array_like, optional
44+
A given matrix `X`, for accumulation of transformations (only if
45+
`jobx`='U')
46+
jobx : {'N', 'U'}, optional
47+
Specifies whether or not the transformations are
48+
accumulated, as follows:
49+
50+
:= 'N': The transformations are not accumulated
51+
:= 'U': The transformations are accumulated in `Xr` (default)
52+
53+
sort : {'N', 'S', 'C', 'B'}, optional
54+
Specifies whether or not the diagonal blocks of the real
55+
Schur form are reordered, as follows:
56+
57+
:= 'N': The diagonal blocks are not reordered (default);
58+
:= 'S': The diagonal blocks are reordered before each
59+
step of reduction, so that clustered eigenvalues
60+
appear in the same block;
61+
:= 'C': The diagonal blocks are not reordered, but the
62+
"closest-neighbour" strategy is used instead of
63+
the standard "closest to the mean" strategy
64+
(see Notes_);
65+
:= 'B': The diagonal blocks are reordered before each
66+
step of reduction, and the "closest-neighbour"
67+
strategy is used (see Notes_).
68+
69+
pmax : float, optional
70+
An upper bound for the infinity norm of elementary
71+
submatrices of the individual transformations used for
72+
reduction (see Notes_). `pmax` >= 1.0
73+
tol : float, optional
74+
The tolerance to be used in the ordering of the diagonal
75+
blocks of the real Schur form matrix.
76+
If the user sets `tol` > 0, then the given value of `tol` is
77+
used as an absolute tolerance: a block `i` and a temporarily
78+
fixed block 1 (the first block of the current trailing
79+
submatrix to be reduced) are considered to belong to the
80+
same cluster if their eigenvalues satisfy
81+
82+
.. math:: | \\lambda_1 - \\lambda_i | <= tol.
83+
84+
If the user sets `tol` < 0, then the given value of tol is
85+
used as a relative tolerance: a block i and a temporarily
86+
fixed block 1 are considered to belong to the same cluster
87+
if their eigenvalues satisfy, for ``j = 1, ..., n``
88+
89+
.. math:: | \\lambda_1 - \\lambda_i | <= | tol | * \\max | \\lambda_j |.
90+
91+
If the user sets `tol` = 0, then an implicitly computed,
92+
default tolerance, defined by ``tol = SQRT( SQRT( EPS ) )``
93+
is used instead, as a relative tolerance, where `EPS` is
94+
the machine precision (see LAPACK Library routine DLAMCH).
95+
If `sort` = 'N' or 'C', this parameter is not referenced.
96+
97+
Returns
98+
-------
99+
Ar : (n, n) ndarray
100+
Contains the computed block-diagonal matrix, in real Schur
101+
canonical form. The non-diagonal blocks are set to zero.
102+
Xr : (n, n) ndarray or None
103+
Contains the product of the given matrix `X` and the
104+
transformation matrix that reduced `A` to block-diagonal
105+
form. The transformation matrix is itself a product of
106+
non-orthogonal similarity transformations having elements
107+
with magnitude less than or equal to `pmax`.
108+
If `jobx` = 'N', this array is returned as None
109+
blsize : (n,) ndarray
110+
The orders of the resulting diagonal blocks of the matrix `Ar`.
111+
W : (n,) complex ndarray
112+
Contains the complex eigenvalues of the matrix `A`.
113+
114+
Notes
115+
-----
116+
**Method**
117+
118+
Consider first that `sort` = 'N'. Let
119+
120+
::
121+
122+
( A A )
123+
( 11 12 )
124+
A = ( ),
125+
( 0 A )
126+
( 22 )
127+
128+
be the given matrix in real Schur form, where initially :math:`A_{11}` is the
129+
first diagonal block of dimension 1-by-1 or 2-by-2. An attempt is
130+
made to compute a transformation matrix `X` of the form
131+
132+
::
133+
134+
( I P )
135+
X = ( ) (1)
136+
( 0 I )
137+
138+
(partitioned as `A`), so that
139+
140+
::
141+
142+
( A 0 )
143+
-1 ( 11 )
144+
X A X = ( ),
145+
( 0 A )
146+
( 22 )
147+
148+
and the elements of `P` do not exceed the value `pmax` in magnitude.
149+
An adaptation of the standard method for solving Sylvester
150+
equations [1]_, which controls the magnitude of the individual
151+
elements of the computed solution [2]_, is used to obtain matrix `P`.
152+
When this attempt failed, an 1-by-1 (or 2-by-2) diagonal block of
153+
:math:`A_{22}` , whose eigenvalue(s) is (are) the closest to the mean of those
154+
of :math:`A_{11}` is selected, and moved by orthogonal similarity
155+
transformations in the leading position of :math:`A_{22}` ; the moved diagonal
156+
block is then added to the block :math:`A_{11}` , increasing its order by 1
157+
(or 2). Another attempt is made to compute a suitable
158+
transformation matrix X with the new definitions of the blocks :math:`A_{11}`
159+
and :math:`A_{22}` . After a successful transformation matrix `X` has been
160+
obtained, it postmultiplies the current transformation matrix
161+
(if `jobx` = 'U'), and the whole procedure is repeated for the
162+
matrix :math:`A_{22}`.
163+
164+
When `sort` = 'S', the diagonal blocks of the real Schur form are
165+
reordered before each step of the reduction, so that each cluster
166+
of eigenvalues, defined as specified in the definition of TOL,
167+
appears in adjacent blocks. The blocks for each cluster are merged
168+
together, and the procedure described above is applied to the
169+
larger blocks. Using the option `sort` = 'S' will usually provide
170+
better efficiency than the standard option (`sort` = 'N'), proposed
171+
in [2]_, because there could be no or few unsuccessful attempts
172+
to compute individual transformation matrices `X` of the form (1).
173+
However, the resulting dimensions of the blocks are usually
174+
larger; this could make subsequent calculations less efficient.
175+
176+
When `sort` = 'C' or 'B', the procedure is similar to that for
177+
`sort` = 'N' or 'S', respectively, but the block of :math:`A_{22}` whose
178+
eigenvalue(s) is (are) the closest to those of :math:`A_{11}` (not to their
179+
mean) is selected and moved to the leading position of :math:`A_{22}`. This
180+
is called the "closest-neighbour" strategy.
181+
182+
**Numerical Aspects**
183+
184+
The algorithm usually requires :math:`\mathcal{O}(N^3)` operations,
185+
but :math:`\mathcal{O}(N^4)` are
186+
possible in the worst case, when all diagonal blocks in the real
187+
Schur form of `A` are 1-by-1, and the matrix cannot be diagonalized
188+
by well-conditioned transformations.
189+
190+
**Further Comments**
191+
192+
The individual non-orthogonal transformation matrices used in the
193+
reduction of `A` to a block-diagonal form have condition numbers
194+
of the order `pmax`*`pmax`. This does not guarantee that their product
195+
is well-conditioned enough. The routine can be easily modified to
196+
provide estimates for the condition numbers of the clusters of
197+
eigenvalues.
198+
199+
**Contributor**
200+
201+
V. Sima, Katholieke Univ. Leuven, Belgium, June 1998.
202+
Partly based on the RASP routine BDIAG by A. Varga, German
203+
Aerospace Center, DLR Oberpfaffenhofen.
204+
205+
**Revisions**
206+
207+
\V. Sima, Research Institute for Informatics, Bucharest, Apr. 2003.
208+
209+
References
210+
----------
211+
.. [1] Bartels, R.H. and Stewart, G.W. T
212+
Solution of the matrix equation A X + XB = C.
213+
Comm. A.C.M., 15, pp. 820-826, 1972.
214+
215+
.. [2] Bavely, C. and Stewart, G.W.
216+
An Algorithm for Computing Reducing Subspaces by Block
217+
Diagonalization.
218+
SIAM J. Numer. Anal., 16, pp. 359-367, 1979.
219+
220+
.. [3] Demmel, J.
221+
The Condition Number of Equivalence Transformations that
222+
Block Diagonalize Matrix Pencils.
223+
SIAM J. Numer. Anal., 20, pp. 599-610, 1983.
224+
225+
"""
226+
hidden = ' (hidden by the wrapper)'
227+
arg_list = ['jobx', 'sort', 'n', 'pmax',
228+
'A', 'lda' + hidden, 'X', 'ldx' + hidden,
229+
'nblcks', 'blsize', 'wr', 'wi', 'tol',
230+
'dwork' + hidden, 'info']
231+
232+
if X is None:
233+
X = np.zeros((1, n))
234+
235+
Ar, Xr, nblcks, blsize, wr, wi, info = _wrapper.mb03rd(
236+
jobx, sort, n, pmax, A, X, tol)
237+
238+
if info < 0:
239+
fmt = "The following argument had an illegal value: '{}'"
240+
e = ValueError(fmt.format(arg_list[-info - 1]))
241+
e.info = info
242+
raise e
243+
if jobx == 'N':
244+
Xr = None
245+
else:
246+
Xr = Xr[:n, :n]
247+
Ar = Ar[:n, :n]
248+
W = wr + 1J*wi
249+
return Ar, Xr, blsize[:nblcks], W
250+
23251

24252
def mb03vd(n, ilo, ihi, A):
25253
""" HQ, Tau = mb03vd(n, ilo, ihi, A)

‎slycot/src/math.pyf‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,24 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f
1212
integer intent(out) :: info
1313
end subroutine mc01td
1414

15+
subroutine mb03rd(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
16+
character intent(in) :: jobx
17+
character intent(in),required :: sort
18+
integer intent(in),required,check(n>=0) :: n
19+
double precision intent(in),required,check(pmax>=1.0) :: pmax
20+
double precision intent(in,out,copy),dimension(lda,n),depend(n) :: a
21+
integer intent(hide),check(lda>=max(1,n)) :: lda=shape(a,0)
22+
double precision intent(in,out,copy),dimension(ldx,n),depend(n) :: x
23+
integer intent(hide),check((*jobx == 'N' && ldx>=1) || (*jobx == 'U' && ldx >= max(1,n))) :: ldx=shape(x,0)
24+
integer intent(out) :: nblcks
25+
integer intent(out),dimension(n) :: blsize
26+
double precision intent(out),dimension(n) :: wr
27+
double precision intent(out),dimension(n) :: wi
28+
double precision intent(in) :: tol
29+
double precision intent(cache,hide),dimension(n) :: dwork
30+
integer intent(out) :: info
31+
end subroutine mb03rd
32+
1533
subroutine mb03vd(n,p,ilo,ihi,a,lda1,lda2,tau,ldtau,dwork,info) ! in MB03VD.f
1634
integer intent(in),check(n>=0) :: n
1735
integer intent(hide),depend(a),check(p>=1) :: p=shape(a,2)
@@ -97,5 +115,3 @@ subroutine mb05nd(n,delta,a,lda,ex,ldex,exint,ldexin,tol,iwork,dwork,ldwork,info
97115
integer intent(out) :: info
98116
end subroutine mb05nd
99117

100-
! This file was auto-generated with f2py (version:2).
101-
! See http://cens.ioc.ee/projects/f2py2e/

‎slycot/src/transform.pyf‎

Lines changed: 7 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ subroutine tb03ad_l(leri,equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,pcoeff,
4545
double precision :: tol = 0
4646
integer intent(hide,cache),dimension(n+max(m,p)) :: iwork
4747
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
48-
integer :: ldwork = max( 2*n + 3*max(m,p), p*(p+2))
48+
integer :: ldwork = max( 2*n + 3*max(m,p), p*(p+2))
4949
integer intent(out) :: info
5050
end subroutine tb03ad_l
5151
subroutine tb03ad_r(leri,equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,vcoeff,ldvco1,ldvco2,tol,iwork,dwork,ldwork,info) ! in :new:TB03AD.f
@@ -77,7 +77,7 @@ subroutine tb03ad_r(leri,equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,pcoeff,
7777
double precision :: tol = 0
7878
integer intent(hide,cache),dimension(n+max(m,p)) :: iwork
7979
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
80-
integer :: ldwork = max( 2*n + 3*max(m,p), m*(m+2))
80+
integer :: ldwork = max( 2*n + 3*max(m,p), m*(m+2))
8181
integer intent(out) :: info
8282
end subroutine tb03ad_r
8383
subroutine tb04ad_r(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,tol1,tol2,iwork,dwork,ldwork,info) ! in TB04AD.f
@@ -99,7 +99,7 @@ subroutine tb04ad_r(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddc
9999
double precision intent(out),dimension(max(1,p),n+1),depend(p,n) :: dcoeff
100100
integer intent(hide),depend(dcoeff) :: lddcoe=shape(dcoeff,0)
101101
double precision intent(out),dimension(max(1,p),max(1,m),n+1),depend(p,m,n) :: ucoeff
102-
integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
102+
integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
103103
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
104104
double precision :: tol1 = 0
105105
double precision :: tol2 = 0
@@ -284,7 +284,7 @@ subroutine tc04ad_l(leri,m,p,index_bn,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,
284284
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
285285
integer depend(m,p) :: ldwork = max(m,p)*(max(m,p)+4)
286286
integer intent(out) :: info
287-
end subroutine tc04ad_l
287+
end subroutine tc04ad_l
288288
subroutine tc04ad_r(leri,m,p,index_bn,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,n,rcond,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in TC04AD.f
289289
fortranname tc04ad
290290
character intent(hide) :: leri = 'R'
@@ -325,7 +325,7 @@ subroutine td04ad_r(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
325325
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
326326
integer intent(in,out) :: nr !=sum(index_bn)
327327
double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
328-
integer intent(hide),depend(a) :: lda = shape(a,0)
328+
integer intent(hide),depend(a) :: lda = shape(a,0)
329329
double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
330330
integer intent(hide),depend(b) :: ldb = shape(b,0)
331331
double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
@@ -351,7 +351,7 @@ subroutine td04ad_c(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
351351
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
352352
integer intent(in,out) :: nr != sum(index_bn)
353353
double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
354-
integer intent(hide),depend(a) :: lda = shape(a,0)
354+
integer intent(hide),depend(a) :: lda = shape(a,0)
355355
double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
356356
integer intent(hide),depend(b) :: ldb = shape(b,0)
357357
double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
@@ -525,43 +525,6 @@ subroutine tg01fd_uu(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ld
525525
double precision intent(in) :: tol
526526
integer intent(cache,hide),dimension(ldwork),depend(ldwork) :: iwork
527527
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
528-
integer required intent(in) :: ldwork
528+
integer required intent(in) :: ldwork
529529
integer intent(out) :: info
530530
end subroutine tg01fd_uu
531-
subroutine mb03rd_n(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
532-
fortranname mb03rd
533-
character intent(hide) :: jobx = 'N'
534-
character intent(in),required :: sort
535-
integer intent(in),required,check(n>0) :: n
536-
double precision intent(in),required,check(pmax>=1.0) :: pmax
537-
double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
538-
integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
539-
double precision intent(cache,hide) :: x
540-
integer intent(in,hide) :: ldx=1
541-
integer intent(out) :: nblcks
542-
integer intent(out),dimension(n),depend(n) :: blsize
543-
double precision intent(out),dimension(n),depend(n) :: wr
544-
double precision intent(out),dimension(n),depend(n) :: wi
545-
double precision intent(in) :: tol
546-
double precision intent(cache,hide),dimension(n),depend(n) :: dwork
547-
integer intent(out) :: info
548-
end subroutine mb03rd_n
549-
subroutine mb03rd_u(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
550-
fortranname mb03rd
551-
character intent(hide) :: jobx = 'U'
552-
character intent(in),required :: sort
553-
integer intent(in),required,check(n>0) :: n
554-
double precision intent(in),required,check(pmax>=1.0) :: pmax
555-
double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
556-
integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
557-
double precision intent(in,out,copy),dimension(n,n),depend(n) :: x
558-
integer intent(hide),depend(x) :: ldx=MAX(shape(x,0),1)
559-
integer intent(out) :: nblcks
560-
integer intent(out),dimension(n),depend(n) :: blsize
561-
double precision intent(out),dimension(n),depend(n) :: wr
562-
double precision intent(out),dimension(n),depend(n) :: wi
563-
double precision intent(in) :: tol
564-
double precision intent(cache,hide),dimension(n),depend(n) :: dwork
565-
integer intent(out) :: info
566-
end subroutine mb03rd_u
567-

0 commit comments

Comments
 (0)