Open kkkmail opened 6 years ago
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);
}
It's been over 8 months. I can submit a PR with the fix if you give me access right to do so. Thanks.
@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.
@kkkmail I'm experiencing this hanging as well. Do you have a fork/dll/nuget that I can try?
@gismofx Here is the PR: https://github.com/mathnet/mathnet-numerics/pull/944 Sorry that it too so long. I gave up on any response and did not look here for a good while.
Given the following complex matrix:
call to Evd() hangs up and never comes back. See enclosed FSX script.
EvdErr_01.zip