Skip to content

Add mb03rd wrapper: schur to block-diagonal transform #116

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
merged 2 commits into from
May 3, 2020
Merged

Add mb03rd wrapper: schur to block-diagonal transform #116

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
de9090a
addition of wrapper for mb03rd
repagh Aug 2, 2019
fcf4396
mb03rd schur to block-diagonal transform
bnavigator May 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
4 changes: 2 additions & 2 deletions slycot/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@

# Identification routines (0/5 wrapped)

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

# Synthesis routines (14/50 wrapped)
from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od
Expand Down
228 changes: 228 additions & 0 deletions slycot/math.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,234 @@
from . import _wrapper
import warnings

import numpy as np


def mb03rd(n, A, X=None, jobx='U', sort='N', pmax=1.0, tol=0.0):
"""Ar, Xr, blsize, W = mb03rd(n, A, [X, jobx, sort, pmax, tol])

To reduce a matrix `A` in real Schur form to a block-diagonal form
using well-conditioned non-orthogonal similarity transformations.
The condition numbers of the transformations used for reduction
are roughly bounded by `pmax`*`pmax`, where `pmax` is a given value.
The transformations are optionally postmultiplied in a given
matrix `X`. The real Schur form is optionally ordered, so that
clustered eigenvalues are grouped in the same block.

Parameters
----------
n : int
The order of the matrices `A` and `X`. `n` >= 0.
A : (n, n) array_like
The matrix `A` to be block-diagonalized, in real Schur form.
X : (n, n) array_like, optional
A given matrix `X`, for accumulation of transformations (only if
`jobx`='U')
jobx : {'N', 'U'}, optional
Specifies whether or not the transformations are
accumulated, as follows:

:= 'N': The transformations are not accumulated
:= 'U': The transformations are accumulated in `Xr` (default)

sort : {'N', 'S', 'C', 'B'}, optional
Specifies whether or not the diagonal blocks of the real
Schur form are reordered, as follows:

:= 'N': The diagonal blocks are not reordered (default);
:= 'S': The diagonal blocks are reordered before each
step of reduction, so that clustered eigenvalues
appear in the same block;
:= 'C': The diagonal blocks are not reordered, but the
"closest-neighbour" strategy is used instead of
the standard "closest to the mean" strategy
(see Notes_);
:= 'B': The diagonal blocks are reordered before each
step of reduction, and the "closest-neighbour"
strategy is used (see Notes_).

pmax : float, optional
An upper bound for the infinity norm of elementary
submatrices of the individual transformations used for
reduction (see Notes_). `pmax` >= 1.0
tol : float, optional
The tolerance to be used in the ordering of the diagonal
blocks of the real Schur form matrix.
If the user sets `tol` > 0, then the given value of `tol` is
used as an absolute tolerance: a block `i` and a temporarily
fixed block 1 (the first block of the current trailing
submatrix to be reduced) are considered to belong to the
same cluster if their eigenvalues satisfy

.. math:: | \\lambda_1 - \\lambda_i | <= tol.

If the user sets `tol` < 0, then the given value of tol is
used as a relative tolerance: a block i and a temporarily
fixed block 1 are considered to belong to the same cluster
if their eigenvalues satisfy, for ``j = 1, ..., n``

.. math:: | \\lambda_1 - \\lambda_i | <= | tol | * \\max | \\lambda_j |.

If the user sets `tol` = 0, then an implicitly computed,
default tolerance, defined by ``tol = SQRT( SQRT( EPS ) )``
is used instead, as a relative tolerance, where `EPS` is
the machine precision (see LAPACK Library routine DLAMCH).
If `sort` = 'N' or 'C', this parameter is not referenced.

Returns
-------
Ar : (n, n) ndarray
Contains the computed block-diagonal matrix, in real Schur
canonical form. The non-diagonal blocks are set to zero.
Xr : (n, n) ndarray or None
Contains the product of the given matrix `X` and the
transformation matrix that reduced `A` to block-diagonal
form. The transformation matrix is itself a product of
non-orthogonal similarity transformations having elements
with magnitude less than or equal to `pmax`.
If `jobx` = 'N', this array is returned as None
blsize : (n,) ndarray
The orders of the resulting diagonal blocks of the matrix `Ar`.
W : (n,) complex ndarray
Contains the complex eigenvalues of the matrix `A`.

Notes
-----
**Method**

Consider first that `sort` = 'N'. Let

::

( A A )
( 11 12 )
A = ( ),
( 0 A )
( 22 )

be the given matrix in real Schur form, where initially :math:`A_{11}` is the
first diagonal block of dimension 1-by-1 or 2-by-2. An attempt is
made to compute a transformation matrix `X` of the form

::

( I P )
X = ( ) (1)
( 0 I )

(partitioned as `A`), so that

::

( A 0 )
-1 ( 11 )
X A X = ( ),
( 0 A )
( 22 )

and the elements of `P` do not exceed the value `pmax` in magnitude.
An adaptation of the standard method for solving Sylvester
equations [1]_, which controls the magnitude of the individual
elements of the computed solution [2]_, is used to obtain matrix `P`.
When this attempt failed, an 1-by-1 (or 2-by-2) diagonal block of
:math:`A_{22}` , whose eigenvalue(s) is (are) the closest to the mean of those
of :math:`A_{11}` is selected, and moved by orthogonal similarity
transformations in the leading position of :math:`A_{22}` ; the moved diagonal
block is then added to the block :math:`A_{11}` , increasing its order by 1
(or 2). Another attempt is made to compute a suitable
transformation matrix X with the new definitions of the blocks :math:`A_{11}`
and :math:`A_{22}` . After a successful transformation matrix `X` has been
obtained, it postmultiplies the current transformation matrix
(if `jobx` = 'U'), and the whole procedure is repeated for the
matrix :math:`A_{22}`.

When `sort` = 'S', the diagonal blocks of the real Schur form are
reordered before each step of the reduction, so that each cluster
of eigenvalues, defined as specified in the definition of TOL,
appears in adjacent blocks. The blocks for each cluster are merged
together, and the procedure described above is applied to the
larger blocks. Using the option `sort` = 'S' will usually provide
better efficiency than the standard option (`sort` = 'N'), proposed
in [2]_, because there could be no or few unsuccessful attempts
to compute individual transformation matrices `X` of the form (1).
However, the resulting dimensions of the blocks are usually
larger; this could make subsequent calculations less efficient.

When `sort` = 'C' or 'B', the procedure is similar to that for
`sort` = 'N' or 'S', respectively, but the block of :math:`A_{22}` whose
eigenvalue(s) is (are) the closest to those of :math:`A_{11}` (not to their
mean) is selected and moved to the leading position of :math:`A_{22}`. This
is called the "closest-neighbour" strategy.

**Numerical Aspects**

The algorithm usually requires :math:`\mathcal{O}(N^3)` operations,
but :math:`\mathcal{O}(N^4)` are
possible in the worst case, when all diagonal blocks in the real
Schur form of `A` are 1-by-1, and the matrix cannot be diagonalized
by well-conditioned transformations.

**Further Comments**

The individual non-orthogonal transformation matrices used in the
reduction of `A` to a block-diagonal form have condition numbers
of the order `pmax`*`pmax`. This does not guarantee that their product
is well-conditioned enough. The routine can be easily modified to
provide estimates for the condition numbers of the clusters of
eigenvalues.

**Contributor**

V. Sima, Katholieke Univ. Leuven, Belgium, June 1998.
Partly based on the RASP routine BDIAG by A. Varga, German
Aerospace Center, DLR Oberpfaffenhofen.

**Revisions**

\V. Sima, Research Institute for Informatics, Bucharest, Apr. 2003.

References
----------
.. [1] Bartels, R.H. and Stewart, G.W. T
Solution of the matrix equation A X + XB = C.
Comm. A.C.M., 15, pp. 820-826, 1972.

.. [2] Bavely, C. and Stewart, G.W.
An Algorithm for Computing Reducing Subspaces by Block
Diagonalization.
SIAM J. Numer. Anal., 16, pp. 359-367, 1979.

.. [3] Demmel, J.
The Condition Number of Equivalence Transformations that
Block Diagonalize Matrix Pencils.
SIAM J. Numer. Anal., 20, pp. 599-610, 1983.

"""
hidden = ' (hidden by the wrapper)'
arg_list = ['jobx', 'sort', 'n', 'pmax',
'A', 'lda' + hidden, 'X', 'ldx' + hidden,
'nblcks', 'blsize', 'wr', 'wi', 'tol',
'dwork' + hidden, 'info']

if X is None:
X = np.zeros((1, n))

Ar, Xr, nblcks, blsize, wr, wi, info = _wrapper.mb03rd(
jobx, sort, n, pmax, A, X, tol)

if info < 0:
fmt = "The following argument had an illegal value: '{}'"
e = ValueError(fmt.format(arg_list[-info - 1]))
e.info = info
raise e
if jobx == 'N':
Xr = None
else:
Xr = Xr[:n, :n]
Ar = Ar[:n, :n]
W = wr + 1J*wi
return Ar, Xr, blsize[:nblcks], W


def mb03vd(n, ilo, ihi, A):
""" HQ, Tau = mb03vd(n, ilo, ihi, A)
Expand Down
20 changes: 18 additions & 2 deletions slycot/src/math.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,24 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f
integer intent(out) :: info
end subroutine mc01td

subroutine mb03rd(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
character intent(in) :: jobx
character intent(in),required :: sort
integer intent(in),required,check(n>=0) :: n
double precision intent(in),required,check(pmax>=1.0) :: pmax
double precision intent(in,out,copy),dimension(lda,n),depend(n) :: a
integer intent(hide),check(lda>=max(1,n)) :: lda=shape(a,0)
double precision intent(in,out,copy),dimension(ldx,n),depend(n) :: x
integer intent(hide),check((*jobx == 'N' && ldx>=1) || (*jobx == 'U' && ldx >= max(1,n))) :: ldx=shape(x,0)
integer intent(out) :: nblcks
integer intent(out),dimension(n) :: blsize
double precision intent(out),dimension(n) :: wr
double precision intent(out),dimension(n) :: wi
double precision intent(in) :: tol
double precision intent(cache,hide),dimension(n) :: dwork
integer intent(out) :: info
end subroutine mb03rd

subroutine mb03vd(n,p,ilo,ihi,a,lda1,lda2,tau,ldtau,dwork,info) ! in MB03VD.f
integer intent(in),check(n>=0) :: n
integer intent(hide),depend(a),check(p>=1) :: p=shape(a,2)
Expand Down Expand Up @@ -97,5 +115,3 @@ subroutine mb05nd(n,delta,a,lda,ex,ldex,exint,ldexin,tol,iwork,dwork,ldwork,info
integer intent(out) :: info
end subroutine mb05nd

! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/
14 changes: 7 additions & 7 deletions slycot/src/transform.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ subroutine tb03ad_l(leri,equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,pcoeff,
double precision :: tol = 0
integer intent(hide,cache),dimension(n+max(m,p)) :: iwork
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
integer :: ldwork = max( 2*n + 3*max(m,p), p*(p+2))
integer :: ldwork = max( 2*n + 3*max(m,p), p*(p+2))
integer intent(out) :: info
end subroutine tb03ad_l
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
Expand Down Expand Up @@ -77,7 +77,7 @@ subroutine tb03ad_r(leri,equil,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,pcoeff,
double precision :: tol = 0
integer intent(hide,cache),dimension(n+max(m,p)) :: iwork
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
integer :: ldwork = max( 2*n + 3*max(m,p), m*(m+2))
integer :: ldwork = max( 2*n + 3*max(m,p), m*(m+2))
integer intent(out) :: info
end subroutine tb03ad_r
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
Expand All @@ -99,7 +99,7 @@ subroutine tb04ad_r(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddc
double precision intent(out),dimension(max(1,p),n+1),depend(p,n) :: dcoeff
integer intent(hide),depend(dcoeff) :: lddcoe=shape(dcoeff,0)
double precision intent(out),dimension(max(1,p),max(1,m),n+1),depend(p,m,n) :: ucoeff
integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
double precision :: tol1 = 0
double precision :: tol2 = 0
Expand Down Expand Up @@ -284,7 +284,7 @@ subroutine tc04ad_l(leri,m,p,index_bn,pcoeff,ldpco1,ldpco2,qcoeff,ldqco1,ldqco2,
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
integer depend(m,p) :: ldwork = max(m,p)*(max(m,p)+4)
integer intent(out) :: info
end subroutine tc04ad_l
end subroutine tc04ad_l
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
fortranname tc04ad
character intent(hide) :: leri = 'R'
Expand Down Expand Up @@ -325,7 +325,7 @@ subroutine td04ad_r(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
integer intent(in,out) :: nr !=sum(index_bn)
double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
integer intent(hide),depend(a) :: lda = shape(a,0)
integer intent(hide),depend(a) :: lda = shape(a,0)
double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
integer intent(hide),depend(b) :: ldb = shape(b,0)
double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
Expand All @@ -351,7 +351,7 @@ subroutine td04ad_c(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
integer intent(in,out) :: nr != sum(index_bn)
double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
integer intent(hide),depend(a) :: lda = shape(a,0)
integer intent(hide),depend(a) :: lda = shape(a,0)
double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
integer intent(hide),depend(b) :: ldb = shape(b,0)
double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
Expand Down Expand Up @@ -525,6 +525,6 @@ subroutine tg01fd_uu(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ld
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 required intent(in) :: ldwork
integer intent(out) :: info
end subroutine tg01fd_uu
Loading