Skip to content

Commit b23ceb8

Browse files
repaghrepagh
authored
Merge pull request #122 from repagh/dgebal-application-order
Dgebal application order
2 parents 3431401 + 6f1d2e4 commit b23ceb8

3 files changed

Lines changed: 79 additions & 24 deletions

File tree

‎slycot/src/TB01TD.f‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -245,7 +245,7 @@ SUBROUTINE TB01TD( N, M, P, A, LDA, B, LDB, C, LDC, D, LDD, LOW,
245245
C to transform B and C.
246246
C
247247
DO 10 K = 1, N
248-
KOLD = K
248+
KOLD = N + 1 - K ! RvP, rabraker, slycot #11
249249
IF ( ( LOW.GT.KOLD ) .OR. ( KOLD.GT.IGH ) ) THEN
250250
IF ( KOLD.LT.LOW ) KOLD = LOW - KOLD
251251
KNEW = INT( SCSTAT(KOLD) )

‎slycot/src/TB05AD.f‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -346,10 +346,10 @@ SUBROUTINE TB05AD( BALEIG, INITA, N, M, P, FREQ, A, LDA, B, LDB,
346346
C defined in the subroutine DGEBAL.
347347
C
348348
DO 10 J = 1, N
349-
JJ = J
349+
JJ = N + 1 - J ! RvP, rabraker, slycot #11
350350
IF ( JJ.LT.LOW .OR. JJ.GT.IGH ) THEN
351351
IF ( JJ.LT.LOW ) JJ = LOW - JJ
352-
JP = DWORK(JJ)
352+
JP = INT( DWORK(JJ) )
353353
IF ( JP.NE.JJ ) THEN
354354
C
355355
C Permute rows of B.

‎slycot/tests/test_tb05ad.py‎

Lines changed: 76 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -3,30 +3,28 @@
33

44
from slycot import transform
55
from slycot.exceptions import SlycotArithmeticError, SlycotParameterError
6-
76
import sys
87
import numpy as np
9-
8+
from scipy.linalg import matrix_balance, eig
109
import unittest
1110
from numpy.testing import assert_almost_equal
1211

13-
1412
# set the random seed so we can get consistent results.
1513
np.random.seed(40)
1614
CASES = {}
1715

18-
# This is a known failure for tb05ad when running job 'AG'
19-
CASES['fail1'] = {'A': np.array([[-0.5, 0., 0., 0. ],
20-
[ 0., -1., 0. , 0. ],
21-
[ 1., 0., -0.5, 0. ],
22-
[ 0., 1., 0., -1. ]]),
23-
'B': np.array([[ 1., 0.],
24-
[ 0., 1.],
25-
[ 0., 0.],
26-
[ 0., 0.]]),
27-
'C': np.array([[ 0., 1., 1., 0.],
28-
[ 0., 1., 0., 1.],
29-
[ 0., 1., 1., 1.]])}
16+
# This was (pre 2020) a known failure for tb05ad when running job 'AG'
17+
CASES['known'] = {'A': np.array([[-0.5, 0., 0., 0.],
18+
[ 0., -1., 0., 0.],
19+
[ 1., 0., -0.5, 0.],
20+
[ 0., 1., 0., -1.]]),
21+
'B': np.array([[ 1., 0.],
22+
[ 0., 1.],
23+
[ 0., 0.],
24+
[ 0., 0.]]),
25+
'C': np.array([[ 0., 1., 1., 0.],
26+
[ 0., 1., 0., 1.],
27+
[ 0., 1., 1., 1.]])}
3028

3129
n = 20
3230
p = 10
@@ -48,12 +46,14 @@ def test_tb05ad_ng(self):
4846
sys = CASES[key]
4947
self.check_tb05ad_AG_NG(sys, 10*1j, 'NG')
5048

51-
@unittest.expectedFailure
52-
def test_tb05ad_ag_failure(self):
53-
""" Test tb05ad and job 'AG' (i.e., balancing enabled) fails
54-
on certain A matrices.
49+
def test_tb05ad_ag(self):
5550
"""
56-
self.check_tb05ad_AG_NG(CASES['fail1'], 10*1j, 'AG')
51+
Test that tb05ad with job 'AG' computes the correct
52+
frequency response.
53+
"""
54+
for key in CASES:
55+
sys = CASES[key]
56+
self.check_tb05ad_AG_NG(sys, 10*1j, 'AG')
5757

5858
def test_tb05ad_nh(self):
5959
"""Test that tb05ad with job = 'NH' computes the correct
@@ -170,7 +170,7 @@ def check_tb05ad_errors(self, sys):
170170
def test_tb05ad_resonance(self):
171171
""" Test tb05ad resonance failure.
172172
173-
Actually test one of the exception messages. For many routines these
173+
Actually test one of the exception messages. These
174174
are parsed from the docstring, tests both the info index and the
175175
message
176176
"""
@@ -187,6 +187,61 @@ def test_tb05ad_resonance(self):
187187
transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH')
188188
assert cm.exception.info == 2
189189

190+
def test_tb05ad_balance(self):
191+
"""Test balancing in tb05ad.
192+
193+
Tests for the cause of the problem reported in issue #11
194+
balancing permutations were not correctly applied to the
195+
C and D matrix.
196+
"""
197+
198+
# find a good test case. Some sparsity,
199+
# some zero eigenvalues, some non-zero eigenvalues,
200+
# and proof that the 1st step, with dgebal, does some
201+
# permutation and some scaling
202+
crit = False
203+
n = 8
204+
while not crit:
205+
A = np.random.randn(n, n)
206+
A[np.random.uniform(size=(n, n)) > 0.35] = 0.0
207+
208+
Aeig = eig(A)[0]
209+
neig0 = np.sum(np.abs(Aeig) == 0)
210+
As, T = matrix_balance(A)
211+
nperm = np.sum(np.diag(T == 0))
212+
nscale = n - np.sum(T == 1.0)
213+
crit = nperm < n and nperm >= n//2 and \
214+
neig0 > 1 and neig0 <= 3 and nscale > 0
215+
216+
# print("number of permutations", nperm, "eigenvalues=0", neig0)
217+
B = np.random.randn(8, 4)
218+
C = np.random.randn(3, 8)
219+
220+
# do a run
221+
jomega = 1.0
222+
At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad(
223+
8, 4, 3, jomega, A, B, C, job='AG')
224+
225+
# remove information on Q, in lower sub-triangle part of A
226+
At = np.triu(At, k=-1)
227+
228+
# now after the balancing in DGEBAL, and conversion to
229+
# upper Hessenberg form:
230+
# At = Q^T * (P^-1 * A * P ) * Q
231+
# with Q orthogonal
232+
# Ct = C * P * Q
233+
# Bt = Q^T * P^-1 * B
234+
# so test with Ct * At * Bt == C * A * B
235+
# and verify that eigenvalues of both A matrices are close
236+
assert_almost_equal(np.dot(np.dot(Ct, At), Bt),
237+
np.dot(np.dot(C, A), B))
238+
# uses a sort, there is no guarantee on the order of eigenvalues
239+
eigAt = eig(At)[0]
240+
idxAt = np.argsort(eigAt)
241+
eigA = eig(A)[0]
242+
idxA = np.argsort(eigA)
243+
assert_almost_equal(eigA[idxA], eigAt[idxAt])
244+
190245

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

0 commit comments

Comments
 (0)