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 all 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
97 changes: 76 additions & 21 deletions slycot/tests/test_tb05ad.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,30 +3,28 @@

from slycot import transform
from slycot.exceptions import SlycotArithmeticError, SlycotParameterError

import sys
import numpy as np

from scipy.linalg import matrix_balance, eig
import unittest
from numpy.testing import assert_almost_equal


# set the random seed so we can get consistent results.
np.random.seed(40)
CASES = {}

# This is a known failure for tb05ad when running job 'AG'
CASES['fail1'] = {'A': np.array([[-0.5, 0., 0., 0. ],
[ 0., -1., 0. , 0. ],
[ 1., 0., -0.5, 0. ],
[ 0., 1., 0., -1. ]]),
'B': np.array([[ 1., 0.],
[ 0., 1.],
[ 0., 0.],
[ 0., 0.]]),
'C': np.array([[ 0., 1., 1., 0.],
[ 0., 1., 0., 1.],
[ 0., 1., 1., 1.]])}
# This was (pre 2020) a known failure for tb05ad when running job 'AG'
CASES['known'] = {'A': np.array([[-0.5, 0., 0., 0.],
[ 0., -1., 0., 0.],
[ 1., 0., -0.5, 0.],
[ 0., 1., 0., -1.]]),
'B': np.array([[ 1., 0.],
[ 0., 1.],
[ 0., 0.],
[ 0., 0.]]),
'C': np.array([[ 0., 1., 1., 0.],
[ 0., 1., 0., 1.],
[ 0., 1., 1., 1.]])}

n = 20
p = 10
Expand All @@ -48,12 +46,14 @@ 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):
"""
self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG')
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
Expand Down Expand Up @@ -170,7 +170,7 @@ def check_tb05ad_errors(self, sys):
def test_tb05ad_resonance(self):
""" Test tb05ad resonance failure.

Actually test one of the exception messages. For many routines these
Actually test one of the exception messages. These
are parsed from the docstring, tests both the info index and the
message
"""
Expand All @@ -187,6 +187,61 @@ 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):
"""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])


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