M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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.
This commit is contained in:
daingun 2019-11-01 23:27:08 +01:00 committed by GitHub
parent ead2360f8e
commit 640b008fa5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 7 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/linalg/schur.rs
View File

@ -309,16 +309,17 @@ where
let hmn = t[(m, n)]; let hmn = t[(m, n)];
let hnn = t[(n, n)]; let hnn = t[(n, n)];
let tra = hnn + hmm; // NOTE: use the same algorithm as in compute_2x2_eigvals.
let det = hnn * hmm - hnm * hmn; let val = (hmm - hnn) * crate::convert(0.5);
let discr = tra * tra * crate::convert(0.25) - det; let discr = hnm * hmn + val * val;
// All 2x2 blocks have negative discriminant because we already decoupled those // 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()); let sqrt_discr = NumComplex::new(N::zero(), (-discr).sqrt());
out[m] = NumComplex::new(tra * crate::convert(0.5), N::zero()) + sqrt_discr; let half_tra = (hnn + hmm) * crate::convert(0.5);
out[m + 1] = NumComplex::new(tra * crate::convert(0.5), N::zero()) - sqrt_discr; out[m] = NumComplex::new(half_tra, N::zero()) + sqrt_discr;
out[m + 1] = NumComplex::new(half_tra, N::zero()) - sqrt_discr;
m += 2; m += 2;
} }
@ -413,7 +414,6 @@ where
let inv_rot = rot.inverse(); let inv_rot = rot.inverse();
inv_rot.rotate(&mut m); inv_rot.rotate(&mut m);
rot.rotate_rows(&mut m); rot.rotate_rows(&mut m);
m[(1, 0)] = N::zero();
if compute_q { if compute_q {
// XXX: we have to build the matrix manually because // XXX: we have to build the matrix manually because