diff --git a/slycot/__init__.py b/slycot/__init__.py index 9d8704a8..b9bde0a1 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -62,14 +62,16 @@ sg02ad, sg03ad, sg03bd) - # Transformation routines (10/77 wrapped) + # Transformation routines (12/77 wrapped) from .transform import (tb01id, tb01pd, tb03ad, tb04ad, tb05ad, - tc01od, tc04ad, + tc01od, + tc04ad, td04ad, - tf01md, tf01rd) + tf01md, tf01rd, + tg01ad, tg01fd) # Utility routines (0/7 wrapped) diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 80d751da..05ecb6c5 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -11,11 +11,18 @@ set(PYSOURCE test_mc.py test_sb.py test_analysis.py - test_transform.py test_sg02ad.py test_sg03ad.py + test_tb01id.py + test_tb01pd.py + test_tb03ad.py + test_tb04ad.py test_tb05ad.py + test_tc01od.py + test_tc04ad.py test_td04ad.py + test_tf01md.py + test_tf01rd.py test_tg01ad.py test_tg01fd.py ) diff --git a/slycot/tests/test_tb01id.py b/slycot/tests/test_tb01id.py new file mode 100644 index 00000000..efb6bcb7 --- /dev/null +++ b/slycot/tests/test_tb01id.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb01id: + + @pytest.mark.skip(reason="not implemented") + def test_tb01id(self): + pass diff --git a/slycot/tests/test_tb01pd.py b/slycot/tests/test_tb01pd.py new file mode 100644 index 00000000..f4018b8a --- /dev/null +++ b/slycot/tests/test_tb01pd.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb01pd: + + @pytest.mark.skip(reason="not implemented") + def test_tb01pd(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb03ad.py b/slycot/tests/test_tb03ad.py new file mode 100644 index 00000000..a9cb6024 --- /dev/null +++ b/slycot/tests/test_tb03ad.py @@ -0,0 +1,21 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb03ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tb03ad, SlycotArithmeticError, 2, {}),)) + def test_tb03ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + @pytest.mark.skip(reason="not implemented") + def test_tb03ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb04ad.py b/slycot/tests/test_tb04ad.py new file mode 100644 index 00000000..7d241333 --- /dev/null +++ b/slycot/tests/test_tb04ad.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb04ad: + + @pytest.mark.skip(reason="not implemented") + def test_tb04ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb05ad.py b/slycot/tests/test_tb05ad.py index 900a49a2..df5996d1 100644 --- a/slycot/tests/test_tb05ad.py +++ b/slycot/tests/test_tb05ad.py @@ -1,6 +1,3 @@ -# =================================================== -# tb05ad tests - import sys import numpy as np @@ -9,8 +6,11 @@ from scipy.linalg import eig, matrix_balance from slycot import transform +from slycot import transform as tf from slycot.exceptions import SlycotArithmeticError, SlycotParameterError +from .test_exceptions import assert_docstring_parse + # set the random seed so we can get consistent results. np.random.seed(40) CASES = {} @@ -36,216 +36,226 @@ 'B': np.random.randn(n, m), 'C': np.random.randn(p, n)} - -def test_tb05ad_ng(): - """ - Test that tb05ad with job 'NG' computes the correct - frequency response. - """ - for key in CASES: - sys = CASES[key] - check_tb05ad_AG_NG(sys, 10*1j, 'NG') - - -def test_tb05ad_ag(): - """ - Test that tb05ad with job 'AG' computes the correct - frequency response. - """ - for key in CASES: - sys = CASES[key] - check_tb05ad_AG_NG(sys, 10*1j, 'AG') - - -def test_tb05ad_nh(): - """Test that tb05ad with job = 'NH' computes the correct - frequency response after conversion to Hessenberg form. - - First call tb05ad with job='NH' to transform to upper Hessenberg - form which outputs the transformed system. - Subsequently, call tb05ad with job='NH' using this transformed system. - """ - jomega = 10*1j - for key in CASES: - sys = CASES[key] - sys_transformed = check_tb05ad_AG_NG(sys, jomega, 'NG') - check_tb05ad_NH(sys_transformed, sys, jomega) - - -def test_tb05ad_errors(): - """ - Test tb05ad error handling. We give wrong inputs and - and check that this raises an error. - """ - check_tb05ad_errors(CASES['pass1']) - - -def check_tb05ad_AG_NG(sys, jomega, job): - """ - Check that tb05ad computes the correct frequency response when - running jobs 'AG' and/or 'NG'. - - Inputs - ------ - - sys: A a dict of system matrices with keys 'A', 'B', and 'C'. - jomega: A complex scalar, which is the frequency we are - evaluating the system at. - job: A string, either 'AG' or 'NH' - - Returns - ------- - sys_transformed: A dict of the system matrices which have been - transformed according the job. - """ - n, m = sys['B'].shape - p = sys['C'].shape[0] - result = transform.tb05ad(n, m, p, jomega, - sys['A'], sys['B'], sys['C'], job=job) - g_i = result[3] - hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) - g_i_solve = sys['C'].dot(hinvb) - assert_almost_equal(g_i_solve, g_i) - sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]} - return sys_transformed - - -def check_tb05ad_NH(sys_transformed, sys, jomega): - """ - Check tb05ad, computes the correct frequency response when - job='NH' and we supply system matrices 'A', 'B', and 'C' - which have been transformed by a previous call to tb05ad. - We check we get the same result as computing C(sI - A)^-1B - with the original system. - - Inputs - ------ - - sys_transformed: A a dict of the transformed (A in upper - hessenberg form) system matrices with keys - 'A', 'B', and 'C'. - - sys: A dict of the original un-transformed system matrices. - - jomega: A complex scalar, which is the frequency to evaluate at. - - """ - - n, m = sys_transformed['B'].shape - p = sys_transformed['C'].shape[0] - result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'], - sys_transformed['B'], sys_transformed['C'], - job='NH') - g_i = result[0] - hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) - g_i_solve = sys['C'].dot(hinvb) - assert_almost_equal(g_i_solve, g_i) - - -def check_tb05ad_errors(sys): - """ - Check the error handling of tb05ad. We give wrong inputs and - and check that this raises an error. - """ - n, m = sys['B'].shape - p = sys['C'].shape[0] - jomega = 10*1j - # test error handling - # wrong size A - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n+1, m, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -7 - # wrong size B - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m+1, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -9 - # wrong size C - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m, p+1, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -11 - # unrecognized job - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m, p, jomega, sys['A'], sys['B'], sys['C'], job='a') - assert cm.value.info == -1 - - -def test_tb05ad_resonance(): - """ Test tb05ad resonance failure. - - Actually test one of the exception messages. These - are parsed from the docstring, tests both the info index and the - message - """ - A = np.array([[0, -1], - [1, 0]]) - B = np.array([[1], - [0]]) - C = np.array([[0, 1]]) - jomega = 1j - with pytest.raises( - SlycotArithmeticError, - match=r"Either `freq`.* is too near to an eigenvalue of A,\n" - r"or `rcond` is less than the machine precision EPS.") as cm: - transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH') - assert cm.value.info == 2 - - -def test_tb05ad_balance(): - """Test balancing in tb05ad. - - Tests for the cause of the problem reported in issue #11 - balancing permutations were not correctly applied to the - C and D matrix. - """ - - # find a good test case. Some sparsity, - # some zero eigenvalues, some non-zero eigenvalues, - # and proof that the 1st step, with dgebal, does some - # permutation and some scaling - crit = False - n = 8 - while not crit: - A = np.random.randn(n, n) - A[np.random.uniform(size=(n, n)) > 0.35] = 0.0 - - Aeig = eig(A)[0] - neig0 = np.sum(np.abs(Aeig) == 0) - As, T = matrix_balance(A) - nperm = np.sum(np.diag(T == 0)) - nscale = n - np.sum(T == 1.0) - crit = nperm < n and nperm >= n//2 and \ - neig0 > 1 and neig0 <= 3 and nscale > 0 - - # print("number of permutations", nperm, "eigenvalues=0", neig0) - B = np.random.randn(8, 4) - C = np.random.randn(3, 8) - - # do a run - jomega = 1.0 - At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad( - 8, 4, 3, jomega, A, B, C, job='AG') - - # remove information on Q, in lower sub-triangle part of A - At = np.triu(At, k=-1) - - # now after the balancing in DGEBAL, and conversion to - # upper Hessenberg form: - # At = Q^T * (P^-1 * A * P ) * Q - # with Q orthogonal - # Ct = C * P * Q - # Bt = Q^T * P^-1 * B - # so test with Ct * At * Bt == C * A * B - # and verify that eigenvalues of both A matrices are close - assert_almost_equal(np.dot(np.dot(Ct, At), Bt), - np.dot(np.dot(C, A), B)) - # uses a sort, there is no guarantee on the order of eigenvalues - eigAt = eig(At)[0] - idxAt = np.argsort(eigAt) - eigA = eig(A)[0] - idxA = np.argsort(eigA) - assert_almost_equal(eigA[idxA], eigAt[idxAt]) +class Test_tb05ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tb05ad, SlycotArithmeticError, 2, {'n30': 90, + 'jomega': 2.0, + 'rcond': 1e-12}),)) + def test_tb05ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + + def test_tb05ad_ng(self): + """ + Test that tb05ad with job 'NG' computes the correct + frequency response. + """ + for key in CASES: + sys = CASES[key] + self.check_tb05ad_AG_NG(sys, 10*1j, 'NG') + + + def test_tb05ad_ag(self): + """ + Test that tb05ad with job 'AG' computes the correct + frequency response. + """ + for key in CASES: + sys = CASES[key] + self.check_tb05ad_AG_NG(sys, 10*1j, 'AG') + + + def test_tb05ad_nh(self): + """Test that tb05ad with job = 'NH' computes the correct + frequency response after conversion to Hessenberg form. + + First call tb05ad with job='NH' to transform to upper Hessenberg + form which outputs the transformed system. + Subsequently, call tb05ad with job='NH' using this transformed system. + """ + jomega = 10*1j + for key in CASES: + sys = CASES[key] + sys_transformed = self.check_tb05ad_AG_NG(sys, jomega, 'NG') + self.check_tb05ad_NH(sys_transformed, sys, jomega) + + + def test_tb05ad_errors(self): + """ + Test tb05ad error handling. We give wrong inputs and + and check that this raises an error. + """ + self.check_tb05ad_errors(CASES['pass1']) + + + def check_tb05ad_AG_NG(self, sys, jomega, job): + """ + Check that tb05ad computes the correct frequency response when + running jobs 'AG' and/or 'NG'. + + Inputs + ------ + + sys: A a dict of system matrices with keys 'A', 'B', and 'C'. + jomega: A complex scalar, which is the frequency we are + evaluating the system at. + job: A string, either 'AG' or 'NH' + + Returns + ------- + sys_transformed: A dict of the system matrices which have been + transformed according the job. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, + sys['A'], sys['B'], sys['C'], job=job) + g_i = result[3] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]} + return sys_transformed + + + def check_tb05ad_NH(self, sys_transformed, sys, jomega): + """ + Check tb05ad, computes the correct frequency response when + job='NH' and we supply system matrices 'A', 'B', and 'C' + which have been transformed by a previous call to tb05ad. + We check we get the same result as computing C(sI - A)^-1B + with the original system. + + Inputs + ------ + + sys_transformed: A a dict of the transformed (A in upper + hessenberg form) system matrices with keys + 'A', 'B', and 'C'. + + sys: A dict of the original un-transformed system matrices. + + jomega: A complex scalar, which is the frequency to evaluate at. + + """ + + n, m = sys_transformed['B'].shape + p = sys_transformed['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'], + sys_transformed['B'], sys_transformed['C'], + job='NH') + g_i = result[0] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + + + def check_tb05ad_errors(self, sys): + """ + Check the error handling of tb05ad. We give wrong inputs and + and check that this raises an error. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + jomega = 10*1j + # test error handling + # wrong size A + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n+1, m, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -7 + # wrong size B + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m+1, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -9 + # wrong size C + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m, p+1, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -11 + # unrecognized job + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m, p, jomega, sys['A'], sys['B'], sys['C'], job='a') + assert cm.value.info == -1 + + + def test_tb05ad_resonance(self): + """ Test tb05ad resonance failure. + + Actually test one of the exception messages. These + are parsed from the docstring, tests both the info index and the + message + """ + A = np.array([[0, -1], + [1, 0]]) + B = np.array([[1], + [0]]) + C = np.array([[0, 1]]) + jomega = 1j + with pytest.raises( + SlycotArithmeticError, + match=r"Either `freq`.* is too near to an eigenvalue of A,\n" + r"or `rcond` is less than the machine precision EPS.") as cm: + transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH') + assert cm.value.info == 2 + + + def test_tb05ad_balance(self): + """Test balancing in tb05ad. + + Tests for the cause of the problem reported in issue #11 + balancing permutations were not correctly applied to the + C and D matrix. + """ + + # find a good test case. Some sparsity, + # some zero eigenvalues, some non-zero eigenvalues, + # and proof that the 1st step, with dgebal, does some + # permutation and some scaling + crit = False + n = 8 + while not crit: + A = np.random.randn(n, n) + A[np.random.uniform(size=(n, n)) > 0.35] = 0.0 + + Aeig = eig(A)[0] + neig0 = np.sum(np.abs(Aeig) == 0) + As, T = matrix_balance(A) + nperm = np.sum(np.diag(T == 0)) + nscale = n - np.sum(T == 1.0) + crit = nperm < n and nperm >= n//2 and \ + neig0 > 1 and neig0 <= 3 and nscale > 0 + + # print("number of permutations", nperm, "eigenvalues=0", neig0) + B = np.random.randn(8, 4) + C = np.random.randn(3, 8) + + # do a run + jomega = 1.0 + At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad( + 8, 4, 3, jomega, A, B, C, job='AG') + + # remove information on Q, in lower sub-triangle part of A + At = np.triu(At, k=-1) + + # now after the balancing in DGEBAL, and conversion to + # upper Hessenberg form: + # At = Q^T * (P^-1 * A * P ) * Q + # with Q orthogonal + # Ct = C * P * Q + # Bt = Q^T * P^-1 * B + # so test with Ct * At * Bt == C * A * B + # and verify that eigenvalues of both A matrices are close + assert_almost_equal(np.dot(np.dot(Ct, At), Bt), + np.dot(np.dot(C, A), B)) + # uses a sort, there is no guarantee on the order of eigenvalues + eigAt = eig(At)[0] + idxAt = np.argsort(eigAt) + eigA = eig(A)[0] + idxA = np.argsort(eigA) + assert_almost_equal(eigA[idxA], eigAt[idxAt]) diff --git a/slycot/tests/test_tc01od.py b/slycot/tests/test_tc01od.py new file mode 100644 index 00000000..9b246ed7 --- /dev/null +++ b/slycot/tests/test_tc01od.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tc01od: + + @pytest.mark.skip(reason="not implemented") + def test_tc01od(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tc04ad.py b/slycot/tests/test_tc04ad.py new file mode 100644 index 00000000..c8aa4fb0 --- /dev/null +++ b/slycot/tests/test_tc04ad.py @@ -0,0 +1,21 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tc04ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tc04ad, SlycotArithmeticError, 1, {'leri': 'L'}),)) + def test_tc04ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + @pytest.mark.skip(reason="not implemented") + def test_tc04ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_td04ad.py b/slycot/tests/test_td04ad.py index 57c557d0..c1947a6a 100644 --- a/slycot/tests/test_td04ad.py +++ b/slycot/tests/test_td04ad.py @@ -1,220 +1,230 @@ -# -# test_td04ad.py - test suite for tf -> ss conversion -# RvP, 04 Jun 2018 - import numpy as np +import pytest from slycot import transform - - -def test_td04ad_c(): - """td04ad: Convert with 'C' option""" - - # for octave: - """ - num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ]; - [3.0, -1.0, 1.0 ], [ 0.0, 1.0 ]; - [0.0, 0.0, 1.0], [ 0.0, 2.0 ] }; - den = { [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; - [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; - [1.0, 0.4, 3.0], [ 1.0, 1.0 ]}; - """ - - m = 2 - p = 3 - d = 3 - num = np.array([ - [ [0.0, 0.0, 1.0], [1.0, 0.0, 0.0] ], - [ [3.0, -1.0, 1.0], [0.0, 1.0, 0.0] ], - [ [0.0, 0.0, 1.0], [0.0, 2.0, 0.0] ] ]) - - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - denc = np.array( - [ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ]) - indc = np.array( - [ 2, 1 ], dtype=int) - - nref = 3 - Aref = np.array([ [-1, 0, 0], - [ 0, -0.4, -0.3], - [ 0, 10, 0] ]) - Bref = np.array([ [0, -1], - [1, 0], - [0, 0] ]) - Cref = np.array([ [1, 0, 0.1], - [-1, -2.2, -0.8], - [-2, 0, 0.1] ]) - Dref = np.array([ [0, 1], - [3, 0], - [0, 0] ]) - - nr, A, B, C, D = transform.td04ad('C', m, p, indc, denc, numc) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - np.testing.assert_equal(nref, nr) - # the returned state space representation is not guaranteed to - # be of one form for all architectures, so we transform back - # to tf and check for equality then - _, _, _, _, _, dcoeff, ucoeff = transform.tb04ad( - nr, m, p, A, B, C, D) - _, _, _, _, _, dcoeffref, ucoeffref = transform.tb04ad( - nref, m, p, Aref, Bref, Cref, Dref) - np.testing.assert_array_almost_equal(dcoeff,dcoeffref) - np.testing.assert_array_almost_equal(ucoeff,ucoeffref) - -def test_td04ad_r(): - """td04ad: Convert with 'R' option - - example program from - http://slicot.org/objects/software/shared/doc/TD04AD.html - """ - - m = 2 - p = 2 - rowcol = 'R' - index = [3, 3] - dcoeff = np.array([ [1.0, 6.0, 11.0, 6.0], [1.0, 6.0, 11.0, 6.0] ]) - - ucoeff = np.array([ [[1.0, 6.0, 12.0, 7.0], [0.0, 1.0, 4.0, 3.0]], - [[0.0, 0.0, 1.0, 1.0], [1.0, 8.0, 20.0, 15.0]] ]) - - nref = 3 - - Aref = np.array([ [ 0.5000, -0.8028, 0.9387], - [ 4.4047, -2.3380, 2.5076], - [-5.5541, 1.6872, -4.1620] ]) - Bref = np.array([ [-0.2000, -1.2500], - [ 0.0000, -0.6097], - [ 0.0000, 2.2217] ]) - Cref = np.array([ [0.0000, -0.8679, 0.2119], - [0.0000, 0.0000, 0.9002] ]) - Dref = np.array([ [1.0000, 0.0000], - [0.0000, 1.0000] ]) - - nr, A, B, C, D = transform.td04ad(rowcol, m, p, index, dcoeff, ucoeff) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - np.testing.assert_equal(nref, nr) - # order of states is not guaranteed, so we reorder the reference - rindex = np.flip(np.argsort(np.diag(A))) - Arref = Aref[rindex, :][:, rindex] - Brref = Bref[rindex, :] - Crref = Cref[:, rindex] - Drref = Dref - np.testing.assert_array_almost_equal(A, Arref,decimal=4) - np.testing.assert_array_almost_equal(B, Brref,decimal=4) - np.testing.assert_array_almost_equal(C, Crref,decimal=4) - np.testing.assert_array_almost_equal(D, Drref,decimal=4) - - -def test_staticgain(): - """td04ad: Convert a transferfunction to SS with only static gain""" - - # 2 inputs, 3 outputs? columns share a denominator - num = np.array([ [ [1.0], [2.0] ], - [ [0.2], [4.3] ], - [ [1.2], [3.2] ] ]) - p, m, d = num.shape - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - - # denc, columns share a common denominator - denc = np.array([ [ 1.0], [0.5] ]) - Dc = (num / denc).reshape((3,2)) - idxc = np.zeros((2,), dtype=int) - - # denr, rows share a common denominator - denr = np.array([ [1.0], [0.5], [3.0] ]) - idxr = np.zeros((3,), dtype=int) - Dr = (num / denr[:, np.newaxis]).reshape((3,2)) - - # fails with: - # On entry to TB01XD parameter number 5 had an illegal value - - n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - assert A.shape == (0,0) - assert B.shape == (0,2) - assert C.shape == (3,0) - np.testing.assert_array_almost_equal(D, Dc) - - n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - assert A.shape == (0,0) - assert B.shape == (0,2) - assert C.shape == (3,0) - np.testing.assert_array_almost_equal(D, Dr) - -def test_td04ad_static(): - """Regression: td04ad (TFM -> SS transformation) for static TFM""" - from itertools import product - for nout, nin, rc in product(range(1, 6), range(1, 6), ['R', 'C']): - Dref = np.zeros((nout, nin)) - if rc == 'R': - num = np.reshape(np.arange(nout * nin), (nout, nin, 1)) - den = np.reshape(np.arange(1, 1 + nout), (nout, 1)) - index = np.repeat(0, nout) - Dref = num[:nout, :nin, 0] / np.broadcast_to(den, (nout, nin)) - else: - maxn = max(nout, nin) - num = np.zeros((maxn, maxn, 1)) - num[:nout, :nin, 0] = np.reshape( - np.arange(nout * nin), (nout, nin)) - den = np.reshape(np.arange(1, 1 + nin), (nin, 1)) - index = np.repeat(0, nin) - Dref = num[:nout, :nin, 0] / np.broadcast_to(den.T, (nout, nin)) - nr, A, B, C, D = transform.td04ad(rc, nin, nout, index, den, num) - np.testing.assert_equal(nr, 0) - for M in [A, B, C]: - np.testing.assert_equal(M, np.zeros_like(M)) - np.testing.assert_almost_equal(D, Dref) - -def test_mixfeedthrough(): - """Test case popping up from control testing - - a mix of feedthrough and dynamics. The problem from the control - package was somewhere else - """ - num = np.array([ [ [ 0.0, 0.0 ], [ 0.0, -0.2 ] ], - [ [ -0.1, 0.0 ], [ 0.0, 0.0 ] ] ]) - p, m, d = num.shape - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - denc = np.array([[1.0, 1.1], - [1.0, 0.0]]) - idxc = np.array([1, 0]) - n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc) - np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]])) - -def test_toandfrom(): - A = np.array([[-3.0]]) - B = np.array([[0.1, 0.0]]) - C = np.array([[1.0], - [0.0]]) - D = np.array([[0.0, 0.0], - [0.0, 1.0]]) - - tfout = transform.tb04ad(1, 2, 2, A, B, C, D) - - num = tfout[6] - den = tfout[5] - idxc = np.array([1, 0]) - n, At, Bt, Ct, Dt = transform.td04ad('R', 2, 2, idxc, den, num) - np.testing.assert_array_almost_equal(D, Dt) - np.testing.assert_array_almost_equal(A, At) - -def test_tfm2ss_6(): - """Python version of Fortran test program from - -- Bug in TD04AD when ROWCOL='C' #6 - This bug was fixed in PR #27""" - m = 1 - p = 1 - index = np.array([0]) - dcoeff = np.array([[0.5]]) - ucoeff = np.array([[[32]]]) - n, A, B, C, D = transform.td04ad('R', m, p, index, dcoeff, ucoeff) - assert n == 0 - np.testing.assert_array_almost_equal(D, np.array([[64]])) - n, A, B, C, D = transform.td04ad('C', m, p, index, dcoeff, ucoeff) - assert n == 0 - np.testing.assert_array_almost_equal(D, np.array([[64]])) +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError, SlycotParameterError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb04ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.td04ad, SlycotArithmeticError, (3,), {}),)) + def test_td04ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + + def test_td04ad_c(self): + """td04ad: Convert with 'C' option""" + + # for octave: + """ + num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ]; + [3.0, -1.0, 1.0 ], [ 0.0, 1.0 ]; + [0.0, 0.0, 1.0], [ 0.0, 2.0 ] }; + den = { [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; + [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; + [1.0, 0.4, 3.0], [ 1.0, 1.0 ]}; + """ + + m = 2 + p = 3 + d = 3 + num = np.array([ + [ [0.0, 0.0, 1.0], [1.0, 0.0, 0.0] ], + [ [3.0, -1.0, 1.0], [0.0, 1.0, 0.0] ], + [ [0.0, 0.0, 1.0], [0.0, 2.0, 0.0] ] ]) + + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + denc = np.array( + [ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ]) + indc = np.array( + [ 2, 1 ], dtype=int) + + nref = 3 + Aref = np.array([ [-1, 0, 0], + [ 0, -0.4, -0.3], + [ 0, 10, 0] ]) + Bref = np.array([ [0, -1], + [1, 0], + [0, 0] ]) + Cref = np.array([ [1, 0, 0.1], + [-1, -2.2, -0.8], + [-2, 0, 0.1] ]) + Dref = np.array([ [0, 1], + [3, 0], + [0, 0] ]) + + nr, A, B, C, D = transform.td04ad('C', m, p, indc, denc, numc) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + np.testing.assert_equal(nref, nr) + # the returned state space representation is not guaranteed to + # be of one form for all architectures, so we transform back + # to tf and check for equality then + _, _, _, _, _, dcoeff, ucoeff = transform.tb04ad( + nr, m, p, A, B, C, D) + _, _, _, _, _, dcoeffref, ucoeffref = transform.tb04ad( + nref, m, p, Aref, Bref, Cref, Dref) + np.testing.assert_array_almost_equal(dcoeff,dcoeffref) + np.testing.assert_array_almost_equal(ucoeff,ucoeffref) + + def test_td04ad_r(self): + """td04ad: Convert with 'R' option + + example program from + http://slicot.org/objects/software/shared/doc/TD04AD.html + """ + + m = 2 + p = 2 + rowcol = 'R' + index = [3, 3] + dcoeff = np.array([ [1.0, 6.0, 11.0, 6.0], [1.0, 6.0, 11.0, 6.0] ]) + + ucoeff = np.array([ [[1.0, 6.0, 12.0, 7.0], [0.0, 1.0, 4.0, 3.0]], + [[0.0, 0.0, 1.0, 1.0], [1.0, 8.0, 20.0, 15.0]] ]) + + nref = 3 + + Aref = np.array([ [ 0.5000, -0.8028, 0.9387], + [ 4.4047, -2.3380, 2.5076], + [-5.5541, 1.6872, -4.1620] ]) + Bref = np.array([ [-0.2000, -1.2500], + [ 0.0000, -0.6097], + [ 0.0000, 2.2217] ]) + Cref = np.array([ [0.0000, -0.8679, 0.2119], + [0.0000, 0.0000, 0.9002] ]) + Dref = np.array([ [1.0000, 0.0000], + [0.0000, 1.0000] ]) + + nr, A, B, C, D = transform.td04ad(rowcol, m, p, index, dcoeff, ucoeff) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + np.testing.assert_equal(nref, nr) + # order of states is not guaranteed, so we reorder the reference + rindex = np.flip(np.argsort(np.diag(A))) + Arref = Aref[rindex, :][:, rindex] + Brref = Bref[rindex, :] + Crref = Cref[:, rindex] + Drref = Dref + np.testing.assert_array_almost_equal(A, Arref,decimal=4) + np.testing.assert_array_almost_equal(B, Brref,decimal=4) + np.testing.assert_array_almost_equal(C, Crref,decimal=4) + np.testing.assert_array_almost_equal(D, Drref,decimal=4) + + + def test_staticgain(self): + """td04ad: Convert a transferfunction to SS with only static gain""" + + # 2 inputs, 3 outputs? columns share a denominator + num = np.array([ [ [1.0], [2.0] ], + [ [0.2], [4.3] ], + [ [1.2], [3.2] ] ]) + p, m, d = num.shape + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + + # denc, columns share a common denominator + denc = np.array([ [ 1.0], [0.5] ]) + Dc = (num / denc).reshape((3,2)) + idxc = np.zeros((2,), dtype=int) + + # denr, rows share a common denominator + denr = np.array([ [1.0], [0.5], [3.0] ]) + idxr = np.zeros((3,), dtype=int) + Dr = (num / denr[:, np.newaxis]).reshape((3,2)) + + # fails with: + # On entry to TB01XD parameter number 5 had an illegal value + + n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + assert A.shape == (0,0) + assert B.shape == (0,2) + assert C.shape == (3,0) + np.testing.assert_array_almost_equal(D, Dc) + + n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + assert A.shape == (0,0) + assert B.shape == (0,2) + assert C.shape == (3,0) + np.testing.assert_array_almost_equal(D, Dr) + + def test_td04ad_static(self): + """Regression: td04ad (TFM -> SS transformation) for static TFM""" + from itertools import product + for nout, nin, rc in product(range(1, 6), range(1, 6), ['R', 'C']): + Dref = np.zeros((nout, nin)) + if rc == 'R': + num = np.reshape(np.arange(nout * nin), (nout, nin, 1)) + den = np.reshape(np.arange(1, 1 + nout), (nout, 1)) + index = np.repeat(0, nout) + Dref = num[:nout, :nin, 0] / np.broadcast_to(den, (nout, nin)) + else: + maxn = max(nout, nin) + num = np.zeros((maxn, maxn, 1)) + num[:nout, :nin, 0] = np.reshape( + np.arange(nout * nin), (nout, nin)) + den = np.reshape(np.arange(1, 1 + nin), (nin, 1)) + index = np.repeat(0, nin) + Dref = num[:nout, :nin, 0] / np.broadcast_to(den.T, (nout, nin)) + nr, A, B, C, D = transform.td04ad(rc, nin, nout, index, den, num) + np.testing.assert_equal(nr, 0) + for M in [A, B, C]: + np.testing.assert_equal(M, np.zeros_like(M)) + np.testing.assert_almost_equal(D, Dref) + + def test_mixfeedthrough(self): + """Test case popping up from control testing + + a mix of feedthrough and dynamics. The problem from the control + package was somewhere else + """ + num = np.array([ [ [ 0.0, 0.0 ], [ 0.0, -0.2 ] ], + [ [ -0.1, 0.0 ], [ 0.0, 0.0 ] ] ]) + p, m, d = num.shape + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + denc = np.array([[1.0, 1.1], + [1.0, 0.0]]) + idxc = np.array([1, 0]) + n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc) + np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]])) + + def test_toandfrom(self): + A = np.array([[-3.0]]) + B = np.array([[0.1, 0.0]]) + C = np.array([[1.0], + [0.0]]) + D = np.array([[0.0, 0.0], + [0.0, 1.0]]) + + tfout = transform.tb04ad(1, 2, 2, A, B, C, D) + + num = tfout[6] + den = tfout[5] + idxc = np.array([1, 0]) + n, At, Bt, Ct, Dt = transform.td04ad('R', 2, 2, idxc, den, num) + np.testing.assert_array_almost_equal(D, Dt) + np.testing.assert_array_almost_equal(A, At) + + def test_tfm2ss_6(self): + """Python version of Fortran test program from + -- Bug in TD04AD when ROWCOL='C' #6 + This bug was fixed in PR #27""" + m = 1 + p = 1 + index = np.array([0]) + dcoeff = np.array([[0.5]]) + ucoeff = np.array([[[32]]]) + n, A, B, C, D = transform.td04ad('R', m, p, index, dcoeff, ucoeff) + assert n == 0 + np.testing.assert_array_almost_equal(D, np.array([[64]])) + n, A, B, C, D = transform.td04ad('C', m, p, index, dcoeff, ucoeff) + assert n == 0 + np.testing.assert_array_almost_equal(D, np.array([[64]])) diff --git a/slycot/tests/test_tf01md.py b/slycot/tests/test_tf01md.py new file mode 100644 index 00000000..d349f5f5 --- /dev/null +++ b/slycot/tests/test_tf01md.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tf01md: + + @pytest.mark.skip(reason="not implemented") + def test_tf01md(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tf01rd.py b/slycot/tests/test_tf01rd.py new file mode 100644 index 00000000..f67226eb --- /dev/null +++ b/slycot/tests/test_tf01rd.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tf01rd: + + @pytest.mark.skip(reason="not implemented") + def test_tf01rd(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tg01ad.py b/slycot/tests/test_tg01ad.py index 864d41be..64f4b84c 100644 --- a/slycot/tests/test_tg01ad.py +++ b/slycot/tests/test_tg01ad.py @@ -1,6 +1,3 @@ -# =================================================== -# tg01ad tests - import numpy as np from numpy.testing import assert_almost_equal, assert_equal, assert_raises diff --git a/slycot/tests/test_tg01fd.py b/slycot/tests/test_tg01fd.py index 902d1ef2..1168cc80 100644 --- a/slycot/tests/test_tg01fd.py +++ b/slycot/tests/test_tg01fd.py @@ -1,6 +1,3 @@ -# =================================================== -# tg01fd tests - import numpy as np from numpy.testing import assert_almost_equal, assert_equal, assert_raises diff --git a/slycot/tests/test_transform.py b/slycot/tests/test_transform.py deleted file mode 100644 index 928ac5d4..00000000 --- a/slycot/tests/test_transform.py +++ /dev/null @@ -1,22 +0,0 @@ -# -# test_transform.py - generic tests for transform routines -# repagh 0. + m : int + The number of columns of matrix B. It represents the dimension of + the input vector. m > 0. + p : int + The number of rows of matrix C. It represents the dimension of + the output vector. p > 0. + maxred : float + The maximum allowed reduction in the 1-norm of S (in an iteration) + if zero rows or columns are encountered. + If maxred > 0.0, maxred must be larger than one (to enable the norm + reduction). + If maxred <= 0.0, then the value 10.0 for maxred is used. + A : (n, n) array_like + The leading n-by-n part of this array must contain the system state + matrix A. + B : (n, m) array_like + The leading n-by-m part of this array must contain the system input + matrix B. + C : (p, n) array_like + The leading p-by-n part of this array must contain the system output + matrix C. + job := {'A', 'B', 'C', 'N'}, optional + Indicates which matrices are involved in balancing, as follows: + = 'A': All matrices are involved in balancing; + = 'B': B and A matrices are involved in balancing; + = 'C': C and A matrices are involved in balancing; + = 'N': B and C matrices are not involved in balancing. + Returns + ------- + s_norm : float + The 1-norm of the given matrix S is non-zero, the ratio between + the 1-norm of the given matrix and the 1-norm of the balanced matrix. + A : (n, n) ndarray + The leading n-by-n part of this array contains the balanced matrix + inv(D)*A*D. + B : (n, m) ndarray + The leading n-by-m part of this array contains the balanced matrix + inv(D)*B. + C : (p ,n) ndarray + The leading p-by-n part of this array contains the balanced matrix C*D. + scale : rank-1 array('d') with bounds (n) + The scaling factors applied to S. If D(j) is the scaling factor + applied to row and column j, then scale(j) = D(j), for j = 1,...,n. - Required arguments: - n : input int - The order of the matrix A, the number of rows of matrix B and - the number of columns of matrix C. It represents the dimension of - the state vector. n > 0. - m : input int - The number of columns of matrix B. It represents the dimension of - the input vector. m > 0. - p : input int - The number of rows of matrix C. It represents the dimension of - the output vector. p > 0. - maxred : input float - The maximum allowed reduction in the 1-norm of S (in an iteration) - if zero rows or columns are encountered. - If maxred > 0.0, maxred must be larger than one (to enable the norm - reduction). - If maxred <= 0.0, then the value 10.0 for maxred is used. - A : input rank-2 array('d') with bounds (n,n) - The leading n-by-n part of this array must contain the system state - matrix A. - B : input rank-2 array('d') with bounds (n,m) - The leading n-by-m part of this array must contain the system input - matrix B. - C : input rank-2 array('d') with bounds (p,n) - The leading p-by-n part of this array must contain the system output - matrix C. - Optional arguments: - job := 'A' input string(len=1) - Indicates which matrices are involved in balancing, as follows: - = 'A': All matrices are involved in balancing; - = 'B': B and A matrices are involved in balancing; - = 'C': C and A matrices are involved in balancing; - = 'N': B and C matrices are not involved in balancing. - Return objects: - s_norm : float - The 1-norm of the given matrix S is non-zero, the ratio between - the 1-norm of the given matrix and the 1-norm of the balanced matrix. - A : rank-2 array('d') with bounds (n,n) - The leading n-by-n part of this array contains the balanced matrix - inv(D)*A*D. - B : rank-2 array('d') with bounds (n,m) - The leading n-by-m part of this array contains the balanced matrix - inv(D)*B. - C : rank-2 array('d') with bounds (p,n) - The leading p-by-n part of this array contains the balanced matrix C*D. - scale : rank-1 array('d') with bounds (n) - The scaling factors applied to S. If D(j) is the scaling factor - applied to row and column j, then scale(j) = D(j), for j = 1,...,n. + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value; """ hidden = ' (hidden by the wrapper)' arg_list = ['job', 'N', 'M', 'P', 'maxred', 'A', 'LDA'+hidden, 'B', @@ -100,6 +105,93 @@ def tb01id(n,m,p,maxred,a,b,c,job='A'): raise_if_slycot_error(out[-1], arg_list) return out[:-1] +def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None): + """Ar, Br, Cr, nr = tb01pd(n,m,p,A,B,C,[job,equil,tol,ldwork]) + + To find a reduced (controllable, observable, or minimal) state- + space representation (Ar,Br,Cr) for any original state-space + representation (A,B,C). The matrix Ar is in upper block + Hessenberg form. + + Parameters + ---------- + n : int + Order of the State-space representation. + m : int + Number of inputs. + p : int + Number of outputs. + A : (n, n) array_like + State dynamics matrix. + B : (n, max(m,p)) array_like + The leading n-by-m part of this array must contain the original + input/state matrix B; the remainder of the leading n-by-max(m,p) + part is used as internal workspace. + C : (p, n) array_like + The leading p-by-n part of this array must contain the original + state/output matrix C; the remainder of the leading max(1,m,p)-by-n + part is used as internal workspace. + job : {'M', 'C', 'O'}, optional + Indicates whether the user wishes to remove the + uncontrollable and/or unobservable parts as follows: + = 'M': Remove both the uncontrollable and unobservable + parts to get a minimal state-space representation; + = 'C': Remove the uncontrollable part only to get a + controllable state-space representation; + = 'O': Remove the unobservable part only to get an + observable state-space representation. + Default is 'M'. + equil : {'S', 'N'}, optional + Specifies whether the user wishes to preliminarily balance + the triplet (A,B,C) as follows: + = 'S': Perform balancing (scaling); + = 'N': Do not perform balancing. + + Returns + ------- + Ar : (nr, nr) ndarray + Contains the upper block Hessenberg state dynamics matrix + Ar of a minimal, controllable, or observable realization + for the original system, depending on the value of JOB, + JOB = 'M', JOB = 'C', or JOB = 'O', respectively. + Br : (nr, m) ndarray + Contains the transformed input/state matrix Br of a + minimal, controllable, or observable realization for the + original system, depending on the value of JOB, JOB = 'M', + JOB = 'C', or JOB = 'O', respectively. If JOB = 'C', only + the first IWORK(1) rows of B are nonzero. + Cr : (p, nr) ndarray + Contains the transformed state/output matrix Cr of a + minimal, C controllable, or observable realization for the + original C system, depending on the value of JOB, JOB = + 'M', C JOB = 'C', or JOB = 'O', respectively. C If JOB = + 'M', or JOB = 'O', only the last IWORK(1) columns C (in + the first NR columns) of C are nonzero. + nr : int + The order of the reduced state-space representation + (Ar,Br,Cr) of a minimal, controllable, or observable + realization for the original system, depending on + JOB = 'M', JOB = 'C', or JOB = 'O'. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value + """ + hidden = ' (hidden by the wrapper)' + arg_list = ['job', 'equil', 'n','m','p','A','lda'+hidden,'B','ldb'+hidden, + 'C','ldc'+hidden,'nr','tol','iwork'+hidden,'dwork'+hidden, + 'ldwork','info'+hidden] + if ldwork is None: + ldwork = max(1, n+max(n,3*m,3*p)) + elif ldwork < max(1, n+max(n,3*m,3*p)): + raise SlycotParameterError("ldwork is too small", -15) + out = _wrapper.tb01pd(n=n,m=m,p=p,a=A,b=B,c=C, + job=job,equil=equil,tol=tol,ldwork=ldwork) + + raise_if_slycot_error(out[-1], arg_list) + return out[:-1] + def tb03ad(n,m,p,A,B,C,D,leri,equil='N',tol=0.0,ldwork=None): """ A_min,b_min,C_min,nr,index,pcoeff,qcoeff,vcoeff = tb03ad_l(n,m,p,A,B,C,D,leri,[equil,tol,ldwork]) @@ -116,92 +208,102 @@ def tb03ad(n,m,p,A,B,C,D,leri,equil='N',tol=0.0,ldwork=None): Additionally a minimal realization (A_min,B_min,C_min) of the original system (A,B,C) is returned. - Required arguments: - n : input int - The order of the state-space representation, i.e. the order of - the original state dynamics matrix A. n > 0. - m : input int - The number of system inputs. m > 0. - p : input int - The number of system outputs. p > 0. - A : input rank-2 array('d') with bounds (n,n) - The leading n-by-n part of this array must contain the original - state dynamics matrix A. - B : input rank-2 array('d') with bounds (n,max(m,p)) - The leading n-by-m part of this array must contain the original - input/state matrix B; the remainder of the leading n-by-max(m,p) - part is used as internal workspace. - C : input rank-2 array('d') with bounds (max(m,p),n) - The leading p-by-n part of this array must contain the original - state/output matrix C; the remainder of the leading max(m,p)-by-n - part is used as internal workspace. - D : input rank-2 array('d') with bounds (max(m,p),max(m,p)) - The leading p-by-m part of this array must contain the original - direct transmission matrix D; the remainder of the leading - max(m,p)-by-max(m,p) part is used as internal workspace. - leri : input string(len=1) - Indicates whether the left polynomial matrix representation or - the right polynomial matrix representation is required. - Optional arguments: - equil := 'N' input string(len=1) - Specifies whether the user wishes to balance the triplet (A,B,C), - before computing a minimal state-space representation, as follows: - = 'S': Perform balancing (scaling); - = 'N': Do not perform balancing. - tol := 0.0 input float - The tolerance to be used in rank determination when transforming - (A, B). If tol <= 0 a default value is used. - ldwork := max(2*n+3*max(m,p), p*(p+2)) input int - The length of the cache array. - ldwork >= max( n + max(n, 3*m, 3*p), pm*(pm + 2)) - where pm = p, if leri = 'L'; - pm = m, if leri = 'R'. - For optimum performance it should be larger. - Return objects: - A_min : rank-2 array('d') with bounds (n,n) - The leading nr-by-nr part of this array contains the upper block - Hessenberg state dynamics matrix A_min of a minimal realization for - the original system. - B_min : rank-2 array('d') with bounds (n,max(m,p)) - The leading nr-by-m part of this array contains the transformed - input/state matrix B_min. - C_min : rank-2 array('d') with bounds (max(m,p),n) - The leading p-by-nr part of this array contains the transformed - state/output matrix C_min. - nr : int - The order of the minimal state-space representation - (A_min,B_min,C_min). - index : rank-1 array('i') with bounds either (p) or (m) - If leri = 'L', index(i), i = 1,2,...,p, contains the maximum degree - of the polynomials in the i-th row of the denominator matrix P(s) - of the left polynomial matrix representation. These elements are - ordered so that index(1) >= index(2) >= ... >= index(p). - If leri = 'R', index(i), i = 1,2,...,m, contains the maximum degree - of the polynomials in the i-th column of the denominator matrix P(s) - of the right polynomial matrix representation. These elements are - ordered so that index(1) >= index(2) >= ... >= index(m). - pcoeff : rank-3 array('d') with bounds either (p,p,n+1) or (m,m,n+1) - If leri = 'L' then porm = p, otherwise porm = m. - The leading porm-by-porm-by-kpcoef part of this array contains - the coefficients of the denominator matrix P(s), where - kpcoef = max(index) + 1. - pcoeff(i,j,k) is the coefficient in s**(index(iorj)-k+1) of - polynomial (i,j) of P(s), where k = 1,2,...,kpcoef; if leri = 'L' - then iorj = I, otherwise iorj = J. Thus for leri = 'L', - P(s) = diag(s**index)*(pcoeff(.,.,1)+pcoeff(.,.,2)/s+...). - qcoeff : rank-3 array('d') with bounds (p,m,n + 1) or (max(m,p),max(m,p)) - If leri = 'L' then porp = m, otherwise porp = p. - If leri = 'L', the leading porm-by-porp-by-kpcoef part of this array - contains the coefficients of the numerator matrix Q(s). - If leri = 'R', the leading porp-by-porm-by-kpcoef part of this array - contains the coefficients of the numerator matrix Q(s). - qcoeff(i,j,k) is defined as for pcoeff(i,j,k). - vcoeff : rank-3 array('d') with bounds (p,n,n+1) or (m,n,n+1) - The leading porm-by-nr-by-kpcoef part of this array contains - the coefficients of the intermediate matrix V(s). - vcoeff(i,j,k) is defined as for pcoeff(i,j,k). + Parameters + ---------- + n : int + The order of the state-space representation, i.e. the order of + the original state dynamics matrix A. n > 0. + m : int + The number of system inputs. m > 0. + p : int + The number of system outputs. p > 0. + A : (n, n) array_like + The leading n-by-n part of this array must contain the original + state dynamics matrix A. + B : (n, max(m,p)) array_like + The leading n-by-m part of this array must contain the original + input/state matrix B; the remainder of the leading n-by-max(m,p) + part is used as internal workspace. + C : (max(m,p), n) + The leading p-by-n part of this array must contain the original + state/output matrix C; the remainder of the leading max(m,p)-by-n + part is used as internal workspace. + D : (max(m,p), max(m,p)) array_like + The leading p-by-m part of this array must contain the original + direct transmission matrix D; the remainder of the leading + max(m,p)-by-max(m,p) part is used as internal workspace. + leri : {'L', 'R'} + Indicates whether the left polynomial matrix representation or + the right polynomial matrix representation is required. + = 'L': A left matrix fraction is required; + = 'R': A right matrix fraction is required. + equil : {'S', 'N'}, optional + Specifies whether the user wishes to balance the triplet (A,B,C), + before computing a minimal state-space representation, as follows: + = 'S': Perform balancing (scaling); + = 'N': Do not perform balancing. + Default is `N`. + tol : float, optional + The tolerance to be used in rank determination when transforming + (A, B). If tol <= 0 a default value is used. + Default is `0.0`. + ldwork : int, optional + The length of the cache array. + ldwork >= max( n + max(n, 3*m, 3*p), pm*(pm + 2)) + where pm = p, if leri = 'L'; + pm = m, if leri = 'R'. + For optimum performance it should be larger. + Default is None. + + Returns + ------- + A_min : (n, n) ndarray + The leading nr-by-nr part of this array contains the upper block + Hessenberg state dynamics matrix A_min of a minimal realization for + the original system. + B_min : (n, max(m,p)) ndarray + The leading nr-by-m part of this array contains the transformed + input/state matrix B_min. + C_min : (max(m,p), n) ndarray + The leading p-by-nr part of this array contains the transformed + state/output matrix C_min. + nr : int + The order of the minimal state-space representation + (A_min,B_min,C_min). + index : (p, ) or (m, ) ndarray + If leri = 'L', index(i), i = 1,2,...,p, contains the maximum degree + of the polynomials in the i-th row of the denominator matrix P(s) + of the left polynomial matrix representation. These elements are + ordered so that index(1) >= index(2) >= ... >= index(p). + If leri = 'R', index(i), i = 1,2,...,m, contains the maximum degree + of the polynomials in the i-th column of the denominator matrix P(s) + of the right polynomial matrix representation. These elements are + ordered so that index(1) >= index(2) >= ... >= index(m). + pcoeff : (p, p, n+1) or (m, m, n+1) ndarray + If leri = 'L' then porm = p, otherwise porm = m. + The leading porm-by-porm-by-kpcoef part of this array contains + the coefficients of the denominator matrix P(s), where + kpcoef = max(index) + 1. + pcoeff(i,j,k) is the coefficient in s**(index(iorj)-k+1) of + polynomial (i,j) of P(s), where k = 1,2,...,kpcoef; if leri = 'L' + then iorj = I, otherwise iorj = J. Thus for leri = 'L', + P(s) = diag(s**index)*(pcoeff(.,.,1)+pcoeff(.,.,2)/s+...). + qcoeff : (p, m, n+1) or (max(m,p), max(m,p), n+1) ndarray + If leri = 'L' then porp = m, otherwise porp = p. + If leri = 'L', the leading porm-by-porp-by-kpcoef part of this array + contains the coefficients of the numerator matrix Q(s). + If leri = 'R', the leading porp-by-porm-by-kpcoef part of this array + contains the coefficients of the numerator matrix Q(s). + qcoeff(i,j,k) is defined as for pcoeff(i,j,k). + vcoeff : (p, n, n+1) or (m, n, n+1) ndarray + The leading porm-by-nr-by-kpcoef part of this array contains + the coefficients of the intermediate matrix V(s). + vcoeff(i,j,k) is defined as for pcoeff(i,j,k). + Raises ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value; SlycotArithmeticError :info == 1: A singular matrix was encountered during the @@ -209,7 +311,6 @@ def tb03ad(n,m,p,A,B,C,D,leri,equil='N',tol=0.0,ldwork=None): :info == 2: A singular matrix was encountered during the computation of P(s). - """ hidden = ' (hidden by the wrapper)' arg_list = ['leri', 'equil', 'n', 'm', 'P', 'A', 'LDA'+hidden, 'B', @@ -233,58 +334,59 @@ def tb03ad(n,m,p,A,B,C,D,leri,equil='N',tol=0.0,ldwork=None): def tb04ad(n,m,p,A,B,C,D,tol1=0.0,tol2=0.0,ldwork=None): """ Ar,Br,Cr,nr,denom_degs,denom_coeffs,num_coeffs = tb04ad(n,m,p,A,B,C,D,[tol1,tol2,ldwork]) - Convert a state-space system to a tranfer function or matrix of transfer functions. + Convert a state-space system to a transfer function or matrix of transfer functions. The transfer function is given as rows over common denominators. - Required arguments - ------------------ - - n : integer - state dimension - m : integer - input dimension - p : integer - output dimension - A : rank-2 array, shape(n,n) - state dynamics matrix. - B : rank-2 array, shape (n,m) - input matrix - C : rank-2 array, shape (p,n) - output matri - D : rank-2 array, shape (p,m) - direct transmission matrix - - Optional arguments - ------------------ - - tol1 = 0.0: double - tolerance in determining the transfer function coefficients, - when set to 0, a default value is used - tol2 = 0.0: double - tolerance in separating out a controllable/observable subsystem - of (A,B,C), when set to 0, a default value is used - ldwork : int - The length of the cache array. The default values is - max(1,n*(n+1)+max(n*m+2*n+max(n,p),max(3*m,p))) - - Returns - ------- - - nr : int - state dimension of the controllable subsystem - Ar : rank-2 array, shape(nr,nr) - state dynamics matrix of the controllable subsystem - Br : rank-2 array, shape (nr,m) - input matrix of the controllable subsystem - Cr : rank-2 array, shape (p,nr) - output matri of the controllable subsystem - index : rank-1 array, shape (p) - array of orders of the denominator polynomials - dcoeff : rank-2 array, shape (p,max(index)+1) - array of denominator coefficients - ucoeff : rank-3 array, shape (p,m,max(index)+1) - array of numerator coefficients + Parameters + ---------- + n : int + state dimension + m : int + input dimension + p : int + output dimension + A : (n, n) array_like + state dynamics matrix. + B : (n, m) array_like + input matrix + C : (p, n) array_like + output matrix + D : (p, m) array_like + direct transmission matrix + tol1 : float, optional + tolerance in determining the transfer function coefficients, + when set to 0, a default value is used + Default is `0.0`. + tol2 : float, optional + tolerance in separating out a controllable/observable subsystem + of (A,B,C), when set to 0, a default value is used + Default is `0.0`. + ldwork : int, optional + The length of the cache array. The default values is + max(1,n*(n+1)+max(n*m+2*n+max(n,p),max(3*m,p))) + Default is None. + + Returns + ------- + nr : int + state dimension of the controllable subsystem + Ar : (nr, nr) ndarray + state dynamics matrix of the controllable subsystem + Br : (nr, m) ndarray + input matrix of the controllable subsystem + Cr : (p, nr) ndarray + output matrix of the controllable subsystem + index : (p, ) ndarray + array of orders of the denominator polynomials + dcoeff : (p, max(index)+1) ndarray + array of denominator coefficients + ucoeff : (p, m, max(index)+1) ndarray + array of numerator coefficients + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)' arg_list = ['rowcol','n','m','p','A','lda'+hidden,'B','ldb'+hidden,'C','ldc'+hidden,'D', 'ldd'+hidden, @@ -340,7 +442,7 @@ def tb05ad(n, m, p, jomega, A, B, C, job='NG'): m : int The number of inputs, i.e. the number of columns in the matrix B. - p : int + p : int The number of outputs, i.e. the number of rows in the matrix C. jomega : complex float @@ -349,14 +451,14 @@ def tb05ad(n, m, p, jomega, A, B, C, job='NG'): systems, this is j*omega, where omega is the frequency to be evaluated. For discrete time systems, freq = exp(j*omega*Ts) - A : (n,n) ndarray + A : (n, n) ndarray On entry, this array must contain the state transition matrix A. - B : (n,m) ndarray + B : (n, m) ndarray On entry, this array must contain the input/state matrix B. - C : (p,n) ndarray + C : (p, n) ndarray On entry, of this array must contain the state/output matrix C. - job : {'AG', 'NG', 'NH'} + job : {'AG', 'NG', 'NH'}, optional If job = 'AG' (i.e., 'all', 'general matrix'), the A matrix is first balanced. The balancing transformation is then appropriately applied to matrices B and C. The A matrix @@ -377,75 +479,66 @@ def tb05ad(n, m, p, jomega, A, B, C, job='NG'): Returns ------- - if job = 'AG': - -------------- - At: The A matrix which has been both balanced and - transformed to upper Hessenberg form. The balancing - transforms A according to - A1 = P^-1 * A * P. - The transformation to upper Hessenberg form then yields - At = Q^T * (P^-1 * A * P ) * Q. - Note that the lower triangle of At is in general not zero. - Rather, it contains information on the orthogonal matrix Q - used to transform A1 to Hessenberg form. See docs for lappack - DGEHRD(): - http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html - However, it does not apparently contain information on P, the - matrix used in the balancing procedure. - - Bt: The matrix B transformed according to - Bt = Q^T * P^-1 * B. - - Ct: The matrix C transformed according to - Ct = C * P * Q - - rcond: RCOND contains an estimate of the reciprocal of the - condition number of matrix H with respect to inversion, where - H = (j*freq * I - A) - - g_jw: complex p-by-m array, which contains the frequency response - matrix G(freq). - - ev: Eigenvalues of the matrix A. - - hinvb : complex n-by-m array, which contains the product - -1 - H B. - - if job = 'NG': - -------------- - At: The matrix A transformed to upper Hessenberg form according - to - At = Q^T * A * Q. - The lower triangle is not zero. It containts info on the - orthoganal transformation. See docs for linpack DGEHRD() - http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html - - Bt: The matrix B transformed according to - Bt = Q^T * B. - - Ct: The matrix C transformed according to - Ct = C * Q - g_jw: complex array with dim p-by-m which contains the frequency - response matrix G(freq). - - hinvb : complex array with dimension p-by-m. - This array contains the - -1 - product H B. - + if job = 'AG' + ------------- + At : The A matrix which has been both balanced and + transformed to upper Hessenberg form. The balancing + transforms A according to + A1 = P^-1 * A * P. + The transformation to upper Hessenberg form then yields + At = Q^T * (P^-1 * A * P ) * Q. + Note that the lower triangle of At is in general not zero. + Rather, it contains information on the orthogonal matrix Q + used to transform A1 to Hessenberg form. See docs for lappack + DGEHRD(): + http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html + However, it does not apparently contain information on P, the + matrix used in the balancing procedure. + Bt : The matrix B transformed according to + Bt = Q^T * P^-1 * B. + Ct : The matrix C transformed according to + Ct = C * P * Q + rcond : RCOND contains an estimate of the reciprocal of the + condition number of matrix H with respect to inversion, where + H = (j*freq * I - A) + g_jw : complex p-by-m array, which contains the frequency response + matrix G(freq). + ev : Eigenvalues of the matrix A. + hinvb : complex n-by-m array, which contains the product + -1 + H B. + + if job = 'NG' + ------------- + At : The matrix A transformed to upper Hessenberg form according + to + At = Q^T * A * Q. + The lower triangle is not zero. It containts info on the + orthoganal transformation. See docs for linpack DGEHRD() + http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html + Bt : The matrix B transformed according to + Bt = Q^T * B. + Ct : The matrix C transformed according to + Ct = C * Q + g_jw : complex array with dim p-by-m which contains the frequency + response matrix G(freq). + hinvb : complex array with dimension p-by-m. + This array contains the + -1 + product H B. if job = 'NH' - -------------- - g_jw: complex p-by-m array which contains the frequency - response matrix G(freq). - - hinvb : complex p-by-m array which contains the - -1 - product H B. + ------------- + g_jw : complex p-by-m array which contains the frequency + response matrix G(freq). + hinvb : complex p-by-m array which contains the + -1 + product H B. Raises ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value SlycotArithmeticError :info = 1: More than {n30} (30*`n`) iterations were required to isolate the @@ -573,6 +666,8 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None): Raises ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value; SlycotArithmeticError :info > 0: i={info} is the first index of `dcoeff` for which @@ -644,67 +739,74 @@ def tc04ad(m,p,index,pcoeff,qcoeff,leri,ldwork=None): respectively. - Required arguments: - m : input int - The number of system inputs. m > 0. - p := len(index) input int - The number of system outputs. p > 0. - index : input rank-1 array('i') with bounds (p) or (m) - If leri = 'L', index(i), i = 1,2,...,p, must contain the maximum - degree of the polynomials in the I-th row of the denominator matrix - P(s) of the given left polynomial matrix representation. - If leri = 'R', index(i), i = 1,2,...,m, must contain the maximum - degree of the polynomials in the I-th column of the denominator - matrix P(s) of the given right polynomial matrix representation. - pcoeff : input rank-3 array('d') with bounds (p,p,*) or (m,m,*) - If leri = 'L' then porm = p, otherwise porm = m. The leading - porm-by-porm-by-kpcoef part of this array must contain - the coefficients of the denominator matrix P(s). pcoeff(i,j,k) is - the coefficient in s**(index(iorj)-K+1) of polynomial (I,J) of P(s), - where k = 1,2,...,kpcoef and kpcoef = max(index) + 1; if leri = 'L' - then iorj = i, otherwise iorj = j. Thus for leri = 'L', - P(s) = diag(s**index)*(pcoeff(.,.,1)+pcoeff(.,.,2)/s+...). - If leri = 'R', pcoeff is modified by the routine but restored on exit. - qcoeff : input rank-3 array('d') with bounds (p,m,*) or (max(m,p),max(m,p),*) - If leri = 'L' then porp = m, otherwise porp = p. The leading - porm-by-porp-by-kpcoef part of this array must contain - the coefficients of the numerator matrix Q(s). - qcoeff(i,j,k) is defined as for pcoeff(i,j,k). - If leri = 'R', qcoeff is modified by the routine but restored on exit. - leri : input string(len=1) - Indicates whether a left polynomial matrix representation or a right - polynomial matrix representation is input as follows: - = 'L': A left matrix fraction is input; - = 'R': A right matrix fraction is input. - Optional arguments: - ldwork := max(m,p)*(max(m,p)+4) input int - The length of the cache array. ldwork >= max(m,p)*(max(m,p)+4) - For optimum performance it should be larger. - Return objects: - n : int - The order of the resulting state-space representation. - That is, n = sum(index). - rcond : float - The estimated reciprocal of the condition number of the leading row - (if leri = 'L') or the leading column (if leri = 'R') coefficient - matrix of P(s). - If rcond is nearly zero, P(s) is nearly row or column non-proper. - A : rank-2 array('d') with bounds (n,n) - The leading n-by-n part of this array contains the state dynamics matrix A. - B : rank-2 array('d') with bounds (n,max(m,p)) - The leading n-by-n part of this array contains the input/state matrix B; - the remainder of the leading n-by-max(m,p) part is used as internal - workspace. - C : rank-2 array('d') with bounds (max(m,p),n) - The leading p-by-n part of this array contains the state/output matrix C; - the remainder of the leading max(m,p)-by-n part is used as internal - workspace. - D : rank-2 array('d') with bounds (max(m,p),max(m,p)) - The leading p-by-m part of this array contains the direct transmission - matrix D; the remainder of the leading max(m,p)-by-max(m,p) part is - used as internal workspace. + Parameters + ---------- + m : int + The number of system inputs. m > 0. + p : int + The number of system outputs. p > 0. + lend(index) + index : (p) or (m) array_like + If leri = 'L', index(i), i = 1,2,...,p, must contain the maximum + degree of the polynomials in the I-th row of the denominator matrix + P(s) of the given left polynomial matrix representation. + If leri = 'R', index(i), i = 1,2,...,m, must contain the maximum + degree of the polynomials in the I-th column of the denominator + matrix P(s) of the given right polynomial matrix representation. + pcoeff : (p,p,*) or (m,m,*) array_like + If leri = 'L' then porm = p, otherwise porm = m. The leading + porm-by-porm-by-kpcoef part of this array must contain + the coefficients of the denominator matrix P(s). pcoeff(i,j,k) is + the coefficient in s**(index(iorj)-K+1) of polynomial (I,J) of P(s), + where k = 1,2,...,kpcoef and kpcoef = max(index) + 1; if leri = 'L' + then iorj = i, otherwise iorj = j. Thus for leri = 'L', + P(s) = diag(s**index)*(pcoeff(.,.,1)+pcoeff(.,.,2)/s+...). + If leri = 'R', pcoeff is modified by the routine but restored on exit. + qcoeff : (p, m, *) or (max(m,p), max(m,p), *) array_like + If leri = 'L' then porp = m, otherwise porp = p. The leading + porm-by-porp-by-kpcoef part of this array must contain + the coefficients of the numerator matrix Q(s). + qcoeff(i,j,k) is defined as for pcoeff(i,j,k). + If leri = 'R', qcoeff is modified by the routine but restored on exit. + leri : {'L', 'R'} + Indicates whether a left polynomial matrix representation or a right + polynomial matrix representation is input as follows: + = 'L': A left matrix fraction is input; + = 'R': A right matrix fraction is input. + ldwork : int, optional + The length of the cache array. ldwork >= max(m,p)*(max(m,p)+4) + For optimum performance it should be larger. + Default is None. + + Returns + ------- + n : int + The order of the resulting state-space representation. + That is, n = sum(index). + rcond : float + The estimated reciprocal of the condition number of the leading row + (if leri = 'L') or the leading column (if leri = 'R') coefficient + matrix of P(s). + If rcond is nearly zero, P(s) is nearly row or column non-proper. + A : (n, n) ndarray + The leading n-by-n part of this array contains the state dynamics matrix A. + B : rank-2 array('d') with bounds (n,max(m,p)) + The leading n-by-n part of this array contains the input/state matrix B; + the remainder of the leading n-by-max(m,p) part is used as internal + workspace. + C : (max(m,p), n) ndarray + The leading p-by-n part of this array contains the state/output matrix C; + the remainder of the leading max(m,p)-by-n part is used as internal + workspace. + D : (max(m,p), max(m,p)) ndarray + The leading p-by-m part of this array contains the direct transmission + matrix D; the remainder of the leading max(m,p)-by-max(m,p) part is + used as internal workspace. + Raises ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value SlycotArithmeticError :info == 1 and leri = 'L': P(s) is not row proper @@ -730,7 +832,6 @@ def tc04ad(m,p,index,pcoeff,qcoeff,leri,ldwork=None): raise_if_slycot_error(out[-1], arg_list, tc04ad.__doc__, locals()) return out[:-1] - def tc01od(m,p,indlin,pcoeff,qcoeff,leri): """ pcoeff,qcoeff = tc01od_l(m,p,indlim,pcoeff,qcoeff,leri) @@ -739,39 +840,46 @@ def tc01od(m,p,indlin,pcoeff,qcoeff,leri): polynomial matrix representations are of the form Q(s)*inv(P(s)) and inv(P(s))*Q(s) respectively. - Required arguments: - m : input int - The number of system inputs. m > 0. - p : input int - The number of system outputs. p > 0. - indlim : input int - The highest value of k for which pcoeff(.,.,k) and qcoeff(.,.,k) - are to be transposed. - k = kpcoef + 1, where kpcoef is the maximum degree of the polynomials - in P(s). indlim > 0. - pcoeff : input rank-3 array('d') with bounds (p,p,indlim) or (m,m,indlim) - If leri = 'L' then porm = p, otherwise porm = m. - On entry, the leading porm-by-porm-by-indlim part of this array - must contain the coefficients of the denominator matrix P(s). - pcoeff(i,j,k) is the coefficient in s**(indlim-k) of polynomial - (i,j) of P(s), where k = 1,2,...,indlim. - qcoeff : input rank-3 array('d') with bounds (max(m,p),max(m,p),indlim) - On entry, the leading p-by-m-by-indlim part of this array must - contain the coefficients of the numerator matrix Q(s). - qcoeff(i,j,k) is the coefficient in s**(indlim-k) of polynomial - (i,j) of Q(s), where k = 1,2,...,indlim. - leri : input string(len=1) - Return objects: - pcoeff : rank-3 array('d') with bounds (p,p,indlim) - On exit, the leading porm-by-porm-by-indlim part of this array - contains the coefficients of the denominator matrix P'(s) of - the dual system. - qcoeff : rank-3 array('d') with bounds (max(m,p),max(m,p),indlim) - On exit, the leading m-by-p-by-indlim part of the array contains - the coefficients of the numerator matrix Q'(s) of the dual system. - info : int - = 0: successful exit; - < 0: if info = -i, the i-th argument had an illegal value. + Parameters + ---------- + m : int + The number of system inputs. m > 0. + p : int + The number of system outputs. p > 0. + indlim : int + The highest value of k for which pcoeff(.,.,k) and qcoeff(.,.,k) + are to be transposed. + k = kpcoef + 1, where kpcoef is the maximum degree of the polynomials + in P(s). indlim > 0. + pcoeff : (p, p, indlim) or (m, m, indlim) array_like + If leri = 'L' then porm = p, otherwise porm = m. + On entry, the leading porm-by-porm-by-indlim part of this array + must contain the coefficients of the denominator matrix P(s). + pcoeff(i,j,k) is the coefficient in s**(indlim-k) of polynomial + (i,j) of P(s), where k = 1,2,...,indlim. + qcoeff : (max(m,p), max(m,p), indlim) array_like + On entry, the leading p-by-m-by-indlim part of this array must + contain the coefficients of the numerator matrix Q(s). + qcoeff(i,j,k) is the coefficient in s**(indlim-k) of polynomial + (i,j) of Q(s), where k = 1,2,...,indlim. + leri : {'L', 'R'} + = 'L': A left matrix fraction is input; + = 'R': A right matrix fraction is input. + + Returns + ------- + pcoeff : (p, p, indlim) ndarray + On exit, the leading porm-by-porm-by-indlim part of this array + contains the coefficients of the denominator matrix P'(s) of + the dual system. + qcoeff : (max(m,p), max(m,p), indlim) ndarray + On exit, the leading m-by-p-by-indlim part of the array contains + the coefficients of the numerator matrix Q'(s) of the dual system. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)' arg_list = ['leri', 'M', 'P', 'indlim', 'pcoeff', 'LDPCO1'+hidden, @@ -793,32 +901,40 @@ def tf01md(n,m,p,N,A,B,C,D,u,x0): To compute the output sequence of a linear time-invariant open-loop system given by its discrete-time state-space model - Required arguments: - n : input int - Order of the State-space representation. - m : input int - Number of inputs. - p : input int - Number of outputs. - N : input int - Number of output samples to be computed. - A : input rank-2 array('d') with bounds (n,n) - State dynamics matrix. - B : input rank-2 array('d') with bounds (n,m) - Input/state matrix. - C : input rank-2 array('d') with bounds (p,n) - State/output matrix. - D : input rank-2 array('d') with bounds (p,m) - Direct transmission matrix. - u : input rank-2 array('d') with bounds (m,N) - Input signal. - x0 : input rank-1 array('d') with bounds (n) - Initial state, at time 0. - Return objects: - xf : rank-1 array('d') with bounds (n) - Final state, at time N+1. - y : rank-2 array('d') with bounds (p,N) - Output signal. + Parameters + ---------- + n : int + Order of the State-space representation. + m : int + Number of inputs. + p : int + Number of outputs. + N : int + Number of output samples to be computed. + A : (n, n) array_like + State dynamics matrix. + B : (n, m) array_like + Input/state matrix. + C : (p, n) array_like + State/output matrix. + D : (p, m) array_like + Direct transmission matrix. + u : (m, n) + Input signal. + x0 : (n, ) array_like + Initial state, at time 0. + + Returns + ------- + xf : (n) ndarray + Final state, at time n+1. + y : (p, n) ndarray + Output signal. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)' arg_list = ['n','m','p','ny','A','lda'+hidden,'B','ldb'+hidden, @@ -839,28 +955,36 @@ def tf01rd(n,m,p,N,A,B,C,ldwork=None): All matrices are treated as dense, and hence TF01RD is not intended for large sparse problems. + Parameters + ---------- + n : int + Order of the State-space representation. + m : int + Number of inputs. + p : int + Number of outputs. + N : int + Number of Markov parameters to be computed. + A : (n, n) array_like + State dynamics matrix. + B : (n, m) array_like + Input/state matrix. + C : (p, n) array_like + State/output matrix. + ldwork : int, optional + The length of the array DWORK. + ldwork >= max(1, 2*n*p). - Required arguments: - n : input int - Order of the State-space representation. - m : input int - Number of inputs. - p : input int - Number of outputs. - N : input int - Number of Markov parameters to be computed. - A : input rank-2 array('d') with bounds (n,n) - State dynamics matrix. - B : input rank-2 array('d') with bounds (n,m) - Input/state matrix. - C : input rank-2 array('d') with bounds (p,n) - State/output matrix. - Optional arguments: - ldwork := 2*na*nc input int - Return objects: - H : rank-2 array('d') with bounds (p,N*m) - H[:,(k-1)*m : k*m] contains the k-th Markov parameter, - for k = 1,2...N. + Returns + ------- + H : (p, N*m) ndarray + H[:,(k-1)*m : k*m] contains the k-th Markov parameter, + for k = 1,2...N. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)' arg_list = ['n','m','p','N','A','lda'+hidden,'B','ldb'+hidden,'C', @@ -873,86 +997,6 @@ def tf01rd(n,m,p,N,A,B,C,ldwork=None): raise_if_slycot_error(out[-1], arg_list) return out[0] -def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None): - """Ar, Br, Cr, nr = tb01pd(n,m,p,A,B,C,[job,equil,tol,ldwork]) - - To find a reduced (controllable, observable, or minimal) state- - space representation (Ar,Br,Cr) for any original state-space - representation (A,B,C). The matrix Ar is in upper block - Hessenberg form. - - Required arguments: - n : input int - Order of the State-space representation. - m : input int - Number of inputs. - p : input int - Number of outputs. - A : input rank-2 array('d') with bounds (n,n) - State dynamics matrix. - B : input rank-2 array('d') with bounds (n,max(m,p)) - The leading n-by-m part of this array must contain the original - input/state matrix B; the remainder of the leading n-by-max(m,p) - part is used as internal workspace. - C : input rank-2 array('d') with bounds (p,n) - The leading p-by-n part of this array must contain the original - state/output matrix C; the remainder of the leading max(1,m,p)-by-n - part is used as internal workspace. - Optional arguments: - job : input char*1 - Indicates whether the user wishes to remove the - uncontrollable and/or unobservable parts as follows: - = 'M': Remove both the uncontrollable and unobservable - parts to get a minimal state-space representation; - = 'C': Remove the uncontrollable part only to get a - controllable state-space representation; - = 'O': Remove the unobservable part only to get an - observable state-space representation. - equil : input char*1 - Specifies whether the user wishes to preliminarily balance - the triplet (A,B,C) as follows: - = 'S': Perform balancing (scaling); - = 'N': Do not perform balancing. - Return objects: - Ar : output rank-2 array('d') with bounds (nr,nr) - Contains the upper block Hessenberg state dynamics matrix - Ar of a minimal, controllable, or observable realization - for the original system, depending on the value of JOB, - JOB = 'M', JOB = 'C', or JOB = 'O', respectively. - Br : output rank-2 array('d') with bounds (nr,m) - Contains the transformed input/state matrix Br of a - minimal, controllable, or observable realization for the - original system, depending on the value of JOB, JOB = 'M', - JOB = 'C', or JOB = 'O', respectively. If JOB = 'C', only - the first IWORK(1) rows of B are nonzero. - Cr : output rank-2 array('d') with bounds (p,nr) - - Contains the transformed state/output matrix Cr of a - minimal, C controllable, or observable realization for the - original C system, depending on the value of JOB, JOB = - 'M', C JOB = 'C', or JOB = 'O', respectively. C If JOB = - 'M', or JOB = 'O', only the last IWORK(1) columns C (in - the first NR columns) of C are nonzero. - nr : output int - The order of the reduced state-space representation - (Ar,Br,Cr) of a minimal, controllable, or observable - realization for the original system, depending on - JOB = 'M', JOB = 'C', or JOB = 'O'. - """ - hidden = ' (hidden by the wrapper)' - arg_list = ['job', 'equil', 'n','m','p','A','lda'+hidden,'B','ldb'+hidden, - 'C','ldc'+hidden,'nr','tol','iwork'+hidden,'dwork'+hidden, - 'ldwork','info'+hidden] - if ldwork is None: - ldwork = max(1, n+max(n,3*m,3*p)) - elif ldwork < max(1, n+max(n,3*m,3*p)): - raise SlycotParameterError("ldwork is too small", -15) - out = _wrapper.tb01pd(n=n,m=m,p=p,a=A,b=B,c=C, - job=job,equil=equil,tol=tol,ldwork=ldwork) - - raise_if_slycot_error(out[-1], arg_list) - return out[:-1] - 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]) @@ -981,64 +1025,74 @@ def tg01ad(l,n,m,p,A,E,B,C,thresh=0.0,job='A'): S = ( A-lambda E ). ( C ) - Required arguments: - l : input int - The number of rows of matrices A, B, and E. l >= 0. - n : input int - The number of columns of matrices A, E, and C. n >= 0. - m : input int - The number of columns of matrix B. m >= 0. - p : input int - The number of rows of matrix C. P >= 0. - A : rank-2 array('d') with bounds (l,n) - The leading L-by-N part of this array must - contain the state dynamics matrix A. - E : rank-2 array('d') with bounds (l,n) - The leading L-by-N part of this array must - contain the descriptor matrix E. - B : rank-2 array('d') with bounds (l,m) - The leading L-by-M part of this array must - contain the input/state matrix B. - The array B is not referenced if M = 0. - C : rank-2 array('d') with bounds (p,n) - The leading P-by-N part of this array must - contain the state/output matrix C. - The array C is not referenced if P = 0. - Optional arguments: - job := 'A' input string(len=1) - Indicates which matrices are involved in balancing, as - follows: - = 'A': All matrices are involved in balancing; - = 'B': B, A and E matrices are involved in balancing; - = 'C': C, A and E matrices are involved in balancing; - = 'N': B and C matrices are not involved in balancing. - thresh := 0.0 input float - Threshold value for magnitude of elements: - elements with magnitude less than or equal to - THRESH are ignored for balancing. THRESH >= 0. - Return objects: - A : rank-2 array('d') with bounds (l,n) - The leading L-by-N part of this array contains - the balanced matrix Dl*A*Dr. - E : rank-2 array('d') with bounds (l,n) - The leading L-by-N part of this array contains - the balanced matrix Dl*E*Dr. - B : rank-2 array('d') with bounds (l,m) - If M > 0, the leading L-by-M part of this array - contains the balanced matrix Dl*B. - The array B is not referenced if M = 0. - C : rank-2 array('d') with bounds (p,n) - If P > 0, the leading P-by-N part of this array - contains the balanced matrix C*Dr. - The array C is not referenced if P = 0. - lscale : rank-1 array('d') with bounds (l) - The scaling factors applied to S from left. If Dl(j) is - the scaling factor applied to row j, then - SCALE(j) = Dl(j), for j = 1,...,L. - rscale : rank-1 array('d') with bounds (n) - The scaling factors applied to S from right. If Dr(j) is - the scaling factor applied to column j, then - SCALE(j) = Dr(j), for j = 1,...,N. + + Parameters + ---------- + l : int + The number of rows of matrices A, B, and E. l >= 0. + n : int + The number of columns of matrices A, E, and C. n >= 0. + m : int + The number of columns of matrix B. m >= 0. + p : int + The number of rows of matrix C. P >= 0. + A : (l, n) array_like + The leading L-by-N part of this array must + contain the state dynamics matrix A. + E : (l, n) array_like + The leading L-by-N part of this array must + contain the descriptor matrix E. + B : (l, m) array_like + The leading L-by-M part of this array must + contain the input/state matrix B. + The array B is not referenced if M = 0. + C : (p, n) array_like + The leading P-by-N part of this array must + contain the state/output matrix C. + The array C is not referenced if P = 0. + job : {'A', 'B', 'C', 'N'}, optional + Indicates which matrices are involved in balancing, as + follows: + = 'A': All matrices are involved in balancing; + = 'B': B, A and E matrices are involved in balancing; + = 'C': C, A and E matrices are involved in balancing; + = 'N': B and C matrices are not involved in balancing. + Default is 'A'. + thresh : float, optional + Threshold value for magnitude of elements: + elements with magnitude less than or equal to + THRESH are ignored for balancing. THRESH >= 0. + Default is `0.0`. + + Returns + ------- + A : (l, n) ndarray + The leading L-by-N part of this array contains + the balanced matrix Dl*A*Dr. + E : (l, n) ndarray + The leading L-by-N part of this array contains + the balanced matrix Dl*E*Dr. + B : (l, m) ndarray + If M > 0, the leading L-by-M part of this array + contains the balanced matrix Dl*B. + The array B is not referenced if M = 0. + C : (p, n) ndarray + If P > 0, the leading P-by-N part of this array + contains the balanced matrix C*Dr. + The array C is not referenced if P = 0. + lscale : (l, ) ndarray + The scaling factors applied to S from left. If Dl(j) is + the scaling factor applied to row j, then + SCALE(j) = Dl(j), for j = 1,...,L. + rscale : (n, ) ndarray + The scaling factors applied to S from right. If Dr(j) is + the scaling factor applied to column j, then + SCALE(j) = Dr(j), for j = 1,...,N. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)' @@ -1072,132 +1126,146 @@ 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 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. - n : input int - The number of columns of matrices A, E, and C. n >= 0. - m : input int - The number of columns of matrix B. m >= 0. - p : input int - The number of rows of matrix C. p >= 0. - A : rank-2 array('d') with bounds (l,n) - The leading l-by-n part of this array must - contain the state dynamics matrix A. - E : rank-2 array('d') with bounds (l,n) - The leading l-by-n part of this array must - contain the descriptor matrix E. - B : rank-2 array('d') with bounds (l,m) - The leading L-by-M part of this array must - contain the input/state matrix B. - C : rank-2 array('d') with bounds (p,n) - The leading P-by-N part of this array must - contain the state/output matrix C. - Optional arguments: - Q : rank-2 array('d') with bounds (l,l) - If COMPQ = 'N': Q is not referenced. - If COMPQ = 'I': Q need not be set. - If COMPQ = 'U': The leading l-by-l part of this - array must contain an orthogonal matrix - Q1. - Z : rank-2 array('d') with bounds (n,n) - If COMPZ = 'N': Z is not referenced. - If COMPZ = 'I': Z need not be set. - If COMPZ = 'U': The leading n-by-n part of this - array must contain an orthogonal matrix - Z1. - compq := 'N' input string(len=1) - = 'N': do not compute Q. - = 'I': Q is initialized to the unit matrix, and the - orthogonal matrix Q is returned. - = 'U': Q must contain an orthogonal matrix Q1 on entry, - and the product Q1*Q is returned. - compz := 'N' input string(len=1) - = 'N': do not compute Z. - = 'I': Z is initialized to the unit matrix, and the - orthogonal matrix Z is returned. - = 'U': Z must contain an orthogonal matrix Z1 on entry, - and the product Z1*Z is returned. - joba := 'N' input string(len=1) - = 'N': do not reduce A22. - = 'R': reduce A22 to a SVD-like upper triangular form. - = 'T': reduce A22 to an upper trapezoidal form. - tol := 0 input float - The tolerance to be used in determining the rank of E - and of A22. If the user sets TOL > 0, then the given - value of TOL is used as a lower bound for the - reciprocal condition numbers of leading submatrices - of R or R22 in the QR decompositions E * P = Q * R of E - or A22 * P22 = Q22 * R22 of A22. - A submatrix whose estimated condition number is less than - 1/TOL is considered to be of full rank. If the user sets - TOL <= 0, then an implicitly computed, default tolerance, - defined by TOLDEF = L*N*EPS, is used instead, where - EPS is the machine precision (see LAPACK Library routine - DLAMCH). TOL < 1. - ldwork : input int - The length of the cache array. - ldwork >= MAX( 1, n+p, MIN(l,n)+MAX(3*n-1,m,l) ). - For optimal performance, ldwork should be larger. - Return objects: - A : rank-2 array('d') with bounds (l,n) - On entry, the leading L-by-N part of this array must - contain the state dynamics matrix A. - 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) - The leading L-by-M part of this array contains - the transformed matrix Q'*B. - C : rank-2 array('d') with bounds (p,n) - The leading P-by-N part of this array contains - the transformed matrix C*Z. - Q : rank-2 array('d') with bounds (l,l) - If COMPQ = 'N': Q is not referenced. - If COMPQ = 'I': The leading L-by-L part of this - array contains the orthogonal matrix Q, - where Q' is the product of Householder - transformations which are applied to A, - E, and B on the left. - If COMPQ = 'U': The leading L-by-L part of this - array contains the orthogonal matrix - Q1*Q. - Z : rank-2 array('d') with bounds (n,n) - If COMPZ = 'N': Z is not referenced. - If COMPZ = 'I': The leading N-by-N part of this - array contains the orthogonal matrix Z, - which is the product of Householder - transformations applied to A, E, and C - on the right. - If COMPZ = 'U': The leading N-by-N part of this - array contains the orthogonal matrix - Z1*Z. - ranke : output int - The estimated rank of matrix E, and thus also the order - of the invertible upper triangular submatrix Er. - rnka22 : output int - If JOBA = 'R' or 'T', then RNKA22 is the estimated rank of - matrix A22, and thus also the order of the invertible - upper triangular submatrix Ar. - If JOBA = 'N', then RNKA22 is not referenced. + Parameters + ---------- + l : int + The number of rows of matrices A, B, and E. l >= 0. + n : int + The number of columns of matrices A, E, and C. n >= 0. + m : int + The number of columns of matrix B. m >= 0. + p : int + The number of rows of matrix C. p >= 0. + A : (l, n) array_like + The leading l-by-n part of this array must + contain the state dynamics matrix A. + E : (l, n) array_like + The leading l-by-n part of this array must + contain the descriptor matrix E. + B : (l, m) array_like + The leading L-by-M part of this array must + contain the input/state matrix B. + C : (p, n) array_like + The leading P-by-N part of this array must + contain the state/output matrix C. + Q : (l, l) array_like, optional + If COMPQ = 'N': Q is not referenced. + If COMPQ = 'I': Q need not be set. + If COMPQ = 'U': The leading l-by-l part of this + array must contain an orthogonal matrix + Q1. + Default is None. + Z : (n, n) array_like, optional + If COMPZ = 'N': Z is not referenced. + If COMPZ = 'I': Z need not be set. + If COMPZ = 'U': The leading n-by-n part of this + array must contain an orthogonal matrix + Z1. + Default is None. + compq : {'N', 'I', 'U'}, optional + = 'N': do not compute Q. + = 'I': Q is initialized to the unit matrix, and the + orthogonal matrix Q is returned. + = 'U': Q must contain an orthogonal matrix Q1 on entry, + and the product Q1*Q is returned. + Default is 'N'. + compz : {'N', 'I', 'U'}, optional + = 'N': do not compute Z. + = 'I': Z is initialized to the unit matrix, and the + orthogonal matrix Z is returned. + = 'U': Z must contain an orthogonal matrix Z1 on entry, + and the product Z1*Z is returned. + Default is 'N'. + joba : {'N', 'R', 'T'}, optional + = 'N': do not reduce A22. + = 'R': reduce A22 to a SVD-like upper triangular form. + = 'T': reduce A22 to an upper trapezoidal form. + Default is 'N'. + tol : float, optional + The tolerance to be used in determining the rank of E + and of A22. If the user sets TOL > 0, then the given + value of TOL is used as a lower bound for the + reciprocal condition numbers of leading submatrices + of R or R22 in the QR decompositions E * P = Q * R of E + or A22 * P22 = Q22 * R22 of A22. + A submatrix whose estimated condition number is less than + 1/TOL is considered to be of full rank. If the user sets + TOL <= 0, then an implicitly computed, default tolerance, + defined by TOLDEF = L*N*EPS, is used instead, where + EPS is the machine precision (see LAPACK Library routine + DLAMCH). TOL < 1. + Default is `0.0`. + ldwork : int, optional + The length of the cache array. + ldwork >= MAX( 1, n+p, MIN(l,n)+MAX(3*n-1,m,l) ). + For optimal performance, ldwork should be larger. + Default is None. + + Returns + ------- + A : (l, n) ndarray + On entry, the leading L-by-N part of this array must + contain the state dynamics matrix A. + 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 : (l, n) ndarray + 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 : (l, m) ndarray + The leading L-by-M part of this array contains + the transformed matrix Q'*B. + C : (p, n) ndarray + The leading P-by-N part of this array contains + the transformed matrix C*Z. + Q : (l, l) ndarray + If COMPQ = 'N': Q is not referenced. + If COMPQ = 'I': The leading L-by-L part of this + array contains the orthogonal matrix Q, + where Q' is the product of Householder + transformations which are applied to A, + E, and B on the left. + If COMPQ = 'U': The leading L-by-L part of this + array contains the orthogonal matrix + Q1*Q. + Z : (n, n) ndarray + If COMPZ = 'N': Z is not referenced. + If COMPZ = 'I': The leading N-by-N part of this + array contains the orthogonal matrix Z, + which is the product of Householder + transformations applied to A, E, and C + on the right. + If COMPZ = 'U': The leading N-by-N part of this + array contains the orthogonal matrix + Z1*Z. + ranke : int + The estimated rank of matrix E, and thus also the order + of the invertible upper triangular submatrix Er. + rnka22 : int + If JOBA = 'R' or 'T', then RNKA22 is the estimated rank of + matrix A22, and thus also the order of the invertible + upper triangular submatrix Ar. + If JOBA = 'N', then RNKA22 is not referenced. + + Raises + ------ + SlycotParameterError + :info = -i: the i-th argument had an illegal value """ hidden = ' (hidden by the wrapper)'