Skip to content

Work on build system #90

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

Closed
repagh wants to merge 2 commits into from
Closed

Work on build system #90

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
009a2b8
adding the mb03rd wrapper
repagh Aug 2, 2019
ba2b2f9
Merge remote-tracking branch 'python-control/master'
repagh Jan 25, 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
37 changes: 37 additions & 0 deletions slycot/src/transform.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

63 changes: 63 additions & 0 deletions slycot/tests/test_mb03rd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()
106 changes: 106 additions & 0 deletions slycot/transform.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1352,4 +1352,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