From c550921bcfe3f3a0fb56863f5e948b8c9234530b 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 640b008fa5b68a8882691529248378e8e57d321c 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