Skip to content

Dgebal application order #122

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
repagh merged 5 commits into from
May 16, 2020
Merged

Dgebal application order #122

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
fca28a3
fix for DGEBAL transformation application, issue #11
repagh May 11, 2020
41c162c
tweak tb05ad test cases
repagh May 12, 2020
b18c7a1
created new test for permutations, removed double runs
repagh May 13, 2020
b67fbf9
order eigenvalues before check
repagh May 15, 2020
6f1d2e4
ensure dgebal does some scaling too on the test case
repagh May 15, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion slycot/src/TB01TD.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ SUBROUTINE TB01TD( N, M, P, A, LDA, B, LDB, C, LDC, D, LDD, LOW,
C to transform B and C.
C
DO 10 K = 1, N
KOLD = K
KOLD = N + 1 - K ! RvP, rabraker, slycot #11
IF ( ( LOW.GT.KOLD ) .OR. ( KOLD.GT.IGH ) ) THEN
IF ( KOLD.LT.LOW ) KOLD = LOW - KOLD
KNEW = INT( SCSTAT(KOLD) )
Expand Down
4 changes: 2 additions & 2 deletions slycot/src/TB05AD.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,10 +346,10 @@ SUBROUTINE TB05AD( BALEIG, INITA, N, M, P, FREQ, A, LDA, B, LDB,
C defined in the subroutine DGEBAL.
C
DO 10 J = 1, N
JJ = J
JJ = N + 1 - J ! RvP, rabraker, slycot #11
IF ( JJ.LT.LOW .OR. JJ.GT.IGH ) THEN
IF ( JJ.LT.LOW ) JJ = LOW - JJ
JP = DWORK(JJ)
JP = INT( DWORK(JJ) )
IF ( JP.NE.JJ ) THEN
C
C Permute rows of B.
Expand Down
41 changes: 37 additions & 4 deletions slycot/tests/test_tb05ad.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
from numpy.testing import assert_almost_equal



@bnavigator bnavigator May 12, 2020

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

E303, too many blank lines

Suggested change

@repagh repagh May 13, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

Should be fixed in the commit.

# set the random seed so we can get consistent results.
np.random.seed(40)
CASES = {}
Expand Down Expand Up @@ -48,10 +49,19 @@ def test_tb05ad_ng(self):
sys = CASES[key]
self.check_tb05ad_AG_NG(sys, 10*1j, 'NG')

@unittest.expectedFailure
def test_tb05ad_ag_failure(self):
""" Test tb05ad and job 'AG' (i.e., balancing enabled) fails
on certain A matrices.
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_ag_fixed_bug_no11(self):
""" Test tb05ad and job 'AG' (i.e., balancing enabled).

Failed on certain A matrices before, issue #11.
"""
self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG')

@bnavigator bnavigator May 12, 2020

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This just runs CASES['fail1'] twice. Probably the name 'fail1' is not appropriate anymore.

Suggested change
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_ag_fixed_bug_no11(self):
""" Test tb05ad and job 'AG' (i.e., balancing enabled).
Failed on certain A matrices before, issue #11.
"""
self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG')
def test_tb05ad_ag(self):
"""
Test that tb05ad with job 'AG' computes the correct
frequency response with balancing A first.
Failed on certain A matrices before, issue #11.
"""
for key in CASES:
sys = CASES[key]
self.check_tb05ad_AG_NG(sys, 10*1j, 'AG')

@repagh repagh May 13, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

renamed, and removed the repeat


Expand Down Expand Up @@ -187,6 +197,29 @@ def test_tb05ad_resonance(self):
transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH')
assert cm.exception.info == 2

def test_tb05ad_balance(self):
"""Specifically check for the balancing output.

Based on https://rdrr.io/rforge/expm/man/balance.html

@bnavigator bnavigator May 12, 2020

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

According to that link, this case has SCALE=[1 1 1 3], ILO=2, IHI=3 and thus only swaps 4 with 3 and 1 with 1. Does not really check the correct order of transformations done, does it?

To check the newly corrected transformation in TB05AD, Bt, Ct need to be checked for correct order as well.

@repagh repagh May 13, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

I created a new test, based on a random matrix, tweaked until that one showed permutation

@bnavigator bnavigator May 13, 2020

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Could this be pinned to a fixed found matrix? Having random data in tests that could fail at times and does not at other times will make debugging hard.

@repagh repagh May 15, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

we're still covered by the random seed at the start of the file

@bnavigator bnavigator May 15, 2020

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Ah okay, I didn't see that.

"""
A = np.array(((-1, -1, 0, 0),
( 0, 0, 10, 10),
( 0, 0, 10, 0),
( 0, 10, 0, 0)), dtype=float)
B = np.eye(4)
C = np.eye(4)
Ar = np.array(((-1, -1, 0, 0),
( 0, 0, 10, 10),
( 0, 10, 0, 0),
( 0, 0, 0, 10)))
jomega = 1.0
At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad(
4, 4, 4, jomega, A, B, C, job='AG')
assert_almost_equal(At, Ar)

At, Bt, Ct, rcond, g_jwb, ev, hinvb, info = transform.tb05ad(
4, 4, 4, jomega, A, B, C, job='AG')


if __name__ == "__main__":
unittest.main()