From de9090ab28aba913b77d93ea75d6cafe62abf349 Mon Sep 17 00:00:00 2001 From: Rene van Paassen Date: Fri, 2 Aug 2019 22:53:42 +0200 Subject: [PATCH 1/2] addition of wrapper for mb03rd --- slycot/src/transform.pyf | 37 +++++++++++++ slycot/tests/test_mb03rd.py | 63 +++++++++++++++++++++ slycot/transform.py | 106 ++++++++++++++++++++++++++++++++++++ 3 files changed, 206 insertions(+) create mode 100644 slycot/tests/test_mb03rd.py diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf index dc7277b1..19d6729c 100644 --- a/slycot/src/transform.pyf +++ b/slycot/src/transform.pyf @@ -528,3 +528,40 @@ subroutine tg01fd_uu(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ld integer required intent(in) :: ldwork integer intent(out) :: info end subroutine tg01fd_uu +subroutine mb03rd_n(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f + fortranname mb03rd + character intent(hide) :: jobx = 'N' + 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(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(cache,hide) :: x + integer intent(in,hide) :: ldx=1 + integer intent(out) :: nblcks + integer intent(out),dimension(n),depend(n) :: blsize + double precision intent(out),dimension(n),depend(n) :: wr + double precision intent(out),dimension(n),depend(n) :: wi + double precision intent(in) :: tol + double precision intent(cache,hide),dimension(n),depend(n) :: dwork + integer intent(out) :: info +end subroutine mb03rd_n +subroutine mb03rd_u(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f + fortranname mb03rd + character intent(hide) :: jobx = 'U' + 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(n,n),depend(n) :: a + integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) + double precision intent(in,out,copy),dimension(n,n),depend(n) :: x + integer intent(hide),depend(x) :: ldx=MAX(shape(x,0),1) + integer intent(out) :: nblcks + integer intent(out),dimension(n),depend(n) :: blsize + double precision intent(out),dimension(n),depend(n) :: wr + double precision intent(out),dimension(n),depend(n) :: wi + double precision intent(in) :: tol + double precision intent(cache,hide),dimension(n),depend(n) :: dwork + integer intent(out) :: info +end subroutine mb03rd_u + diff --git a/slycot/tests/test_mb03rd.py b/slycot/tests/test_mb03rd.py new file mode 100644 index 00000000..6c820871 --- /dev/null +++ b/slycot/tests/test_mb03rd.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python +# +# test_mb03rd.py - test suite for Shur form reduction +# RvP, 31 Jul 2019 +import unittest +from slycot import transform +import numpy as np +from numpy.testing import assert_raises, assert_almost_equal, assert_equal +from scipy.linalg import schur + +test1_A = np.array([ + [ 1., -1., 1., 2., 3., 1., 2., 3.], + [ 1., 1., 3., 4., 2., 3., 4., 2.], + [ 0., 0., 1., -1., 1., 5., 4., 1.], + [ 0., 0., 0., 1., -1., 3., 1., 2.], + [ 0., 0., 0., 1., 1., 2., 3., -1.], + [ 0., 0., 0., 0., 0., 1., 5., 1.], + [ 0., 0., 0., 0., 0., 0., 0.99999999, -0.99999999 ], + [ 0., 0., 0., 0., 0., 0., 0.99999999, 0.99999999 ] + ]) +test1_n = test1_A.shape[0] + +test1_Ar = np.array([ + [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000 ], + [ 1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ], + ]) + +test1_Xr = np.array([ + [ 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957 ], + [ 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755 ], + [ 0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148 ], + [ 0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534 ], + [ 0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000 ] + ]) + +test1_pmax = 1e3 +test1_tol = 0.01 +class test_mb03rd(unittest.TestCase): + def test1(self): + # create schur form with scipy + A, X = schur(test1_A) + Ah, Xh = np.copy(A), np.copy(X) + # on this basis, get the transform + Ar, Xr, blks, eig = transform.mb03rd( + test1_n, A, X, 'U', 'S', test1_pmax, test1_tol) + # ensure X and A are unchanged + assert_equal(A, Ah) + assert_equal(X, Xh) + # compare to test case results + assert_almost_equal(Ar, test1_Ar, decimal=4) + assert_almost_equal(Xr, test1_Xr, decimal=4) + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/transform.py b/slycot/transform.py index b9a72fa3..52357d72 100644 --- a/slycot/transform.py +++ b/slycot/transform.py @@ -1350,4 +1350,110 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld return A,E,B,C,ranke,rnka22,Q,Z +def mb03rd(n,A,X=None,jobx='U',sort='N',pmax=1.0,tol=0.0): + """ A,X,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='U' + A,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='N' + + 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. + + Required arguments: + n : input int + The order of the matrices A and X. n >= 0. + A : input rank-2 array('d') with bounds (n,n) + the matrix A to be block-diagonalized, in real Schur form. + Optional arguments: + X : input rank-2 array('d') with bounds (n,n) + a given matrix X, for accumulation of transformations (only if + jobx='U' + jobx : input char*1 + Specifies whether or not the transformations are + accumulated, as follows: + = 'N': The transformations are not accumulated; + = 'U': The transformations are accumulated in X (the + given matrix X is updated). + sort : input char*1 + Specifies whether or not the diagonal blocks of the real + Schur form are reordered, as follows: + = 'N': The diagonal blocks are not reordered; + = '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 METHOD); + = 'B': The diagonal blocks are reordered before each + step of reduction, and the "closest-neighbour" + strategy is used (see METHOD). + pmax : input float + An upper bound for the infinity norm of elementary + submatrices of the individual transformations used for + reduction (see METHOD). PMAX >= 1.0D0. + tol : input float + 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 + + | 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, + + | 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. + Return objects: + Ar : output rank-2 array('d') with bounds (n,n) + Contains the computed block-diagonal matrix, in real Schur + canonical form. The non-diagonal blocks are set to zero. + Xr : output rank-2 array('d') with bounds (n,n) + 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 not referenced, and not returned + blksize : output rank-1 array('i') with bounds (n) + The orders of the resulting diagonal blocks of the matrix Ar. + W : output rank-1 array('c') size (n) + This arrays contain the eigenvalues of the matrix A. +""" + hidden = ' (hidden by the wrapper)' + arg_list = ('jobx', 'sort', 'n', 'pmax', 'A', 'LDA'+hidden, + 'X', 'LDX'+hidden, 'nblks'+hidden, 'blsize'+hidden, + 'WR'+hidden, 'WI'+hidden, 'tol', + 'DWORK'+hidden, 'INFO'+hidden) + if jobx == 'N': + out = _wrapper.mb03rd_n(sort, n, pmax, A, tol) + else: + if X is None: + X = _np.eye(n) + out = _wrapper.mb03rd_u(sort, n, pmax, A, X, tol) + + if out[-1] < 0: + error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] + e = ValueError(error_text) + e.info = out[-1] + raise e + if jobx == 'N': + return out[0], out[2][:out[1]], out[-3] + out[-2]*1j + else: + return out[0], out[1], out[3][:out[2]], out[-3] + out[-2]*1j + # to be replaced by python wrappers From fcf4396a5cbe82ceb332cffd95a0720ceedab9ce Mon Sep 17 00:00:00 2001 From: bnavigator Date: Sun, 3 May 2020 18:15:14 +0200 Subject: [PATCH 2/2] mb03rd schur to block-diagonal transform Based on PR #73 by @repagh Moved from transform to math single wrapper for all jobx parameter valies docstring in numpydoc (#100) run all jobx and sort parameter values in test check the returned complex eigenvalues in test --- slycot/__init__.py | 4 +- slycot/math.py | 228 ++++++++++++++++++++++++++++++++++++ slycot/src/math.pyf | 20 +++- slycot/src/transform.pyf | 51 ++------ slycot/tests/test_mb.py | 70 ++++++++++- slycot/tests/test_mb03rd.py | 63 ---------- slycot/transform.py | 158 ++++--------------------- 7 files changed, 350 insertions(+), 244 deletions(-) delete mode 100644 slycot/tests/test_mb03rd.py diff --git a/slycot/__init__.py b/slycot/__init__.py index 715a6214..981f43ac 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -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 diff --git a/slycot/math.py b/slycot/math.py index c029b6cc..94a9e98b 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -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) diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index ecf94f2c..8cc06ace 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -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) @@ -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/ diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf index 19d6729c..524996aa 100644 --- a/slycot/src/transform.pyf +++ b/slycot/src/transform.pyf @@ -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 @@ -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 @@ -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 @@ -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' @@ -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 @@ -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 @@ -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 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 -subroutine mb03rd_n(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f - fortranname mb03rd - character intent(hide) :: jobx = 'N' - 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(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) - double precision intent(cache,hide) :: x - integer intent(in,hide) :: ldx=1 - integer intent(out) :: nblcks - integer intent(out),dimension(n),depend(n) :: blsize - double precision intent(out),dimension(n),depend(n) :: wr - double precision intent(out),dimension(n),depend(n) :: wi - double precision intent(in) :: tol - double precision intent(cache,hide),dimension(n),depend(n) :: dwork - integer intent(out) :: info -end subroutine mb03rd_n -subroutine mb03rd_u(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f - fortranname mb03rd - character intent(hide) :: jobx = 'U' - 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(n,n),depend(n) :: a - integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1) - double precision intent(in,out,copy),dimension(n,n),depend(n) :: x - integer intent(hide),depend(x) :: ldx=MAX(shape(x,0),1) - integer intent(out) :: nblcks - integer intent(out),dimension(n),depend(n) :: blsize - double precision intent(out),dimension(n),depend(n) :: wr - double precision intent(out),dimension(n),depend(n) :: wi - double precision intent(in) :: tol - double precision intent(cache,hide),dimension(n),depend(n) :: dwork - integer intent(out) :: info -end subroutine mb03rd_u - diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index 7dc81d6a..4ea494a6 100755 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -5,14 +5,82 @@ import unittest import numpy as np +from scipy.linalg import schur -from slycot import mb03vd, mb03vy, mb03wd, mb05md, mb05nd +from slycot import mb03rd, mb03vd, mb03vy, mb03wd, mb05md, mb05nd from numpy.testing import assert_allclose class test_mb(unittest.TestCase): + def test_mb03rd(self): + """ Test for Schur form reduction. + + RvP, 31 Jul 2019""" + + test1_A = np.array([ + [ 1., -1., 1., 2., 3., 1., 2., 3.], + [ 1., 1., 3., 4., 2., 3., 4., 2.], + [ 0., 0., 1., -1., 1., 5., 4., 1.], + [ 0., 0., 0., 1., -1., 3., 1., 2.], + [ 0., 0., 0., 1., 1., 2., 3., -1.], + [ 0., 0., 0., 0., 0., 1., 5., 1.], + [ 0., 0., 0., 0., 0., 0., 0.99999999, -0.99999999 ], + [ 0., 0., 0., 0., 0., 0., 0.99999999, 0.99999999 ] + ]) + test1_n = test1_A.shape[0] + + test1_Ar = np.array([ + [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000 ], + [ 1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ], + ]) + + test1_Xr = np.array([ + [ 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957 ], + [ 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755 ], + [ 0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148 ], + [ 0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534 ], + [ 0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000 ], + [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000 ] + ]) + + test1_W = np.array([1+1j, 1-1j, + 1+1j, 1-1j, + 0.99999+0.99999j, 0.99999-0.99999j, + 1., 1.]) + + test1_pmax = 1e3 + test1_tol = 0.01 + # create schur form with scipy + A, X = schur(test1_A) + Ah, Xh = np.copy(A), np.copy(X) + # on this basis, get the transform + Ar, Xr, blsize, W = mb03rd( + test1_n, A, X, 'U', 'S', test1_pmax, test1_tol) + # ensure X and A are unchanged + assert_allclose(A, Ah) + assert_allclose(X, Xh) + # compare to test case results + assert_allclose(Ar, test1_Ar, atol=0.0001) + assert_allclose(Xr, test1_Xr, atol=0.0001) + assert_allclose(W, test1_W, atol=0.0001) + + # Test that the non sorting options do not throw errors and that Xr is + # returned as None for jobx='N' + for sort in ['N', 'C', 'B']: + Ar, Xr, blsize, W = mb03rd( + test1_n, A, X, 'N', sort, test1_pmax, test1_tol) + assert Xr is None + def test_mb03vd_mb03vy_ex(self): """Test MB03VD and MB03VY with the example given in the MB03VD SLICOT documentation""" diff --git a/slycot/tests/test_mb03rd.py b/slycot/tests/test_mb03rd.py deleted file mode 100644 index 6c820871..00000000 --- a/slycot/tests/test_mb03rd.py +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env python -# -# test_mb03rd.py - test suite for Shur form reduction -# RvP, 31 Jul 2019 -import unittest -from slycot import transform -import numpy as np -from numpy.testing import assert_raises, assert_almost_equal, assert_equal -from scipy.linalg import schur - -test1_A = np.array([ - [ 1., -1., 1., 2., 3., 1., 2., 3.], - [ 1., 1., 3., 4., 2., 3., 4., 2.], - [ 0., 0., 1., -1., 1., 5., 4., 1.], - [ 0., 0., 0., 1., -1., 3., 1., 2.], - [ 0., 0., 0., 1., 1., 2., 3., -1.], - [ 0., 0., 0., 0., 0., 1., 5., 1.], - [ 0., 0., 0., 0., 0., 0., 0.99999999, -0.99999999 ], - [ 0., 0., 0., 0., 0., 0., 0.99999999, 0.99999999 ] - ]) -test1_n = test1_A.shape[0] - -test1_Ar = np.array([ - [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186, 1.4577, 0.0000, 0.0000 ], - [ 1.0000, 1.0000, 0.0000, 1.4142, -5.1390, 3.1637, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 1.0000, -1.7321, -0.0016, 2.0701, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.5774, 1.0000, 0.7516, 1.1379, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -5.8606, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.1706, 1.0000, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, -0.8850 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ], - ]) - -test1_Xr = np.array([ - [ 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.9045, 0.1957 ], - [ 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3015, 0.9755 ], - [ 0.0000, 0.0000, 0.8165, 0.0000, -0.5768, -0.0156, -0.3015, 0.0148 ], - [ 0.0000, 0.0000, -0.4082, 0.7071, -0.5768, -0.0156, 0.0000, -0.0534 ], - [ 0.0000, 0.0000, -0.4082, -0.7071, -0.5768, -0.0156, 0.0000, 0.0801 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, -0.0276, 0.9805, 0.0000, 0.0267 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0332, -0.0066, 0.0000, 0.0000 ], - [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0011, 0.1948, 0.0000, 0.0000 ] - ]) - -test1_pmax = 1e3 -test1_tol = 0.01 -class test_mb03rd(unittest.TestCase): - def test1(self): - # create schur form with scipy - A, X = schur(test1_A) - Ah, Xh = np.copy(A), np.copy(X) - # on this basis, get the transform - Ar, Xr, blks, eig = transform.mb03rd( - test1_n, A, X, 'U', 'S', test1_pmax, test1_tol) - # ensure X and A are unchanged - assert_equal(A, Ah) - assert_equal(X, Xh) - # compare to test case results - assert_almost_equal(Ar, test1_Ar, decimal=4) - assert_almost_equal(Xr, test1_Xr, decimal=4) - -if __name__ == "__main__": - unittest.main() diff --git a/slycot/transform.py b/slycot/transform.py index 52357d72..0e758f60 100644 --- a/slycot/transform.py +++ b/slycot/transform.py @@ -1049,30 +1049,30 @@ def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None): def tg01ad(l,n,m,p,A,E,B,C,thresh=0.0,job='A'): """ A,E,B,C,lscale,rscale = tg01ad(l,n,m,p,A,E,B,C,[thresh,job]) - + To balance the matrices of the system pencil - + S = ( A B ) - lambda ( E 0 ) := Q - lambda Z, ( C 0 ) ( 0 0 ) - + corresponding to the descriptor triple (A-lambda E,B,C), by balancing. This involves diagonal similarity transformations (Dl*A*Dr - lambda Dl*E*Dr, Dl*B, C*Dr) applied to the system (A-lambda E,B,C) to make the rows and columns of system pencil matrices - + diag(Dl,I) * S * diag(Dr,I) - + as close in norm as possible. Balancing may reduce the 1-norms of the matrices of the system pencil S. - + The balancing can be performed optionally on the following particular system pencils - + S = A-lambda E, - + S = ( A-lambda E B ), or - + S = ( A-lambda E ). ( C ) Required arguments: @@ -1134,15 +1134,15 @@ def tg01ad(l,n,m,p,A,E,B,C,thresh=0.0,job='A'): the scaling factor applied to column j, then SCALE(j) = Dr(j), for j = 1,...,N. """ - + hidden = ' (hidden by the wrapper)' arg_list = ['job', 'l', 'n', 'm', 'p', 'thresh', 'A', 'lda'+hidden, 'E','lde'+hidden,'B','ldb'+hidden,'C','ldc'+hidden, 'lscale', 'rscale', 'dwork'+hidden, 'info'] - + if job != 'A' and job != 'B' and job != 'C' and job != 'N': raise ValueError('Parameter job had an illegal value') A,E,B,C,lscale,rscale,info = _wrapper.tg01ad(job,l,n,m,p,thresh,A,E,B,C) - + if info < 0: error_text = "The following argument had an illegal value: "+arg_list[-info-1] e = ValueError(error_text) @@ -1152,7 +1152,7 @@ def tg01ad(l,n,m,p,A,E,B,C,thresh=0.0,job='A'): e = ArithmeticError('tg01ad failed') e.info = info raise e - + return A,E,B,C,lscale,rscale def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ldwork=None): @@ -1162,23 +1162,23 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld the orthogonal transformation matrices Q and Z such that the transformed system (Q'*A*Z-lambda Q'*E*Z, Q'*B, C*Z) is in a SVD-like coordinate form with - + ( A11 A12 ) ( Er 0 ) Q'*A*Z = ( ) , Q'*E*Z = ( ) , ( A21 A22 ) ( 0 0 ) - + where Er is an upper triangular invertible matrix. Optionally, the A22 matrix can be further reduced to the form - + ( Ar X ) A22 = ( ) , ( 0 0 ) - + with Ar an upper triangular invertible matrix, and X either a full or a zero matrix. The left and/or right orthogonal transformations performed to reduce E and A22 can be optionally accumulated. - + Required arguments: l : input int The number of rows of matrices A, B, and E. l >= 0. @@ -1253,22 +1253,22 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld On exit, the leading L-by-N part of this array contains the transformed matrix Q'*A*Z. If JOBA = 'T', this matrix is in the form - + ( A11 * * ) Q'*A*Z = ( * Ar X ) , ( * 0 0 ) - + where A11 is a RANKE-by-RANKE matrix and Ar is a RNKA22-by-RNKA22 invertible upper triangular matrix. If JOBA = 'R' then A has the above form with X = 0. E : rank-2 array('d') with bounds (l,n) The leading L-by-N part of this array contains the transformed matrix Q'*E*Z. - + ( Er 0 ) Q'*E*Z = ( ) , ( 0 0 ) - + where Er is a RANKE-by-RANKE upper triangular invertible matrix. B : rank-2 array('d') with bounds (l,m) @@ -1310,7 +1310,7 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld hidden = ' (hidden by the wrapper)' arg_list = ['compq', 'compz', 'joba', 'l', 'n', 'm', 'p', 'A', 'lda'+hidden, 'E','lde'+hidden,'B','ldb'+hidden,'C','ldc'+hidden,'Q','ldq'+hidden,'Z','ldz'+hidden,'ranke','rnka22','tol','iwork'+hidden, 'dwork'+hidden, 'ldwork', 'info'] - + if compq != 'N' and compq != 'I' and compq != 'U': raise ValueError('Parameter compq had an illegal value') @@ -1334,7 +1334,7 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld A,E,B,C,Q,Z,ranke,rnka22,info = _wrapper.tg01fd_uu(joba,l,n,m,p,A,E,B,C,Q,Z,tol,ldwork) else: raise ValueError("The combination of compq and compz in not implemented") - + if info < 0: error_text = "The following argument had an illegal value: "+arg_list[-info-1] e = ValueError(error_text) @@ -1344,116 +1344,10 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld e = ArithmeticError('tg01fd failed') e.info = info raise e - + if joba == 'N': rnka22 = None - - return A,E,B,C,ranke,rnka22,Q,Z - -def mb03rd(n,A,X=None,jobx='U',sort='N',pmax=1.0,tol=0.0): - """ A,X,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='U' - A,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='N' - - 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. - Required arguments: - n : input int - The order of the matrices A and X. n >= 0. - A : input rank-2 array('d') with bounds (n,n) - the matrix A to be block-diagonalized, in real Schur form. - Optional arguments: - X : input rank-2 array('d') with bounds (n,n) - a given matrix X, for accumulation of transformations (only if - jobx='U' - jobx : input char*1 - Specifies whether or not the transformations are - accumulated, as follows: - = 'N': The transformations are not accumulated; - = 'U': The transformations are accumulated in X (the - given matrix X is updated). - sort : input char*1 - Specifies whether or not the diagonal blocks of the real - Schur form are reordered, as follows: - = 'N': The diagonal blocks are not reordered; - = '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 METHOD); - = 'B': The diagonal blocks are reordered before each - step of reduction, and the "closest-neighbour" - strategy is used (see METHOD). - pmax : input float - An upper bound for the infinity norm of elementary - submatrices of the individual transformations used for - reduction (see METHOD). PMAX >= 1.0D0. - tol : input float - 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 - - | 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, - - | 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. - Return objects: - Ar : output rank-2 array('d') with bounds (n,n) - Contains the computed block-diagonal matrix, in real Schur - canonical form. The non-diagonal blocks are set to zero. - Xr : output rank-2 array('d') with bounds (n,n) - 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 not referenced, and not returned - blksize : output rank-1 array('i') with bounds (n) - The orders of the resulting diagonal blocks of the matrix Ar. - W : output rank-1 array('c') size (n) - This arrays contain the eigenvalues of the matrix A. -""" - hidden = ' (hidden by the wrapper)' - arg_list = ('jobx', 'sort', 'n', 'pmax', 'A', 'LDA'+hidden, - 'X', 'LDX'+hidden, 'nblks'+hidden, 'blsize'+hidden, - 'WR'+hidden, 'WI'+hidden, 'tol', - 'DWORK'+hidden, 'INFO'+hidden) - if jobx == 'N': - out = _wrapper.mb03rd_n(sort, n, pmax, A, tol) - else: - if X is None: - X = _np.eye(n) - out = _wrapper.mb03rd_u(sort, n, pmax, A, X, tol) - - if out[-1] < 0: - error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] - e = ValueError(error_text) - e.info = out[-1] - raise e - if jobx == 'N': - return out[0], out[2][:out[1]], out[-3] + out[-2]*1j - else: - return out[0], out[1], out[3][:out[2]], out[-3] + out[-2]*1j + return A,E,B,C,ranke,rnka22,Q,Z # to be replaced by python wrappers