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 Sébastien Crozet
parent 1103d49b80
commit 50417494ec
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 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