Skip to content

Commit

Permalink
FIX: Better detecting ill-conditioned pencils
Browse files Browse the repository at this point in the history
  • Loading branch information
ilayn committed Nov 15, 2017
1 parent 74e4586 commit fdb3889
Showing 1 changed file with 6 additions and 4 deletions.
10 changes: 6 additions & 4 deletions scipy/linalg/_solvers.py
Original file line number Diff line number Diff line change
Expand Up @@ -513,11 +513,12 @@ def solve_continuous_are(a, b, q, r, e=None, s=None, balanced=True):
if balanced:
x *= sca[:m, None] * sca[:m]

# Check the deviation from symmetry for success
# Check the deviation from symmetry for lack of success
# See proof of Thm.5 item 3 in [2]
u_sym = u00.conj().T.dot(u10)
n_u_sym = norm(u_sym, 1)
u_sym = u_sym - u_sym.conj().T
sym_threshold = np.max([np.spacing(1000.), n_u_sym])
sym_threshold = np.max([np.spacing(1000.), 0.1*n_u_sym])

if norm(u_sym, 1) > sym_threshold:
raise LinAlgError('The associated Hamiltonian pencil has eigenvalues '
Expand Down Expand Up @@ -720,11 +721,12 @@ def solve_discrete_are(a, b, q, r, e=None, s=None, balanced=True):
if balanced:
x *= sca[:m, None] * sca[:m]

# Check the deviation from symmetry for success
# Check the deviation from symmetry for lack of success
# See proof of Thm.5 item 3 in [2]
u_sym = u00.conj().T.dot(u10)
n_u_sym = norm(u_sym, 1)
u_sym = u_sym - u_sym.conj().T
sym_threshold = np.max([np.spacing(1000.), n_u_sym])
sym_threshold = np.max([np.spacing(1000.), 0.1*n_u_sym])

if norm(u_sym, 1) > sym_threshold:
raise LinAlgError('The associated symplectic pencil has eigenvalues'
Expand Down

0 comments on commit fdb3889

Please sign in to comment.