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

Changes from 1 commit
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
Prev Previous commit
Next Next commit
created new test for permutations, removed double runs
  • Loading branch information
@repagh
repagh committed May 13, 2020
commit b18c7a16bf7c772de7c5feee95a72290878b3488
98 changes: 57 additions & 41 deletions slycot/tests/test_tb05ad.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,31 +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 Down Expand Up @@ -58,13 +55,6 @@ def test_tb05ad_ag(self):
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_nh(self):
"""Test that tb05ad with job = 'NH' computes the correct
frequency response after conversion to Hessenberg form.
Expand Down Expand Up @@ -180,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 @@ -198,27 +188,53 @@ def test_tb05ad_resonance(self):
assert cm.exception.info == 2

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

Based on https://rdrr.io/rforge/expm/man/balance.html
Tests for the cause of the problem reported in issue #11
balancing permutations were not correctly applied to the
C and D matrix.
"""
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)))

# 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
crit = False
while not crit:
A = np.random.randn(8, 8)
A[np.random.uniform(size=(8, 8)) > 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))

crit = nperm < 8 and nperm >= 4 and \
neig0 > 1 and neig0 <= 3

# 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(
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')
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))
assert_almost_equal(eig(At)[0], eig(A)[0])

@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.

Is the order and scaling of eigenvalues guaranteed? Could run into assertion failures on some build systems. Or do scipy.linalg.eig() and MB05AD call the same LAPACK code and thus always yield the same results?

@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.

I'll put in some sorting.

@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.

What about scaling, including different signs? I had to deal with that in test_mb05md() (#79)

@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.

Now you're confusing me. Do you mean scaling of the eigenvalues? But then the frequency should also be scaled to calculate the response, and I don't see anything on that in the TB05AD.f code. Eigenvectors may be differently scaled, but I did not do a test on that.

Or do you mean the scaling in the matrix_balance / dgebal step. I now added some code for that, got to the same testing case as before, on the 10th matrix tried.

@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.

Okay strike that. Brainfart by me. Mixed up eigenvalue properties with eigenvectors.



if __name__ == "__main__":
Expand Down