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

Merge pull request #669 from daingun/patch-2 

Use same algorithm to solve 2x2 eigenvalue problem
This commit is contained in:
Sébastien Crozet 2019-11-05 23:02:13 +01:00
parent 2f77d73226 a5b2a4ceb7
commit d78309b1fd
1 changed files with 7 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
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;
} }