Infinite loop when computing eigen values

[Makogan]

There seems to be a bug with the eigen value finder.

The following code produces a tri-diagonal matrix with real eigenvalues:

pub fn legendre_nodes(size: usize) -> Vec<f64>
{
    let n = size;
    let mut matrix = nalgebra::DMatrix::<f64>::zeros(n, n);

    let a_fn = |i| (2. * i as f64 - 1.) / i as f64;
    let c_fn = |i| (i as f64 - 1.) / i as f64;
    for i in 1..=n
    {
        let a = a_fn(i);
        let an = a_fn(i + 1);
        let cn = c_fn(i + 1);

        let beta = (cn / (a * an)).sqrt();

        if i < size
        {
            matrix[(i, i - 1)] = beta;
            matrix[(i - 1, i)] = beta;
        }
    }

    println!("{:.2}", matrix);

    let evals = matrix.eigenvalues().unwrap();
    let mut res: Vec<_> = evals.iter().copied().collect();
    for e in &evals
    {
        println!("{}", e);
    }

    res
}

Or alternatively:

  ┌                     ┐
  │ 0.00 0.58 0.00 0.00 │
  │ 0.58 0.00 0.52 0.00 │
  │ 0.00 0.52 0.00 0.51 │
  │ 0.00 0.00 0.51 0.00 │
  └                     ┘

I traced the problem down to delimit_subproblem . The start and the end get stuck at values 0 and 3 respectively and never update beyond that.

[alexhroom]

Edit: just to note this issue has been reported in #764, and #611, the former of which also finds it's an issue with delimit_subproblem; the issue below is a different one #1291 which has been fixed

I'm finding an infinite hang on a 3x3 zeroes matrix:

use nalgebra::Matrix3;

fn main() {
    let mat = Matrix3::<f64>::zeros();
    let _ = mat.eigenvalues(); 
}

also freezes for Matrix4, Matrix5, Matrix6 (but not for Matrix2) - although I don't know how you traced down the problem so I don't know for sure whether this is the same issue

[alexhroom]

https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=d82ff99fa31c5bdffabb81ce428abe99 for posterity here's an interactive version of your Rust code - you can set the size of the matrix with the n variable. interestingly I found it runs just fine if you add 0.00001 to each non-zero value.

[alexhroom]

with the example matrix in the issue, after 30 iterations the matrix t is (hard to look at):

  ┌                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ┐
  │                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                                      1.4915318439233651                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                      0.00000000000000007979735930427193 │
  │                                                                                                                                                                                                                                                                                                                      1.4915318439233636                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                     -0.00000000000000010408359384706564                                                                                                                                                                                                                                                                                                                                       0 │
  │                                                                                                                                                                                                                                                                                                                                       0 -0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000005                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                                      0.5888644410992601 │
  │                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                                                       0                                                                                                                                                                                                                                                                                                                      0.5888644410992601                                                                                                                                                                                                                                                                                                                                       0 │
  └                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ┘

and remains stuck there for infinite iterations afterwards.

If I arbitrarily set that very small float to zero (at iteration 30), the algorithm returns the eigenvalues +/- 0.8611..., -0.33998... and 0. The fourth eigenvalue should be +0.33998...