From 640b008fa5b68a8882691529248378e8e57d321c Mon Sep 17 00:00:00 2001 From: daingun Date: Fri, 1 Nov 2019 23:27:08 +0100 Subject: [PATCH] Use same algorithm to solve 2x2 eigenvalue problem The eigenvalue problem is solved in two different method that use different methods to calculate the discriminant of the solution to the quadratic equation. Use the method whose computation is considered more stable. --- src/linalg/schur.rs | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index 2a2bb250..5cd90a94 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -309,16 +309,17 @@ where let hmn = t[(m, n)]; let hnn = t[(n, n)]; - let tra = hnn + hmm; - let det = hnn * hmm - hnm * hmn; - let discr = tra * tra * crate::convert(0.25) - det; + // NOTE: use the same algorithm as in compute_2x2_eigvals. + let val = (hmm - hnn) * crate::convert(0.5); + let discr = hnm * hmn + val * val; // All 2x2 blocks have negative discriminant because we already decoupled those - // with positive eigenvalues.. + // with positive eigenvalues. let sqrt_discr = NumComplex::new(N::zero(), (-discr).sqrt()); - out[m] = NumComplex::new(tra * crate::convert(0.5), N::zero()) + sqrt_discr; - out[m + 1] = NumComplex::new(tra * crate::convert(0.5), N::zero()) - sqrt_discr; + let half_tra = (hnn + hmm) * crate::convert(0.5); + out[m] = NumComplex::new(half_tra, N::zero()) + sqrt_discr; + out[m + 1] = NumComplex::new(half_tra, N::zero()) - sqrt_discr; m += 2; } @@ -413,7 +414,6 @@ where let inv_rot = rot.inverse(); inv_rot.rotate(&mut m); rot.rotate_rows(&mut m); - m[(1, 0)] = N::zero(); if compute_q { // XXX: we have to build the matrix manually because