LA: optional argument on EVD decomposition if matrix is known to be symmetric #209

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -1303,9 +1303,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return DenseSvd.Create(this, computeVectors);
         }
 
-        public override Evd<Complex> Evd()
+        public override Evd<Complex> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return DenseEvd.Create(this);
+            return DenseEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs

@@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static DenseEvd Create(DenseMatrix matrix)
+        public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -78,14 +79,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             var blockDiagonal = new DenseMatrix(order);
             var eigenValues = new DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
-                }
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsConjugateSymmetric();
+                    break;
             }
 
             Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values);

src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs

@@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static UserEvd Create(Matrix<Complex> matrix)
+        public static UserEvd Create(Matrix<Complex> matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -78,14 +79,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             var blockDiagonal = Matrix<Complex>.Build.SameAs(matrix, order, order);
             var eigenValues = new DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
-                }
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsConjugateSymmetric();
+                    break;
             }
 
             if (isSymmetric)

src/Numerics/LinearAlgebra/Complex/Matrix.cs

@@ -710,9 +710,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return UserSvd.Create(this, computeVectors);
         }
 
-        public override Evd<Complex> Evd()
+        public override Evd<Complex> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return UserEvd.Create(this);
+            return UserEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -1300,9 +1300,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return DenseSvd.Create(this, computeVectors);
         }
 
-        public override Evd<Complex32> Evd()
+        public override Evd<Complex32> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return DenseEvd.Create(this);
+            return DenseEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs

@@ -63,9 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static DenseEvd Create(DenseMatrix matrix)
+        public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -79,14 +80,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             var blockDiagonal = new DenseMatrix(order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
-                }
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsConjugateSymmetric();
+                    break;
             }
 
             Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values);

src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs

@@ -61,9 +61,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static UserEvd Create(Matrix<Complex32> matrix)
+        public static UserEvd Create(Matrix<Complex32> matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -77,14 +78,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             var blockDiagonal = Matrix<Complex32>.Build.SameAs(matrix, order, order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
-                }
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsConjugateSymmetric();
+                    break;
             }
 
             if (isSymmetric)

src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@@ -704,9 +704,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return UserSvd.Create(this, computeVectors);
         }
 
-        public override Evd<Complex32> Evd()
+        public override Evd<Complex32> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return UserEvd.Create(this);
+            return UserEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -1230,9 +1230,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return DenseSvd.Create(this, computeVectors);
         }
 
-        public override Evd<double> Evd()
+        public override Evd<double> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return DenseEvd.Create(this);
+            return DenseEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs

@@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static DenseEvd Create(DenseMatrix matrix)
+        public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             var blockDiagonal = new DenseMatrix(order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
-                }
+                case Symmetricity.Symmetric:
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsSymmetric();
+                    break;
             }
 
             Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values);

src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs

@@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static UserEvd Create(Matrix<double> matrix)
+        public static UserEvd Create(Matrix<double> matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             var blockDiagonal = Matrix<double>.Build.SameAs(matrix, order, order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
-                }
+                case Symmetricity.Symmetric:
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsSymmetric();
+                    break;
             }
 
             var d = new double[order];

src/Numerics/LinearAlgebra/Double/Matrix.cs

@@ -689,9 +689,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return UserSvd.Create(this, computeVectors);
         }
 
-        public override Evd<double> Evd()
+        public override Evd<double> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return UserEvd.Create(this);
+            return UserEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Matrix.Solve.cs

@@ -78,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// Computes the EVD decomposition for a matrix.
         /// </summary>
         /// <returns>The EVD decomposition object.</returns>
-        public abstract Evd<T> Evd();
+        public abstract Evd<T> Evd(Symmetricity symmetricity = Symmetricity.Unknown);
 
 
 

src/Numerics/LinearAlgebra/Options.cs

@@ -58,4 +58,27 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         Include = 1
     }
+
+    public enum Symmetricity
+    {
+        /// <summary>
+        /// It is not known yet whether a matrix is symmetric or not.
+        /// </summary>
+        Unknown = 0,
+
+        /// <summary>
+        /// A matrix is symmetric
+        /// </summary>
+        Symmetric = 1,
+
+        /// <summary>
+        /// A matrix is complex conjugate symmetric.
+        /// </summary>
+        ConjugateSymmetric = 2,
+
+        /// <summary>
+        /// A matrix is not symmetric
+        /// </summary>
+        Asymmetric = 3
+    }
 }

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -1230,9 +1230,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return DenseSvd.Create(this, computeVectors);
         }
 
-        public override Evd<float> Evd()
+        public override Evd<float> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return DenseEvd.Create(this);
+            return DenseEvd.Create(this, symmetricity);
         }
     }
 }

src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs

@@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static DenseEvd Create(DenseMatrix matrix)
+        public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             var blockDiagonal = new DenseMatrix(order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
-                }
+                case Symmetricity.Symmetric:
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsSymmetric();
+                    break;
             }
 
             Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values);

src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs

@@ -61,9 +61,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
         /// </summary>
         /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="symmetricity">If it is known whether the matrix is symmetric or not the routine can skip checking it itself.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
         /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
-        public static UserEvd Create(Matrix<float> matrix)
+        public static UserEvd Create(Matrix<float> matrix, Symmetricity symmetricity)
         {
             if (matrix.RowCount != matrix.ColumnCount)
             {
@@ -77,14 +78,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             var blockDiagonal = Matrix<float>.Build.SameAs(matrix, order, order);
             var eigenValues = new LinearAlgebra.Complex.DenseVector(order);
 
-            var isSymmetric = true;
-
-            for (var i = 0; isSymmetric && i < order; i++)
+            bool isSymmetric;
+            switch (symmetricity)
             {
-                for (var j = 0; isSymmetric && j < order; j++)
-                {
-                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
-                }
+                case Symmetricity.Symmetric:
+                case Symmetricity.ConjugateSymmetric:
+                    isSymmetric = true;
+                    break;
+                case Symmetricity.Asymmetric:
+                    isSymmetric = false;
+                    break;
+                default:
+                    isSymmetric = matrix.IsSymmetric();
+                    break;
             }
 
             var d = new float[order];

src/Numerics/LinearAlgebra/Single/Matrix.cs

@@ -689,9 +689,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return UserSvd.Create(this, computeVectors);
         }
 
-        public override Evd<float> Evd()
+        public override Evd<float> Evd(Symmetricity symmetricity = Symmetricity.Unknown)
         {
-            return UserEvd.Create(this);
+            return UserEvd.Create(this, symmetricity);
         }
     }
 }