From c71b4bd5a6bc64d64a83bc0d71bea1d5f3790c00 Mon Sep 17 00:00:00 2001 From: kkkmail Date: Sun, 28 Aug 2022 11:00:01 -0400 Subject: [PATCH] Fix EVD --- .../Complex/Factorization/EvdTests.cs | 27 +++++++++++++++++++ .../Complex32/Factorization/EvdTests.cs | 27 +++++++++++++++++++ .../ManagedLinearAlgebraProvider.Complex.cs | 19 ++++++++++--- .../ManagedLinearAlgebraProvider.Complex32.cs | 19 ++++++++++--- 4 files changed, 84 insertions(+), 8 deletions(-) diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs index 59a6a70c0..7ce7861d3 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs @@ -259,5 +259,32 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// See: https://github.com/mathnet/mathnet-numerics/issues/595 + /// + [Test] + public void CanFactorizeMatrixWithZeroInternalNorm() + { + 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.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); + AssertHelpers.AlmostEqualRelative(Av, Lv, 8); + } } } diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs index ff1cf8efd..54bf2b7e5 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs @@ -257,5 +257,32 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// See: https://github.com/mathnet/mathnet-numerics/issues/595 + /// + [Test] + public void CanFactorizeMatrixWithZeroInternalNorm() + { + Complex32[,] data = + { + {new Complex32(0.0f, 0.0f), new Complex32(1.0f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f)}, + {new Complex32(2.25f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f)}, + {new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(1.0f, 0.0f)}, + {new Complex32(0.0f, 0.0f), new Complex32(0.0f, 0.0f), new Complex32(2.25f, 0.0f), new Complex32(0.0f, 0.0f)} + }; + + var A = Matrix.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); + AssertHelpers.AlmostEqualRelative(Av, Lv, 8); + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 97ca54148..fb5e7db35 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -3064,10 +3064,21 @@ internal static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, C var im1Oim1 = im1O + im1; s = matrixH[im1O + i].Real; norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); - x = matrixH[im1Oim1]/norm; - vectorV[i - 1] = x; - matrixH[im1Oim1] = norm; - matrixH[im1O + i] = new Complex(0.0, s.Real/norm); + + 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); + } for (var j = i; j < order; j++) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index fbbc96f08..7ea62cf69 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -3066,10 +3066,21 @@ internal static void NonsymmetricReduceHessenberToRealSchur(Complex32[] vectorV, var im1Oim1 = im1O + im1; s = matrixH[im1O + i].Real; norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); - x = matrixH[im1Oim1]/norm; - vectorV[i - 1] = x; - matrixH[im1Oim1] = norm; - matrixH[im1O + i] = new Complex32(0.0f, s.Real/norm); + + if (norm != 0.0) + { + x = matrixH[im1Oim1] / norm; + vectorV[i - 1] = x; + matrixH[im1Oim1] = norm; + matrixH[im1O + i] = new Complex32(0.0f, s.Real / norm); + } + else + { + x = 1.0f; + vectorV[i - 1] = x; + matrixH[im1Oim1] = norm; + matrixH[im1O + i] = new Complex32(0.0f, 0.0f); + } for (var j = i; j < order; j++) {