From 667c49d0e18e9bd6d561dc71b73004daee3faae1 Mon Sep 17 00:00:00 2001 From: daingun Date: Fri, 1 Nov 2019 22:12:59 +0100 Subject: [PATCH 1/2] Correct Schur decomposition for 2x2 matrices Due to rounding and possible loss of precision the lower left element of the 2x2 matrix may be different from zero. --- src/linalg/schur.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index b31be9f6..2a2bb250 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -413,6 +413,7 @@ 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 From 50417494ecae77b2b47e1abd389e7d68d995c511 Mon Sep 17 00:00:00 2001 From: daingun Date: Fri, 1 Nov 2019 23:27:08 +0100 Subject: [PATCH 2/2] 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