Matrix.Evd() hangs up

[kkkmail]

Given the following complex matrix:

   (0., 0.)  (1., 0.)     (0., 0.)  (0., 0.)
(2.25, 0.)  (0., 0.)     (0., 0.)  (0., 0.)
   (0., 0.)  (0., 0.)     (0., 0.)  (1., 0.)
   (0., 0.)  (0., 0.)  (2.25, 0.)  (0., 0.)

call to Evd() hangs up and never comes back. See enclosed FSX script.

#r "System.Numerics.dll"
#r "../packages/MathNet.Numerics.4.5.1/lib/net461/MathNet.Numerics.dll"
#r "../packages/MathNet.Numerics.FSharp.4.5.1/lib/net45/MathNet.Numerics.FSharp.dll"

open System.Numerics
open MathNet.Numerics.LinearAlgebra

let cplx a = Complex(a, 0.0)

let c = 
    [
        [cplx 0.0; cplx 1.0; cplx 0.0; cplx 0.0]
        [cplx 2.25; cplx 0.0; cplx 0.0; cplx 0.0]
        [cplx 0.0; cplx 0.0; cplx 0.0; cplx 1.0]
        [cplx 0.0; cplx 0.0; cplx 2.25; cplx 0.0]
    ]
    |> matrix

printfn "c = %A" c

#time
printfn "Calculating Evd"
let ec = c.Evd()
printfn "ec = %A" ec
#time

EvdErr_01.zip

[kkkmail]

The bug was after line 668 here: https://github.com/mathnet/mathnet-numerics/blob/master/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs#L668

 norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real);

If norm comes out as exact zero (which is the case for the matrix above), then the next lines blow up:

x = matrixH[im1Oim1]/norm;

and

matrixH[im1O + i] = new Complex(0.0, s.Real/norm);

A solution turned out to be to replace that block:

x = matrixH[im1Oim1]/norm;
vectorV[i - 1] = x;
matrixH[im1Oim1] = norm;
matrixH[im1O + i] = new Complex(0.0, s.Real/norm);

with this:

if (norm != 0.0)
{
    x = matrixH[im1Oim1] / norm;
    vectorV[i - 1] = x;
    matrixH[im1Oim1] = norm;
    matrixH[im1O + i] = new Complex(0.0, s.Real / norm);
}
else
{
    x = 1.0;
    vectorV[i - 1] = x;
    matrixH[im1Oim1] = norm;
    matrixH[im1O + i] = new Complex(0.0, 0.0);
}

The following test shows that Evd() still works correctly (at least for this matrix).

[Test]
public void CanFactorizeDoesNotHangWhenComplex()
{
    Complex[,] data =
    {
        {new Complex(0.0, 0.0), new Complex(1.0, 0.0), new Complex(0.0, 0.0), new Complex(0.0, 0.0)},
        {new Complex(2.25, 0.0), new Complex(0.0, 0.0), new Complex(0.0, 0.0), new Complex(0.0, 0.0)},
        {new Complex(0.0, 0.0), new Complex(0.0, 0.0), new Complex(0.0, 0.0), new Complex(1.0, 0.0)},
        {new Complex(0.0, 0.0), new Complex(0.0, 0.0), new Complex(2.25, 0.0), new Complex(0.0, 0.0)}
    };

    var A = Matrix<Complex>.Build.DenseOfArray(data);

    var factorEvd = A.Evd();
    var V = factorEvd.EigenVectors;
    var λ = factorEvd.D;

    // Verify A*V = λ*V
    var Av = A * V;
    var Lv = V * λ;
    AssertHelpers.AlmostEqual(Av, Lv, 4);
}

[kkkmail]

It's been over 8 months. I can submit a PR with the fix if you give me access right to do so. Thanks.

[kkkmail]

@cdrnet It's been over two years and the bug is now in at least 4 places as far as I see from the codebase. I'd be glad to submit the PR, provided that you give me the access rights to do so. Thanks.

[gismofx]

@kkkmail I'm experiencing this hanging as well. Do you have a fork/dll/nuget that I can try?

[kkkmail]

@gismofx Here is the PR: #944 Sorry that it too so long. I gave up on any response and did not look here for a good while.