Skip to content

Commit 0b4022a

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 0b4022a

File tree

7 files changed

+351
-244
lines changed

7 files changed

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

24253
def mb03vd(n, ilo, ihi, A):
25254
""" 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)