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);
         }
+
+        /// <summary>
+        /// See: https://github.com/mathnet/mathnet-numerics/issues/595
+        /// </summary>
+        [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<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);
+            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);
         }
+
+        /// <summary>
+        /// See: https://github.com/mathnet/mathnet-numerics/issues/595
+        /// </summary>
+        [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<Complex32>.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++)
                         {