Skip to content

addition of wrapper for mb03rd #73

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

addition of wrapper for mb03rd #73

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
1d5a7e1
addition of wrapper for mb03rd
repagh Aug 2, 2019
89fa119
need to install scipy for the mb03rd test
repagh Aug 3, 2019
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 .travis.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,12 +150,12 @@ script:
# Local unit tests
# TODO: replace with nose?
- cd ..
# Additional packages required for testing
- conda install scipy matplotlib
Copy link
Collaborator

@bnavigator bnavigator Feb 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

adressed in #93

- python Slycot/runtests.py --coverage --no-build
#
# As a deeper set of tests, get test against python-control as well
#
# Additional packages required for python-control
- conda install scipy matplotlib
# Get python-control from source and install
- git clone --depth 1 https://github.com/python-control/python-control.git control
- cd control
Expand Down
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
Copy link
Collaborator

@bnavigator bnavigator Feb 1, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It may be just my personal opinion, but I don't think it is a good idea to introduce a specific wrapper for different job options with different parameter and return signatures when there is no use case and no unit test for it. Why not just keep one wrapper for the function, test that with the default job='U' and deal with 'N' as soon as it becomes relevant?

Copy link
Collaborator

@bnavigator bnavigator Feb 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Move to math.pyf?

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(
Copy link
Collaborator

@bnavigator bnavigator Feb 8, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Assert correct eig as well (complex?)

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 @@ -1353,4 +1353,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):
Copy link
Collaborator

@bnavigator bnavigator Feb 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Belongs into math.py

""" A,X,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='U'
Copy link
Collaborator

@bnavigator bnavigator Feb 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually returns _,_,_,Wr,Wi not one EIG. Would suggest converting to complex array as in #88

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)
Copy link
Collaborator

@bnavigator bnavigator Feb 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does not fit the variable names in the call signature blcks,EIG

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