diff --git a/slycot/CMakeLists.txt b/slycot/CMakeLists.txt index 892e16de..9819de59 100644 --- a/slycot/CMakeLists.txt +++ b/slycot/CMakeLists.txt @@ -104,7 +104,7 @@ set(FSOURCES set(F2PYSOURCE src/_wrapper.pyf) set(F2PYSOURCE_DEPS - src/analysis.pyf src/math.pyf src/mathematical.pyf + src/analysis.pyf src/math.pyf src/transform.pyf src/synthesis.pyf) configure_file(version.py.in version.py @ONLY) diff --git a/slycot/math.py b/slycot/math.py index 7044e2bb..8ad0b020 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -18,6 +18,7 @@ # MA 02110-1301, USA. from . import _wrapper +import warnings def mc01td(dico,dp,p): """ dp,stable,nz = mc01td(dico,dp,p) @@ -61,16 +62,16 @@ def mc01td(dico,dp,p): e.info = out[-1] raise e if out[-1] == 1: - warings.warn('entry P(x) is the zero polynomial.') + warnings.warn('entry P(x) is the zero polynomial.') if out[-1] == 2: - warings.warn('P(x) may have zeros very close to stability boundary.') + warnings.warn('P(x) may have zeros very close to stability boundary.') if out[-2] > 0: - warnings.warn('The degree of P(x) has been reduced to %i' %(dp-k)) + warnings.warn('The degree of P(x) has been reduced to %i' %(dp-out[-2])) return out[:-2] def mb05md(a, delta, balanc='N'): - """Ar, Vr, Yr, VALRr, VALDr = mb05md(a, delta, balanc='N') + """Ar, Vr, Yr, VAL = mb05md(a, delta, balanc='N') Matrix exponential for a real non-defective matrix @@ -88,7 +89,7 @@ def mb05md(a, delta, balanc='N'): Square matrix delta : input 'd' The scalar value delta of the problem. - + Optional arguments: balanc : input char*1 Indicates how the input matrix should be diagonally scaled @@ -114,7 +115,7 @@ def mb05md(a, delta, balanc='N'): (k+1)-th columns of the eigenvector matrix, respectively, then the eigenvector corresponding to the complex eigenvalue with positive (negative) imaginary value is - given by + given by p + q*j (p - q*j), where j^2 = -1. Yr : output rank-2 array('d') with bounds (n,n) contains an intermediate result for computing the matrix @@ -126,29 +127,39 @@ def mb05md(a, delta, balanc='N'): the (right) eigenvector matrix of A, where Lambda is the diagonal matrix of eigenvalues. - VALr : output rank-1 array('c') with bounds (n) + VAL : output rank-1 array('c') with bounds (n) Contains the eigenvalues of the matrix A. The eigenvalues are unordered except that complex conjugate pairs of values appear consecutively with the eigenvalue having positive imaginary part first. """ hidden = ' (hidden by the wrapper)' - arg_list = [ 'balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden, - 'y','ldy'+hidden,'valr','vali', - 'iwork'+hidden,'dwork'+hidden,'ldwork'+hidden,'info'+hidden] - out = _wrapper.mb05md(balanc=balanc,n=min(a.shape),delta=delta,a=a) - if out[-1] == 0: - return out[:-1] - elif out[-1] < 0: - error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] - elif out[-1] > 0 and out[-1] <= n: - error_text = "Incomplete eigenvalue calculation, missing %i eigenvalues" % out[-1] - elif out[-1] == n+1: + arg_list = ['balanc', 'n', 'delta', 'a', 'lda'+hidden, 'v', 'ldv'+hidden, + 'y', 'ldy'+hidden, 'valr', 'vali', + 'iwork'+hidden, 'dwork'+hidden, 'ldwork'+hidden, + 'info'+hidden] + n = min(a.shape) + (Ar, Vr, Yr, VALr, VALi, INFO) = _wrapper.mb05md(balanc=balanc, + n=n, + delta=delta, + a=a) + if INFO == 0: + if not all(VALi == 0): + VAL = VALr + 1J*VALi + else: + VAL = VALr + return (Ar, Vr, Yr, VAL) + elif INFO < 0: + error_text = "The following argument had an illegal value: " \ + + arg_list[-INFO-1] + elif INFO > 0 and INFO <= n: + error_text = "Incomplete eigenvalue calculation, missing %i eigenvalues" % INFO + elif INFO == n+1: error_text = "Eigenvector matrix singular" - elif out[-1] == n+2: + elif INFO == n+2: error_text = "A matrix defective" e = ValueError(error_text) - e.info = out[-1] + e.info = INFO raise e """ @@ -182,21 +193,21 @@ def mb05nd(a, delta, tol=1e-7): H : Int[F(s) ds] from s = 0 to s = delta, """ hidden = ' (hidden by the wrapper)' - arg_list = [ 'n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden, - 'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden, - 'dwork'+hidden, 'ldwork'+hidden] - out = _wrapper.mb05nd(n=min(a.shape), delta=delta, a=a, tol=tol) + arg_list = ['n', 'delta', 'a', 'lda'+hidden, 'ex', 'ldex'+hidden, + 'exint', 'ldexin'+hidden, 'tol', 'iwork'+hidden, + 'dwork'+hidden, 'ldwork'+hidden] + n = min(a.shape) + out = _wrapper.mb05nd(n=n, delta=delta, a=a, tol=tol) if out[-1] == 0: return out[:-1] elif out[-1] < 0: error_text = "The following argument had an illegal value: " \ - +arg_list[-out[-1]-1] + + arg_list[-out[-1]-1] elif out[-1] == n+1: error_text = "Delta too large" e = ValueError(error_text) e.info = out[-1] raise e - # to be replaced by python wrappers diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 2a290939..101e20e2 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -3,7 +3,7 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f character :: dico - integer intent(in,out),check(dp>0) :: dp + integer intent(in,out),check(dp>=0) :: dp double precision intent(in),check(shape(p,0)==dp+1),dimension(dp+1),depend(dp) :: p logical intent(out) :: stable integer intent(out) :: nz diff --git a/slycot/src/mathematical.pyf b/slycot/src/mathematical.pyf deleted file mode 100644 index ab318846..00000000 --- a/slycot/src/mathematical.pyf +++ /dev/null @@ -1,23 +0,0 @@ -! -*- f90 -*- -! Note: the context of this file is case sensitive. - -subroutine mb05md(balanc,n,delta,a,lda,v,ldv,y,ldy,valr,vali,iwork,dwork,ldwork,info) ! in MB05MD.f - character check(balanc=='N' || balanc=='S'):: balanc - integer check(n>=0) :: n - double precision check(delta>=0.0):: delta - double precision intent(in,out,copy),dimension(lda,*) :: a - integer, intent(hide),depend(a) :: lda=shape(a,0) - double precision intent(out,copy),dimension(ldv,*) :: v - integer, intent(hide),check(ldv),depend(v) :: ldv=shape(v,0) - double precision dimension(ldy,*) :: y - integer, optional,check(shape(y,0)==ldy),depend(y) :: ldy=shape(y,0) - double precision dimension(*) :: valr - double precision dimension(*) :: vali - integer dimension(*) :: iwork - double precision dimension(*) :: dwork - integer :: ldwork - integer :: info -end subroutine mb05md - -! This file was auto-generated with f2py (version:2). -! See http://cens.ioc.ee/projects/f2py2e/ diff --git a/slycot/synthesis.py b/slycot/synthesis.py index b197c3b7..f80c22c5 100644 --- a/slycot/synthesis.py +++ b/slycot/synthesis.py @@ -689,20 +689,20 @@ def sb02od(n,m,A,B,Q,R,dico,p=None,L=None,fact='N',uplo='U',sort='S',tol=0.0,ldw out = _wrapper.sb02od_n(dico,n,m,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork) if fact == 'C': if p is None: - p = shape(Q)[0] + p = _np.shape(Q)[0] out = _wrapper.sb02od_c(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork) if fact == 'D': if p is None: - p = shape(R)[0] + p = _np.shape(R)[0] out = _wrapper.sb02od_d(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork) if fact == 'B': if p is None: - p = shape(Q)[0] + p = _np.shape(Q)[0] out = _wrapper.sb02od_b(dico,n,m,p,A,B,Q,R,L,uplo=uplo,jobl=jobl,sort=sort,tol=tol,ldwork=ldwork) if out[-1] < 0: error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1] e = ValueError(error_text) - e.info = info + e.info = out[-1] raise e if out[-1] == 1: e = ValueError('the computed extended matrix pencil is singular, possibly due to rounding errors') diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 7bdd4dc0..8c7f26d1 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -3,6 +3,8 @@ set(PYSOURCE __init__.py test.py test_ag08bd.py + test_mb.py + test_mc.py test_sb10jd.py test_sg02ad.py test_sg03ad.py @@ -11,4 +13,8 @@ set(PYSOURCE test_tg01ad.py test_tg01fd.py ) -install(FILES ${PYSOURCE} DESTINATION slycot/tests) +install(FILES ${PYSOURCE} + PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE + GROUP_READ GROUP_EXECUTE + WORLD_READ WORLD_EXECUTE + DESTINATION slycot/tests) diff --git a/slycot/tests/test.py b/slycot/tests/test.py index ec0b0720..19f7dfc8 100644 --- a/slycot/tests/test.py +++ b/slycot/tests/test.py @@ -1,6 +1,5 @@ import unittest from slycot import synthesis -from slycot import math from slycot import transform class Test(unittest.TestCase): @@ -11,11 +10,6 @@ def setUp(self): def test_1(self): synthesis.sb02mt(1,1,1,1) - def test_2(self): - from numpy import array - a = array([[-2, 0.5], [-1.6, -5]]) - Ar, Vr, Yr, VALRr, VALDr = math.mb05md(a, 0.1) - def test_sb02ad(self): "Test sb10ad, Hinf synthesis" import numpy as np @@ -64,9 +58,5 @@ def test_td04ad_static(self): nr,a,b,c,d = transform.td04ad(rc,nin,nout,index,den,num) -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestConvert) - - if __name__ == "__main__": unittest.main() diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py new file mode 100755 index 00000000..2f8e3750 --- /dev/null +++ b/slycot/tests/test_mb.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python +# +# test_mb.py - test suite for linear algebra commands +# bnavigator , Aug 2019 + +import unittest +import numpy as np + +from slycot import mb05md, mb05nd + +from numpy.testing import assert_allclose + + +class test_mb(unittest.TestCase): + + def test_mb05md(self): + """ test_mb05md: verify Matrix exponential with slicot doc example + data from http://slicot.org/objects/software/shared/doc/MB05MD.html + """ + A = np.array([[ 0.5, 0., 2.3, -2.6], + [ 0., 0.5, -1.4, -0.7], + [ 2.3, -1.4, 0.5, 0.0], + [-2.6, -0.7, 0.0, 0.5]]) + delta = 1.0 + Ar_ref = np.array([[ 26.8551, -3.2824, 18.7409, -19.4430], + [ -3.2824, 4.3474, -5.1848, 0.2700], + [ 18.7409, -5.1848, 15.6012, -11.7228], + [ -19.4430, 0.2700, -11.7228, 15.6012]]) + Vr_ref = np.array([[-0.7, 0.7, 0.1, -0.1], + [ 0.1, -0.1, 0.7, -0.7], + [ 0.5, 0.5, 0.5, 0.5], + [-0.5, -0.5, 0.5, 0.5]]) + Yr_ref = np.array([[ -0.0349, 0.0050, 0.0249, -0.0249], + [ 38.2187, -5.4598, 27.2991, -27.2991], + [ 0.0368, 0.2575, 0.1839, 0.1839], + [ -0.7389, -5.1723, 3.6945, 3.6945]]) + VAL_ref = np.array([-3., 4., -1., 2.]) + (Ar, Vr, Yr, VAL) = mb05md(A, delta) + + assert_allclose(Ar, Ar_ref, atol=0.0001) + + # Order of eigenvalues is not guaranteed, so we check them one by one. + for i, e in enumerate(VAL): + erow = np.ones(VAL.shape)*e + i_ref = np.isclose(erow, VAL_ref) + self.assertTrue(any(i_ref), + msg="eigenvalue {} not expected".format(e)) + # Eigenvectors can have different scaling. + vr_ref = Vr_ref[:, i_ref]*Vr[0, i]/Vr_ref[0, i_ref][0] + assert_allclose(Vr[:, (i,)], vr_ref, atol=0.0001) + + assert_allclose(np.dot(Vr, Yr), np.dot(Vr_ref, Yr_ref), atol=0.0001) + + def test_mb05nd(self): + """ test_mb05nd: verify Matrix exponential and integral + data from http://slicot.org/objects/software/shared/doc/MB05ND.html + """ + A = np.array([[5.0, 4.0, 3.0, 2.0, 1.0], + [1.0, 6.0, 0.0, 4.0, 3.0], + [2.0, 0.0, 7.0, 6.0, 5.0], + [1.0, 3.0, 1.0, 8.0, 7.0], + [2.0, 5.0, 7.0, 1.0, 9.0]]) + delta = 0.1 + F_ref = np.array([[1.8391, 0.9476, 0.7920, 0.8216, 0.7811], + [0.3359, 2.2262, 0.4013, 1.0078, 1.0957], + [0.6335, 0.6776, 2.6933, 1.6155, 1.8502], + [0.4804, 1.1561, 0.9110, 2.7461, 2.0854], + [0.7105, 1.4244, 1.8835, 1.0966, 3.4134]]) + H_ref = np.array([[0.1347, 0.0352, 0.0284, 0.0272, 0.0231], + [0.0114, 0.1477, 0.0104, 0.0369, 0.0368], + [0.0218, 0.0178, 0.1624, 0.0580, 0.0619], + [0.0152, 0.0385, 0.0267, 0.1660, 0.0732], + [0.0240, 0.0503, 0.0679, 0.0317, 0.1863]]) + + (F, H) = mb05nd(A, delta) + + assert_allclose(F, F_ref, atol=0.0001) + assert_allclose(H, H_ref, atol=0.0001) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/tests/test_mc.py b/slycot/tests/test_mc.py new file mode 100644 index 00000000..5e703f87 --- /dev/null +++ b/slycot/tests/test_mc.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python +# +# test_mc.py - test suite for polynomial and rational function manipulation +# bnavigator , Aug 2019 + +import unittest +import warnings + +from slycot import mc01td + + +class test_mc(unittest.TestCase): + + def test_mc01td(self): + """ test_mc01td: doc example + data from http://slicot.org/objects/software/shared/doc/MC01TD.html + """ + (dp, stable, nz) = mc01td('C', 4, [2, 0, 1, -1, 1]) + self.assertEqual(dp, 4) + self.assertEqual(stable, 0) + self.assertEqual(nz, 2) + + def test_mc01td_D(self): + """ test_mc01td_D: test discrete option """ + (dp, stable, nz) = mc01td('D', 3, [1, 2, 3, 4]) + self.assertEqual(dp, 3) + self.assertEqual(stable, 1) + self.assertEqual(nz, 0) + (dp, stable, nz) = mc01td('D', 3, [4, 3, 2, 1]) + self.assertEqual(dp, 3) + self.assertEqual(stable, 0) + self.assertEqual(nz, 3) + + def test_mc01td_warnings(self): + """ test_mc01td_warnings: Test warnings """ + T = [([0, 0], "entry P(x) is the zero polynomial."), + ([0, 1], "P(x) may have zeros very close to stability boundary."), + ([1, 0], "The degree of P(x) has been reduced to 0")] + for P, m in T: + with warnings.catch_warnings(record=True) as w: + (dp, stable, nz) = mc01td('C', len(P)-1, P) + self.assertEqual(str(w[0].message), m) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/transform.py b/slycot/transform.py index 77b31715..d492a075 100644 --- a/slycot/transform.py +++ b/slycot/transform.py @@ -649,7 +649,6 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None): kdcoef = max(index)+1 if rowcol == 'R': - porm = p if ucoeff.ndim != 3: e = ValueError("The numerator is not a 3D array!") e.info = -7 @@ -664,7 +663,6 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None): raise e out = _wrapper.td04ad_r(m,p,index,dcoeff,ucoeff,n,tol,ldwork) elif rowcol == 'C': - porm = m if ucoeff.ndim != 3: e = ValueError("The numerator is not a 3D array!") e.info = -7