diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index 49ef76c5..cb6b4391 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -1036,7 +1036,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Evd Evd() { - return new DenseEvd(this); + return DenseEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs index 66aedcc4..fe90924c 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { @@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization #else using System.Numerics; #endif - + /// /// Eigenvalues and eigenvectors of a complex matrix. /// @@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public sealed class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,13 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public DenseEvd(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -79,22 +74,28 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var order = matrix.RowCount; // Initialize matrices for eigenvalues and eigenvectors - EigenVectors = DenseMatrix.Identity(order); - D = matrix.CreateMatrix(order, order); - EigenValues = new DenseVector(order); + var eigenVectors = DenseMatrix.Identity(order); + var blockDiagonal = new DenseMatrix(order); + var eigenValues = new DenseVector(order); - IsSymmetric = true; + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); } } - Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) EigenVectors).Values, - ((DenseVector) EigenValues).Values, ((DenseMatrix) D).Values); + Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); + + return new DenseEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + DenseEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -586,7 +587,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var eps = Precision.DoubleMachinePrecision; double norm; - Complex x, y, z, exshift =Complex.Zero; + Complex x, y, z, exshift = Complex.Zero; // Outer loop over eigenvalue index var iter = 0; @@ -822,17 +823,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -899,16 +889,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs index f4b2b4d0..05f99df7 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -28,7 +32,8 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { - #if NOSYSNUMERICS + +#if NOSYSNUMERICS using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; @@ -51,6 +56,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// public abstract class Evd : Evd { + protected Evd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { + } + /// /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs index 0098db2b..6ce6b342 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { @@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public sealed class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,13 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public UserEvd(Matrix matrix) + public static UserEvd Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -79,21 +74,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var order = matrix.RowCount; // Initialize matricies for eigenvalues and eigenvectors - EigenVectors = DenseMatrix.Identity(order); - D = matrix.CreateMatrix(order, order); - EigenValues = new DenseVector(order); - - IsSymmetric = true; + var eigenVectors = DenseMatrix.Identity(order); + var blockDiagonal = matrix.CreateMatrix(order, order); + var eigenValues = new DenseVector(order); + + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); } } - if (IsSymmetric) + if (isSymmetric) { var matrixCopy = matrix.ToArray(); var tau = new Complex[order]; @@ -101,22 +96,29 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var e = new double[order]; SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - SymmetricDiagonalize(d, e, order); - SymmetricUntridiagonalize(matrixCopy, tau, order); + SymmetricDiagonalize(eigenVectors, d, e, order); + SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order); for (var i = 0; i < order; i++) { - EigenValues[i] = new Complex(d[i], e[i]); + eigenValues[i] = new Complex(d[i], e[i]); } } else { var matrixH = matrix.ToArray(); - NonsymmetricReduceToHessenberg(matrixH, order); - NonsymmetricReduceHessenberToRealSchur(matrixH, order); + NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order); } - D.SetDiagonal(EigenValues); + blockDiagonal.SetDiagonal(eigenValues); + + return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + UserEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -131,7 +133,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private static void SymmetricTridiagonalize(Complex[,] matrixA, double[] d, double[] e, Complex[] tau, int order) + static void SymmetricTridiagonalize(Complex[,] matrixA, double[] d, double[] e, Complex[] tau, int order) { double hh; tau[order - 1] = Complex.One; @@ -167,15 +169,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } Complex g = Math.Sqrt(h); - e[i] = scale * g.Real; + e[i] = scale*g.Real; Complex temp; var f = matrixA[i, i - 1]; if (f.Magnitude != 0) { - temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude; - h += f.Magnitude * g.Real; - g = 1.0 + (g / f.Magnitude); + temp = -(matrixA[i, i - 1].Conjugate()*tau[i].Conjugate())/f.Magnitude; + h += f.Magnitude*g.Real; + g = 1.0 + (g/f.Magnitude); matrixA[i, i - 1] *= g; } else @@ -194,31 +196,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // Form element of A*U. for (var k = 0; k <= j; k++) { - tmp += matrixA[j, k] * matrixA[i, k].Conjugate(); + tmp += matrixA[j, k]*matrixA[i, k].Conjugate(); } for (var k = j + 1; k <= i - 1; k++) { - tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate(); + tmp += matrixA[k, j].Conjugate()*matrixA[i, k].Conjugate(); } // Form element of P - tau[j] = tmp / h; - f += (tmp / h) * matrixA[i, j]; + tau[j] = tmp/h; + f += (tmp/h)*matrixA[i, j]; } - hh = f.Real / (h + h); + hh = f.Real/(h + h); // Form the reduced A. for (var j = 0; j < i; j++) { f = matrixA[i, j].Conjugate(); - g = tau[j] - (hh * f); + g = tau[j] - (hh*f); tau[j] = g.Conjugate(); for (var k = 0; k <= j; k++) { - matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]); + matrixA[j, k] -= (f*tau[k]) + (g*matrixA[i, k]); } } } @@ -233,7 +235,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization hh = d[i]; d[i] = matrixA[i, i].Real; - matrixA[i, i] = new Complex(hh, scale * Math.Sqrt(h)); + matrixA[i, i] = new Complex(hh, scale*Math.Sqrt(h)); } hh = d[0]; @@ -253,7 +255,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. /// - private void SymmetricDiagonalize(double[] d, double[] e, int order) + static void SymmetricDiagonalize(Matrix eigenVectors, double[] d, double[] e, int order) { const int maxiter = 1000; @@ -274,7 +276,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var m = l; while (m < order) { - if (Math.Abs(e[m]) <= eps * tst1) + if (Math.Abs(e[m]) <= eps*tst1) { break; } @@ -293,15 +295,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // Compute implicit shift var g = d[l]; - var p = (d[l + 1] - g) / (2.0 * e[l]); + var p = (d[l + 1] - g)/(2.0*e[l]); var r = SpecialFunctions.Hypotenuse(p, 1.0); if (p < 0) { r = -r; } - d[l] = e[l] / (p + r); - d[l + 1] = e[l] * (p + r); + d[l] = e[l]/(p + r); + d[l + 1] = e[l]*(p + r); var dl1 = d[l + 1]; var h = g - d[l]; @@ -325,27 +327,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization c3 = c2; c2 = c; s2 = s; - g = c * e[i]; - h = c * p; + g = c*e[i]; + h = c*p; r = SpecialFunctions.Hypotenuse(p, e[i]); - e[i + 1] = s * r; - s = e[i] / r; - c = p / r; - p = (c * d[i]) - (s * g); - d[i + 1] = h + (s * ((c * g) + (s * d[i]))); + e[i + 1] = s*r; + s = e[i]/r; + c = p/r; + p = (c*d[i]) - (s*g); + d[i + 1] = h + (s*((c*g) + (s*d[i]))); // Accumulate transformation. for (var k = 0; k < order; k++) { - h = EigenVectors.At(k, i + 1).Real; - EigenVectors.At(k, i + 1, (s * EigenVectors.At(k, i).Real) + (c * h)); - EigenVectors.At(k, i, (c * EigenVectors.At(k, i).Real) - (s * h)); + h = eigenVectors.At(k, i + 1).Real; + eigenVectors.At(k, i + 1, (s*eigenVectors.At(k, i).Real) + (c*h)); + eigenVectors.At(k, i, (c*eigenVectors.At(k, i).Real) - (s*h)); } } - p = (-s) * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; + p = (-s)*s2*c3*el1*e[l]/dl1; + e[l] = s*p; + d[l] = c*p; // Check for convergence. If too many iterations have been performed, // throw exception that Convergence Failed @@ -353,8 +355,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { throw new NonConvergenceException(); } - } - while (Math.Abs(e[l]) > eps * tst1); + } while (Math.Abs(e[l]) > eps*tst1); } d[l] = d[l] + f; @@ -381,9 +382,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization d[i] = p; for (var j = 0; j < order; j++) { - p = EigenVectors.At(j, i).Real; - EigenVectors.At(j, i, EigenVectors.At(j, k)); - EigenVectors.At(j, k, p); + p = eigenVectors.At(j, i).Real; + eigenVectors.At(j, i, eigenVectors.At(j, k)); + eigenVectors.At(j, k, p); } } } @@ -399,13 +400,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void SymmetricUntridiagonalize(Complex[,] matrixA, Complex[] tau, int order) + static void SymmetricUntridiagonalize(Matrix eigenVectors, Complex[,] matrixA, Complex[] tau, int order) { for (var i = 0; i < order; i++) { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, EigenVectors.At(i, j).Real * tau[i].Conjugate()); + eigenVectors.At(i, j, eigenVectors.At(i, j).Real*tau[i].Conjugate()); } } @@ -420,14 +421,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var s = Complex.Zero; for (var k = 0; k < i; k++) { - s += EigenVectors.At(k, j) * matrixA[i, k]; + s += eigenVectors.At(k, j)*matrixA[i, k]; } - s = (s / h) / h; + s = (s/h)/h; for (var k = 0; k < i; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - s * matrixA[i, k].Conjugate()); + eigenVectors.At(k, j, eigenVectors.At(k, j) - s*matrixA[i, k].Conjugate()); } } } @@ -443,7 +444,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutines in EISPACK. - private void NonsymmetricReduceToHessenberg(Complex[,] matrixH, int order) + static void NonsymmetricReduceToHessenberg(Matrix eigenVectors, Complex[,] matrixH, int order) { var ort = new Complex[order]; @@ -462,16 +463,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var h = 0.0; for (var i = order - 1; i >= m; i--) { - ort[i] = matrixH[i, m - 1] / scale; + ort[i] = matrixH[i, m - 1]/scale; h += ort[i].MagnitudeSquared(); } var g = Math.Sqrt(h); if (ort[m].Magnitude != 0) { - h = h + (ort[m].Magnitude * g); + h = h + (ort[m].Magnitude*g); g /= ort[m].Magnitude; - ort[m] = (1.0 + g) * ort[m]; + ort[m] = (1.0 + g)*ort[m]; } else { @@ -486,13 +487,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var f = Complex.Zero; for (var i = order - 1; i >= m; i--) { - f += ort[i].Conjugate() * matrixH[i, j]; + f += ort[i].Conjugate()*matrixH[i, j]; } - f = f / h; + f = f/h; for (var i = m; i < order; i++) { - matrixH[i, j] -= f * ort[i]; + matrixH[i, j] -= f*ort[i]; } } @@ -501,17 +502,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var f = Complex.Zero; for (var j = order - 1; j >= m; j--) { - f += ort[j] * matrixH[i, j]; + f += ort[j]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var j = m; j < order; j++) { - matrixH[i, j] -= f * ort[j].Conjugate(); + matrixH[i, j] -= f*ort[j].Conjugate(); } } - ort[m] = scale * ort[m]; + ort[m] = scale*ort[m]; matrixH[m, m - 1] *= -g; } } @@ -521,7 +522,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, i == j ? Complex.One : Complex.Zero); + eigenVectors.At(i, j, i == j ? Complex.One : Complex.Zero); } } @@ -529,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { if (matrixH[m, m - 1] != Complex.Zero && ort[m] != Complex.Zero) { - var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary); + var norm = (matrixH[m, m - 1].Real*ort[m].Real) + (matrixH[m, m - 1].Imaginary*ort[m].Imaginary); for (var i = m + 1; i < order; i++) { @@ -541,25 +542,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var g = Complex.Zero; for (var i = m; i < order; i++) { - g += ort[i].Conjugate() * EigenVectors.At(i, j); + g += ort[i].Conjugate()*eigenVectors.At(i, j); } // Double division avoids possible underflow g /= norm; for (var i = m; i < order; i++) { - EigenVectors.At(i, j, EigenVectors.At(i, j) + g * ort[i]); + eigenVectors.At(i, j, eigenVectors.At(i, j) + g*ort[i]); } } } } - + // Create real subdiagonal elements. for (var i = 1; i < order; i++) { if (matrixH[i, i - 1].Imaginary != 0.0) { - var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude; + var y = matrixH[i, i - 1]/matrixH[i, i - 1].Magnitude; matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude; for (var j = i; j < order; j++) { @@ -573,7 +574,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (var j = 0; j < order; j++) { - EigenVectors.At(j, i, EigenVectors.At(j, i) * y); + eigenVectors.At(j, i, eigenVectors.At(j, i)*y); } } } @@ -588,7 +589,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void NonsymmetricReduceHessenberToRealSchur(Complex[,] matrixH, int order) + static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVectors, Vector eigenValues, Complex[,] matrixH, int order) { // Initialize var n = order - 1; @@ -606,7 +607,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization while (l > 0) { var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary); - if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1) + if (Math.Abs(matrixH[l, l - 1].Real) < eps*tst1) { break; } @@ -619,7 +620,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization if (l == n) { matrixH[n, n] += exshift; - EigenValues[n] = matrixH[n, n]; + eigenValues[n] = matrixH[n, n]; n--; iter = 0; } @@ -630,18 +631,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization if (iter != 10 && iter != 20) { s = matrixH[n, n]; - x = matrixH[n - 1, n] * matrixH[n, n - 1].Real; + x = matrixH[n - 1, n]*matrixH[n, n - 1].Real; if (x.Real != 0.0 || x.Imaginary != 0.0) { - y = (matrixH[n - 1, n - 1] - s) / 2.0; - z = ((y * y) + x).SquareRoot(); - if ((y.Real * z.Real) + (y.Imaginary * z.Imaginary) < 0.0) + y = (matrixH[n - 1, n - 1] - s)/2.0; + z = ((y*y) + x).SquareRoot(); + if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0) { z *= -1.0; } - x /= y + z; + x /= y + z; s = s - x; } } @@ -664,17 +665,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { s = matrixH[i, i - 1].Real; norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real); - x = matrixH[i - 1, i - 1] / norm; - EigenValues[i - 1] = x; + x = matrixH[i - 1, i - 1]/norm; + eigenValues[i - 1] = x; matrixH[i - 1, i - 1] = norm; - matrixH[i, i - 1] = new Complex(0.0, s.Real / norm); + matrixH[i, i - 1] = new Complex(0.0, s.Real/norm); for (var j = i; j < order; j++) { y = matrixH[i - 1, j]; z = matrixH[i, j]; - matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z); - matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y); + matrixH[i - 1, j] = (x.Conjugate()*y) + (matrixH[i, i - 1].Imaginary*z); + matrixH[i, j] = (x*z) - (matrixH[i, i - 1].Imaginary*y); } } @@ -693,30 +694,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // Inverse operation (columns). for (var j = l + 1; j <= n; j++) { - x = EigenValues[j - 1]; + x = eigenValues[j - 1]; for (var i = 0; i <= j; i++) { z = matrixH[i, j]; if (i != j) { y = matrixH[i, j - 1]; - matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z); + matrixH[i, j - 1] = (x*y) + (matrixH[j, j - 1].Imaginary*z); } else { y = matrixH[i, j - 1].Real; - matrixH[i, j - 1] = new Complex((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary); + matrixH[i, j - 1] = new Complex((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[j, j - 1].Imaginary*z.Real), matrixH[i, j - 1].Imaginary); } - matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); + matrixH[i, j] = (x.Conjugate()*z) - (matrixH[j, j - 1].Imaginary*y); } for (var i = 0; i < order; i++) { - y = EigenVectors.At(i, j - 1); - z = EigenVectors.At(i, j); - EigenVectors.At(i, j - 1, (x * y) + (matrixH[j, j - 1].Imaginary * z)); - EigenVectors.At(i, j, (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y)); + y = eigenVectors.At(i, j - 1); + z = eigenVectors.At(i, j); + eigenVectors.At(i, j - 1, (x*y) + (matrixH[j, j - 1].Imaginary*z)); + eigenVectors.At(i, j, (x.Conjugate()*z) - (matrixH[j, j - 1].Imaginary*y)); } } @@ -729,7 +730,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (var i = 0; i < order; i++) { - EigenVectors.At(i, n, EigenVectors.At(i, n) * s); + eigenVectors.At(i, n, eigenVectors.At(i, n)*s); } } } @@ -758,7 +759,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (n = order - 1; n > 0; n--) { - x = EigenValues[n]; + x = eigenValues[n]; matrixH[n, n] = 1.0; for (var i = n - 1; i >= 0; i--) @@ -766,24 +767,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization z = 0.0; for (var j = i + 1; j <= n; j++) { - z += matrixH[i, j] * matrixH[j, n]; + z += matrixH[i, j]*matrixH[j, n]; } - y = x - EigenValues[i]; + y = x - eigenValues[i]; if (y.Real == 0.0 && y.Imaginary == 0.0) { - y = eps * norm; + y = eps*norm; } - matrixH[i, n] = z / y; + matrixH[i, n] = z/y; // Overflow control var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary); - if ((eps * tr) * tr > 1) + if ((eps*tr)*tr > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n] = matrixH[j, n] / tr; + matrixH[j, n] = matrixH[j, n]/tr; } } } @@ -797,14 +798,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization z = Complex.Zero; for (var k = 0; k <= j; k++) { - z += EigenVectors.At(i, k) * matrixH[k, j]; + z += eigenVectors.At(i, k)*matrixH[k, j]; } - EigenVectors.At(i, j, z); + eigenVectors.At(i, j, z); } } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -812,17 +813,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -855,7 +845,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j).Conjugate() * input.At(i, k); + value += EigenVectors.At(i, j).Conjugate()*input.At(i, k); } value /= EigenValues[j].Real; @@ -869,7 +859,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization Complex value = 0.0; for (var i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result.At(j, k, value); @@ -878,7 +868,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } else { - throw new ArgumentException(Resources.ArgumentMatrixSymmetric); + throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } @@ -889,16 +879,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) @@ -926,7 +906,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j).Conjugate() * input[i]; + value += EigenVectors.At(i, j).Conjugate()*input[i]; } value /= EigenValues[j].Real; @@ -940,7 +920,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization value = 0; for (int i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result[j] = value; @@ -952,4 +932,4 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } } } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index 4095c001..96b1090d 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -485,7 +485,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Evd Evd() { - return new UserEvd(this); + return UserEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index 12a7c0bb..8c36d9cf 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -1031,7 +1031,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Evd Evd() { - return new DenseEvd(this); + return DenseEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs index 29bfeae3..7accef78 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs @@ -28,11 +28,19 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { + using Numerics; + +#if NOSYSNUMERICS + using Complex = Numerics.Complex; +#else + using Complex = System.Numerics.Complex; +#endif + /// /// Eigenvalues and eigenvectors of a complex matrix. /// @@ -48,7 +56,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public sealed class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,13 +65,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public DenseEvd(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -72,22 +75,28 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var order = matrix.RowCount; // Initialize matrices for eigenvalues and eigenvectors - EigenVectors = DenseMatrix.Identity(order); - D = matrix.CreateMatrix(order, order); - EigenValues = new Complex.DenseVector(order); + var eigenVectors = DenseMatrix.Identity(order); + var blockDiagonal = new DenseMatrix(order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - IsSymmetric = true; + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); } } - Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) EigenVectors).Values, - ((Complex.DenseVector)EigenValues).Values, ((DenseMatrix)D).Values); + Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); + + return new DenseEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + DenseEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -815,17 +824,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -861,7 +859,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization value += ((DenseMatrix) EigenVectors).Values[(j*order) + i].Conjugate()*input.At(i, k); } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; @@ -892,16 +890,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) @@ -932,7 +920,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization value += ((DenseMatrix) EigenVectors).Values[(j*order) + i].Conjugate()*input[i]; } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs index 2e6cb373..3658c3ef 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -24,11 +28,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { - using System; using Numerics; #if NOSYSNUMERICS @@ -54,6 +58,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// public abstract class Evd : Evd { + protected Evd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { + } + /// /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. /// @@ -66,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { det *= EigenValues[i]; - if (((Complex32)EigenValues[i]).AlmostEqual(Complex32.Zero)) + if (((Complex32) EigenValues[i]).AlmostEqual(Complex32.Zero)) { return 0; } @@ -87,7 +96,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var rank = 0; for (var i = 0; i < EigenValues.Count; i++) { - if (((Complex32)EigenValues[i]).AlmostEqual(Complex32.Zero)) + if (((Complex32) EigenValues[i]).AlmostEqual(Complex32.Zero)) { continue; } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs index c4e99684..e52383d2 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { @@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public sealed class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -63,13 +63,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public UserEvd(Matrix matrix) + public static UserEvd Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -78,21 +73,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var order = matrix.RowCount; // Initialize matricies for eigenvalues and eigenvectors - EigenVectors = DenseMatrix.Identity(order); - D = matrix.CreateMatrix(order, order); - EigenValues = new LinearAlgebra.Complex.DenseVector(order); - - IsSymmetric = true; + var eigenVectors = DenseMatrix.Identity(order); + var blockDiagonal = matrix.CreateMatrix(order, order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); + + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); } } - if (IsSymmetric) + if (isSymmetric) { var matrixCopy = matrix.ToArray(); var tau = new Complex32[order]; @@ -100,25 +95,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var e = new float[order]; SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - SymmetricDiagonalize(d, e, order); - SymmetricUntridiagonalize(matrixCopy, tau, order); + SymmetricDiagonalize(eigenVectors, d, e, order); + SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order); for (var i = 0; i < order; i++) { - EigenValues[i] = new Complex(d[i], e[i]); + eigenValues[i] = new Complex(d[i], e[i]); } } else { var matrixH = matrix.ToArray(); - NonsymmetricReduceToHessenberg(matrixH, order); - NonsymmetricReduceHessenberToRealSchur(matrixH, order); + NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order); } - for (var i = 0; i < EigenValues.Count; i++) + for (var i = 0; i < eigenValues.Count; i++) { - D.At(i, i, (Complex32)EigenValues[i]); + blockDiagonal.At(i, i, (Complex32) eigenValues[i]); } + + return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + UserEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -133,7 +135,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private static void SymmetricTridiagonalize(Complex32[,] matrixA, float[] d, float[] e, Complex32[] tau, int order) + static void SymmetricTridiagonalize(Complex32[,] matrixA, float[] d, float[] e, Complex32[] tau, int order) { float hh; tau[order - 1] = Complex32.One; @@ -168,16 +170,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization h += matrixA[i, k].MagnitudeSquared; } - Complex32 g = (float)Math.Sqrt(h); - e[i] = scale * g.Real; + Complex32 g = (float) Math.Sqrt(h); + e[i] = scale*g.Real; Complex32 temp; var f = matrixA[i, i - 1]; if (f.Magnitude != 0) { - temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude; - h += f.Magnitude * g.Real; - g = 1.0f + (g / f.Magnitude); + temp = -(matrixA[i, i - 1].Conjugate()*tau[i].Conjugate())/f.Magnitude; + h += f.Magnitude*g.Real; + g = 1.0f + (g/f.Magnitude); matrixA[i, i - 1] *= g; } else @@ -196,31 +198,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // Form element of A*U. for (var k = 0; k <= j; k++) { - tmp += matrixA[j, k] * matrixA[i, k].Conjugate(); + tmp += matrixA[j, k]*matrixA[i, k].Conjugate(); } for (var k = j + 1; k <= i - 1; k++) { - tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate(); + tmp += matrixA[k, j].Conjugate()*matrixA[i, k].Conjugate(); } // Form element of P - tau[j] = tmp / h; - f += (tmp / h) * matrixA[i, j]; + tau[j] = tmp/h; + f += (tmp/h)*matrixA[i, j]; } - hh = f.Real / (h + h); + hh = f.Real/(h + h); // Form the reduced A. for (var j = 0; j < i; j++) { f = matrixA[i, j].Conjugate(); - g = tau[j] - (hh * f); + g = tau[j] - (hh*f); tau[j] = g.Conjugate(); for (var k = 0; k <= j; k++) { - matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]); + matrixA[j, k] -= (f*tau[k]) + (g*matrixA[i, k]); } } } @@ -235,7 +237,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization hh = d[i]; d[i] = matrixA[i, i].Real; - matrixA[i, i] = new Complex32(hh, scale * (float)Math.Sqrt(h)); + matrixA[i, i] = new Complex32(hh, scale*(float) Math.Sqrt(h)); } hh = d[0]; @@ -255,7 +257,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. /// - private void SymmetricDiagonalize(float[] d, float[] e, int order) + static void SymmetricDiagonalize(Matrix eigenVectors, float[] d, float[] e, int order) { const int maxiter = 1000; @@ -276,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var m = l; while (m < order) { - if (Math.Abs(e[m]) <= eps * tst1) + if (Math.Abs(e[m]) <= eps*tst1) { break; } @@ -295,15 +297,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // Compute implicit shift var g = d[l]; - var p = (d[l + 1] - g) / (2.0f * e[l]); + var p = (d[l + 1] - g)/(2.0f*e[l]); var r = SpecialFunctions.Hypotenuse(p, 1.0f); if (p < 0) { r = -r; } - d[l] = e[l] / (p + r); - d[l + 1] = e[l] * (p + r); + d[l] = e[l]/(p + r); + d[l + 1] = e[l]*(p + r); var dl1 = d[l + 1]; var h = g - d[l]; @@ -327,27 +329,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization c3 = c2; c2 = c; s2 = s; - g = c * e[i]; - h = c * p; + g = c*e[i]; + h = c*p; r = SpecialFunctions.Hypotenuse(p, e[i]); - e[i + 1] = s * r; - s = e[i] / r; - c = p / r; - p = (c * d[i]) - (s * g); - d[i + 1] = h + (s * ((c * g) + (s * d[i]))); + e[i + 1] = s*r; + s = e[i]/r; + c = p/r; + p = (c*d[i]) - (s*g); + d[i + 1] = h + (s*((c*g) + (s*d[i]))); // Accumulate transformation. for (var k = 0; k < order; k++) { - h = EigenVectors.At(k, i + 1).Real; - EigenVectors.At(k, i + 1, (s * EigenVectors.At(k, i).Real) + (c * h)); - EigenVectors.At(k, i, (c * EigenVectors.At(k, i).Real) - (s * h)); + h = eigenVectors.At(k, i + 1).Real; + eigenVectors.At(k, i + 1, (s*eigenVectors.At(k, i).Real) + (c*h)); + eigenVectors.At(k, i, (c*eigenVectors.At(k, i).Real) - (s*h)); } } - p = (-s) * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; + p = (-s)*s2*c3*el1*e[l]/dl1; + e[l] = s*p; + d[l] = c*p; // Check for convergence. If too many iterations have been performed, // throw exception that Convergence Failed @@ -355,8 +357,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { throw new NonConvergenceException(); } - } - while (Math.Abs(e[l]) > eps * tst1); + } while (Math.Abs(e[l]) > eps*tst1); } d[l] = d[l] + f; @@ -383,9 +384,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization d[i] = p; for (var j = 0; j < order; j++) { - p = EigenVectors.At(j, i).Real; - EigenVectors.At(j, i, EigenVectors.At(j, k)); - EigenVectors.At(j, k, p); + p = eigenVectors.At(j, i).Real; + eigenVectors.At(j, i, eigenVectors.At(j, k)); + eigenVectors.At(j, k, p); } } } @@ -401,13 +402,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void SymmetricUntridiagonalize(Complex32[,] matrixA, Complex32[] tau, int order) + static void SymmetricUntridiagonalize(Matrix eigenVectors, Complex32[,] matrixA, Complex32[] tau, int order) { for (var i = 0; i < order; i++) { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, EigenVectors.At(i, j).Real * tau[i].Conjugate()); + eigenVectors.At(i, j, eigenVectors.At(i, j).Real*tau[i].Conjugate()); } } @@ -422,14 +423,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var s = Complex32.Zero; for (var k = 0; k < i; k++) { - s += EigenVectors.At(k, j) * matrixA[i, k]; + s += eigenVectors.At(k, j)*matrixA[i, k]; } - s = (s / h) / h; + s = (s/h)/h; for (var k = 0; k < i; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - s * matrixA[i, k].Conjugate()); + eigenVectors.At(k, j, eigenVectors.At(k, j) - s*matrixA[i, k].Conjugate()); } } } @@ -445,7 +446,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutines in EISPACK. - private void NonsymmetricReduceToHessenberg(Complex32[,] matrixH, int order) + static void NonsymmetricReduceToHessenberg(Matrix eigenVectors, Complex32[,] matrixH, int order) { var ort = new Complex32[order]; @@ -464,16 +465,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var h = 0.0f; for (var i = order - 1; i >= m; i--) { - ort[i] = matrixH[i, m - 1] / scale; + ort[i] = matrixH[i, m - 1]/scale; h += ort[i].MagnitudeSquared; } - var g = (float)Math.Sqrt(h); + var g = (float) Math.Sqrt(h); if (ort[m].Magnitude != 0) { - h = h + (ort[m].Magnitude * g); + h = h + (ort[m].Magnitude*g); g /= ort[m].Magnitude; - ort[m] = (1.0f + g) * ort[m]; + ort[m] = (1.0f + g)*ort[m]; } else { @@ -488,13 +489,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var f = Complex32.Zero; for (var i = order - 1; i >= m; i--) { - f += ort[i].Conjugate() * matrixH[i, j]; + f += ort[i].Conjugate()*matrixH[i, j]; } - f = f / h; + f = f/h; for (var i = m; i < order; i++) { - matrixH[i, j] -= f * ort[i]; + matrixH[i, j] -= f*ort[i]; } } @@ -503,17 +504,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var f = Complex32.Zero; for (var j = order - 1; j >= m; j--) { - f += ort[j] * matrixH[i, j]; + f += ort[j]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var j = m; j < order; j++) { - matrixH[i, j] -= f * ort[j].Conjugate(); + matrixH[i, j] -= f*ort[j].Conjugate(); } } - ort[m] = scale * ort[m]; + ort[m] = scale*ort[m]; matrixH[m, m - 1] *= -g; } } @@ -523,7 +524,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, i == j ? Complex32.One : Complex32.Zero); + eigenVectors.At(i, j, i == j ? Complex32.One : Complex32.Zero); } } @@ -531,7 +532,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { if (matrixH[m, m - 1] != Complex32.Zero && ort[m] != Complex32.Zero) { - var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary); + var norm = (matrixH[m, m - 1].Real*ort[m].Real) + (matrixH[m, m - 1].Imaginary*ort[m].Imaginary); for (var i = m + 1; i < order; i++) { @@ -543,25 +544,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var g = Complex32.Zero; for (var i = m; i < order; i++) { - g += ort[i].Conjugate() * EigenVectors.At(i, j); + g += ort[i].Conjugate()*eigenVectors.At(i, j); } // Double division avoids possible underflow g /= norm; for (var i = m; i < order; i++) { - EigenVectors.At(i, j, EigenVectors.At(i, j) + g * ort[i]); + eigenVectors.At(i, j, eigenVectors.At(i, j) + g*ort[i]); } } } } - + // Create real subdiagonal elements. for (var i = 1; i < order; i++) { if (matrixH[i, i - 1].Imaginary != 0.0f) { - var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude; + var y = matrixH[i, i - 1]/matrixH[i, i - 1].Magnitude; matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude; for (var j = i; j < order; j++) { @@ -575,7 +576,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (var j = 0; j < order; j++) { - EigenVectors.At(j, i, EigenVectors.At(j, i) * y); + eigenVectors.At(j, i, eigenVectors.At(j, i)*y); } } } @@ -590,11 +591,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void NonsymmetricReduceHessenberToRealSchur(Complex32[,] matrixH, int order) + static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVectors, Vector eigenValues, Complex32[,] matrixH, int order) { // Initialize var n = order - 1; - var eps = (float)Precision.SingleMachinePrecision; + var eps = (float) Precision.SingleMachinePrecision; float norm; Complex32 x, y, z, exshift = Complex32.Zero; @@ -608,7 +609,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization while (l > 0) { var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary); - if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1) + if (Math.Abs(matrixH[l, l - 1].Real) < eps*tst1) { break; } @@ -621,7 +622,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization if (l == n) { matrixH[n, n] += exshift; - EigenValues[n] = matrixH[n, n].ToComplex(); + eigenValues[n] = matrixH[n, n].ToComplex(); n--; iter = 0; } @@ -632,18 +633,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization if (iter != 10 && iter != 20) { s = matrixH[n, n]; - x = matrixH[n - 1, n] * matrixH[n, n - 1].Real; + x = matrixH[n - 1, n]*matrixH[n, n - 1].Real; if (x.Real != 0.0f || x.Imaginary != 0.0f) { - y = (matrixH[n - 1, n - 1] - s) / 2.0f; - z = ((y * y) + x).SquareRoot(); - if ((y.Real * z.Real) + (y.Imaginary * z.Imaginary) < 0.0f) + y = (matrixH[n - 1, n - 1] - s)/2.0f; + z = ((y*y) + x).SquareRoot(); + if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0f) { z *= -1.0f; } - x /= y + z; + x /= y + z; s = s - x; } } @@ -666,17 +667,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { s = matrixH[i, i - 1].Real; norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real); - x = matrixH[i - 1, i - 1] / norm; - EigenValues[i - 1] = x.ToComplex(); + x = matrixH[i - 1, i - 1]/norm; + eigenValues[i - 1] = x.ToComplex(); matrixH[i - 1, i - 1] = norm; - matrixH[i, i - 1] = new Complex32(0.0f, s.Real / norm); + matrixH[i, i - 1] = new Complex32(0.0f, s.Real/norm); for (var j = i; j < order; j++) { y = matrixH[i - 1, j]; z = matrixH[i, j]; - matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z); - matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y); + matrixH[i - 1, j] = (x.Conjugate()*y) + (matrixH[i, i - 1].Imaginary*z); + matrixH[i, j] = (x*z) - (matrixH[i, i - 1].Imaginary*y); } } @@ -695,30 +696,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // Inverse operation (columns). for (var j = l + 1; j <= n; j++) { - x = (Complex32)EigenValues[j - 1]; + x = (Complex32) eigenValues[j - 1]; for (var i = 0; i <= j; i++) { z = matrixH[i, j]; if (i != j) { y = matrixH[i, j - 1]; - matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z); + matrixH[i, j - 1] = (x*y) + (matrixH[j, j - 1].Imaginary*z); } else { y = matrixH[i, j - 1].Real; - matrixH[i, j - 1] = new Complex32((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary); + matrixH[i, j - 1] = new Complex32((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[j, j - 1].Imaginary*z.Real), matrixH[i, j - 1].Imaginary); } - matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); + matrixH[i, j] = (x.Conjugate()*z) - (matrixH[j, j - 1].Imaginary*y); } for (var i = 0; i < order; i++) { - y = EigenVectors.At(i, j - 1); - z = EigenVectors.At(i, j); - EigenVectors.At(i, j - 1, (x * y) + (matrixH[j, j - 1].Imaginary * z)); - EigenVectors.At(i, j, (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y)); + y = eigenVectors.At(i, j - 1); + z = eigenVectors.At(i, j); + eigenVectors.At(i, j - 1, (x*y) + (matrixH[j, j - 1].Imaginary*z)); + eigenVectors.At(i, j, (x.Conjugate()*z) - (matrixH[j, j - 1].Imaginary*y)); } } @@ -731,7 +732,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (var i = 0; i < order; i++) { - EigenVectors.At(i, n, EigenVectors.At(i, n) * s); + eigenVectors.At(i, n, eigenVectors.At(i, n)*s); } } } @@ -760,7 +761,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (n = order - 1; n > 0; n--) { - x = (Complex32)EigenValues[n]; + x = (Complex32) eigenValues[n]; matrixH[n, n] = 1.0f; for (var i = n - 1; i >= 0; i--) @@ -768,24 +769,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization z = 0.0f; for (var j = i + 1; j <= n; j++) { - z += matrixH[i, j] * matrixH[j, n]; + z += matrixH[i, j]*matrixH[j, n]; } - y = x - (Complex32)EigenValues[i]; + y = x - (Complex32) eigenValues[i]; if (y.Real == 0.0f && y.Imaginary == 0.0f) { - y = eps * norm; + y = eps*norm; } - matrixH[i, n] = z / y; + matrixH[i, n] = z/y; // Overflow control var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary); - if ((eps * tr) * tr > 1) + if ((eps*tr)*tr > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n] = matrixH[j, n] / tr; + matrixH[j, n] = matrixH[j, n]/tr; } } } @@ -799,14 +800,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization z = Complex32.Zero; for (var k = 0; k <= j; k++) { - z += EigenVectors.At(i, k) * matrixH[k, j]; + z += eigenVectors.At(i, k)*matrixH[k, j]; } - EigenVectors.At(i, j, z); + eigenVectors.At(i, j, z); } } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -814,17 +815,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -857,10 +847,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j).Conjugate() * input.At(i, k); + value += EigenVectors.At(i, j).Conjugate()*input.At(i, k); } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; @@ -871,7 +861,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization Complex32 value = 0.0f; for (var i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result.At(j, k, value); @@ -880,7 +870,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } else { - throw new ArgumentException(Resources.ArgumentMatrixSymmetric); + throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } @@ -891,16 +881,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) @@ -928,10 +908,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j).Conjugate() * input[i]; + value += EigenVectors.At(i, j).Conjugate()*input[i]; } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; @@ -942,7 +922,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization value = 0; for (int i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result[j] = value; @@ -954,4 +934,4 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } } } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index 0028ac73..4ba57d64 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -480,7 +480,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Evd Evd() { - return new UserEvd(this); + return UserEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index a628faaa..645650d3 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -1062,7 +1062,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Evd Evd() { - return new DenseEvd(this); + return DenseEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs index 59134b6a..d355f8c4 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { @@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public sealed class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,13 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public DenseEvd(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -79,22 +74,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var order = matrix.RowCount; // Initialize matrices for eigenvalues and eigenvectors - EigenVectors = matrix.CreateMatrix(order, order); - D = matrix.CreateMatrix(order, order); - EigenValues = new LinearAlgebra.Complex.DenseVector(order); + var eigenVectors = new DenseMatrix(order); + var blockDiagonal = new DenseMatrix(order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - IsSymmetric = true; + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i); } } - Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) EigenVectors).Values, - ((LinearAlgebra.Complex.DenseVector)EigenValues).Values, ((DenseMatrix)D).Values); + Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); + + return new DenseEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + DenseEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -1114,17 +1115,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -1191,16 +1181,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs index 49486dc7..82737f82 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -28,7 +32,8 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { - #if NOSYSNUMERICS + +#if NOSYSNUMERICS using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; @@ -51,6 +56,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// public abstract class Evd : Evd { + protected Evd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { + } + /// /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs index 3908643b..eaedffb1 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { @@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public sealed class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,13 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public UserEvd(Matrix matrix) + public static UserEvd Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -79,57 +74,64 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var order = matrix.RowCount; // Initialize matricies for eigenvalues and eigenvectors - EigenVectors = matrix.CreateMatrix(order, order); - D = matrix.CreateMatrix(order, order); - EigenValues = new LinearAlgebra.Complex.DenseVector(order); - - IsSymmetric = true; + var eigenVectors = matrix.CreateMatrix(order, order); + var blockDiagonal = matrix.CreateMatrix(order, order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); + + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i); } } var d = new double[order]; var e = new double[order]; - if (IsSymmetric) + if (isSymmetric) { - matrix.CopyTo(EigenVectors); - d = EigenVectors.Row(order - 1).ToArray(); + matrix.CopyTo(eigenVectors); + d = eigenVectors.Row(order - 1).ToArray(); - SymmetricTridiagonalize(d, e, order); - SymmetricDiagonalize(d, e, order); + SymmetricTridiagonalize(eigenVectors, d, e, order); + SymmetricDiagonalize(eigenVectors, d, e, order); } else { var matrixH = matrix.ToArray(); - NonsymmetricReduceToHessenberg(matrixH, order); - NonsymmetricReduceHessenberToRealSchur(matrixH, d, e, order); + NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(eigenVectors, matrixH, d, e, order); } for (var i = 0; i < order; i++) { - D.At(i, i, d[i]); + blockDiagonal.At(i, i, d[i]); if (e[i] > 0) { - D.At(i, i + 1, e[i]); + blockDiagonal.At(i, i + 1, e[i]); } else if (e[i] < 0) { - D.At(i, i - 1, e[i]); + blockDiagonal.At(i, i - 1, e[i]); } } for (var i = 0; i < order; i++) { - EigenValues[i] = new Complex(d[i], e[i]); + eigenValues[i] = new Complex(d[i], e[i]); } + + return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + UserEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -142,7 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void SymmetricTridiagonalize(double[] d, double[] e, int order) + static void SymmetricTridiagonalize(Matrix eigenVectors, double[] d, double[] e, int order) { // Householder reduction to tridiagonal form. for (var i = order - 1; i > 0; i--) @@ -161,9 +163,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization e[i] = d[i - 1]; for (var j = 0; j < i; j++) { - d[j] = EigenVectors.At(i - 1, j); - EigenVectors.At(i, j, 0.0); - EigenVectors.At(j, i, 0.0); + d[j] = eigenVectors.At(i - 1, j); + eigenVectors.At(i, j, 0.0); + eigenVectors.At(j, i, 0.0); } } else @@ -172,7 +174,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var k = 0; k < i; k++) { d[k] /= scale; - h += d[k] * d[k]; + h += d[k]*d[k]; } var f = d[i - 1]; @@ -182,8 +184,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization g = -g; } - e[i] = scale * g; - h = h - (f * g); + e[i] = scale*g; + h = h - (f*g); d[i - 1] = f - g; for (var j = 0; j < i; j++) @@ -195,13 +197,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var j = 0; j < i; j++) { f = d[j]; - EigenVectors.At(j, i, f); - g = e[j] + (EigenVectors.At(j, j) * f); + eigenVectors.At(j, i, f); + g = e[j] + (eigenVectors.At(j, j)*f); for (var k = j + 1; k <= i - 1; k++) { - g += EigenVectors.At(k, j) * d[k]; - e[k] += EigenVectors.At(k, j) * f; + g += eigenVectors.At(k, j)*d[k]; + e[k] += eigenVectors.At(k, j)*f; } e[j] = g; @@ -212,14 +214,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var j = 0; j < i; j++) { e[j] /= h; - f += e[j] * d[j]; + f += e[j]*d[j]; } - var hh = f / (h + h); + var hh = f/(h + h); for (var j = 0; j < i; j++) { - e[j] -= hh * d[j]; + e[j] -= hh*d[j]; } for (var j = 0; j < i; j++) @@ -229,11 +231,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var k = j; k <= i - 1; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - (f * e[k]) - (g * d[k])); + eigenVectors.At(k, j, eigenVectors.At(k, j) - (f*e[k]) - (g*d[k])); } - d[j] = EigenVectors.At(i - 1, j); - EigenVectors.At(i, j, 0.0); + d[j] = eigenVectors.At(i - 1, j); + eigenVectors.At(i, j, 0.0); } } @@ -243,14 +245,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Accumulate transformations. for (var i = 0; i < order - 1; i++) { - EigenVectors.At(order - 1, i, EigenVectors.At(i, i)); - EigenVectors.At(i, i, 1.0); + eigenVectors.At(order - 1, i, eigenVectors.At(i, i)); + eigenVectors.At(i, i, 1.0); var h = d[i + 1]; if (h != 0.0) { for (var k = 0; k <= i; k++) { - d[k] = EigenVectors.At(k, i + 1) / h; + d[k] = eigenVectors.At(k, i + 1)/h; } for (var j = 0; j <= i; j++) @@ -258,29 +260,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var g = 0.0; for (var k = 0; k <= i; k++) { - g += EigenVectors.At(k, i + 1) * EigenVectors.At(k, j); + g += eigenVectors.At(k, i + 1)*eigenVectors.At(k, j); } for (var k = 0; k <= i; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - g * d[k]); + eigenVectors.At(k, j, eigenVectors.At(k, j) - g*d[k]); } } } for (var k = 0; k <= i; k++) { - EigenVectors.At(k, i + 1, 0.0); + eigenVectors.At(k, i + 1, 0.0); } } for (var j = 0; j < order; j++) { - d[j] = EigenVectors.At(order - 1, j); - EigenVectors.At(order - 1, j, 0.0); + d[j] = eigenVectors.At(order - 1, j); + eigenVectors.At(order - 1, j, 0.0); } - EigenVectors.At(order - 1, order - 1, 1.0); + eigenVectors.At(order - 1, order - 1, 1.0); e[0] = 0.0; } @@ -295,7 +297,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. /// - private void SymmetricDiagonalize(double[] d, double[] e, int order) + static void SymmetricDiagonalize(Matrix eigenVectors, double[] d, double[] e, int order) { const int maxiter = 1000; @@ -316,7 +318,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var m = l; while (m < order) { - if (Math.Abs(e[m]) <= eps * tst1) + if (Math.Abs(e[m]) <= eps*tst1) { break; } @@ -335,15 +337,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Compute implicit shift var g = d[l]; - var p = (d[l + 1] - g) / (2.0 * e[l]); + var p = (d[l + 1] - g)/(2.0*e[l]); var r = SpecialFunctions.Hypotenuse(p, 1.0); if (p < 0) { r = -r; } - d[l] = e[l] / (p + r); - d[l + 1] = e[l] * (p + r); + d[l] = e[l]/(p + r); + d[l + 1] = e[l]*(p + r); var dl1 = d[l + 1]; var h = g - d[l]; @@ -367,27 +369,27 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization c3 = c2; c2 = c; s2 = s; - g = c * e[i]; - h = c * p; + g = c*e[i]; + h = c*p; r = SpecialFunctions.Hypotenuse(p, e[i]); - e[i + 1] = s * r; - s = e[i] / r; - c = p / r; - p = (c * d[i]) - (s * g); - d[i + 1] = h + (s * ((c * g) + (s * d[i]))); + e[i + 1] = s*r; + s = e[i]/r; + c = p/r; + p = (c*d[i]) - (s*g); + d[i + 1] = h + (s*((c*g) + (s*d[i]))); // Accumulate transformation. for (var k = 0; k < order; k++) { - h = EigenVectors.At(k, i + 1); - EigenVectors.At(k, i + 1, (s * EigenVectors.At(k, i)) + (c * h)); - EigenVectors.At(k, i, (c * EigenVectors.At(k, i)) - (s * h)); + h = eigenVectors.At(k, i + 1); + eigenVectors.At(k, i + 1, (s*eigenVectors.At(k, i)) + (c*h)); + eigenVectors.At(k, i, (c*eigenVectors.At(k, i)) - (s*h)); } } - p = (-s) * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; + p = (-s)*s2*c3*el1*e[l]/dl1; + e[l] = s*p; + d[l] = c*p; // Check for convergence. If too many iterations have been performed, // throw exception that Convergence Failed @@ -395,8 +397,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { throw new NonConvergenceException(); } - } - while (Math.Abs(e[l]) > eps * tst1); + } while (Math.Abs(e[l]) > eps*tst1); } d[l] = d[l] + f; @@ -423,9 +424,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization d[i] = p; for (var j = 0; j < order; j++) { - p = EigenVectors.At(j, i); - EigenVectors.At(j, i, EigenVectors.At(j, k)); - EigenVectors.At(j, k, p); + p = eigenVectors.At(j, i); + eigenVectors.At(j, i, eigenVectors.At(j, k)); + eigenVectors.At(j, k, p); } } } @@ -440,7 +441,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutines in EISPACK. - private void NonsymmetricReduceToHessenberg(double[,] matrixH, int order) + static void NonsymmetricReduceToHessenberg(Matrix eigenVectors, double[,] matrixH, int order) { var ort = new double[order]; @@ -459,8 +460,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var h = 0.0; for (var i = order - 1; i >= m; i--) { - ort[i] = matrixH[i, m - 1] / scale; - h += ort[i] * ort[i]; + ort[i] = matrixH[i, m - 1]/scale; + h += ort[i]*ort[i]; } var g = Math.Sqrt(h); @@ -469,7 +470,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization g = -g; } - h = h - (ort[m] * g); + h = h - (ort[m]*g); ort[m] = ort[m] - g; // Apply Householder similarity transformation @@ -479,13 +480,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var f = 0.0; for (var i = order - 1; i >= m; i--) { - f += ort[i] * matrixH[i, j]; + f += ort[i]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var i = m; i < order; i++) { - matrixH[i, j] -= f * ort[i]; + matrixH[i, j] -= f*ort[i]; } } @@ -494,18 +495,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var f = 0.0; for (var j = order - 1; j >= m; j--) { - f += ort[j] * matrixH[i, j]; + f += ort[j]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var j = m; j < order; j++) { - matrixH[i, j] -= f * ort[j]; + matrixH[i, j] -= f*ort[j]; } } - ort[m] = scale * ort[m]; - matrixH[m, m - 1] = scale * g; + ort[m] = scale*ort[m]; + matrixH[m, m - 1] = scale*g; } } @@ -514,7 +515,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, i == j ? 1.0 : 0.0); + eigenVectors.At(i, j, i == j ? 1.0 : 0.0); } } @@ -532,14 +533,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var g = 0.0; for (var i = m; i < order; i++) { - g += ort[i] * EigenVectors.At(i, j); + g += ort[i]*eigenVectors.At(i, j); } // Double division avoids possible underflow - g = (g / ort[m]) / matrixH[m, m - 1]; + g = (g/ort[m])/matrixH[m, m - 1]; for (var i = m; i < order; i++) { - EigenVectors.At(i, j, EigenVectors.At(i, j) + g * ort[i]); + eigenVectors.At(i, j, eigenVectors.At(i, j) + g*ort[i]); } } } @@ -557,7 +558,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void NonsymmetricReduceHessenberToRealSchur(double[,] matrixH, double[] d, double[] e, int order) + static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVectors, double[,] matrixH, double[] d, double[] e, int order) { // Initialize var n = order - 1; @@ -590,7 +591,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization s = norm; } - if (Math.Abs(matrixH[l, l - 1]) < eps * s) + if (Math.Abs(matrixH[l, l - 1]) < eps*s) { break; } @@ -612,9 +613,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } else if (l == n - 1) { - w = matrixH[n, n - 1] * matrixH[n - 1, n]; - p = (matrixH[n - 1, n - 1] - matrixH[n, n]) / 2.0; - q = (p * p) + w; + w = matrixH[n, n - 1]*matrixH[n - 1, n]; + p = (matrixH[n - 1, n - 1] - matrixH[n, n])/2.0; + q = (p*p) + w; z = Math.Sqrt(Math.Abs(q)); matrixH[n, n] = matrixH[n, n] + exshift; matrixH[n - 1, n - 1] = matrixH[n - 1, n - 1] + exshift; @@ -637,41 +638,41 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization d[n] = d[n - 1]; if (z != 0.0) { - d[n] = x - (w / z); + d[n] = x - (w/z); } e[n - 1] = 0.0; e[n] = 0.0; x = matrixH[n, n - 1]; s = Math.Abs(x) + Math.Abs(z); - p = x / s; - q = z / s; - r = Math.Sqrt((p * p) + (q * q)); - p = p / r; - q = q / r; + p = x/s; + q = z/s; + r = Math.Sqrt((p*p) + (q*q)); + p = p/r; + q = q/r; // Row modification for (var j = n - 1; j < order; j++) { z = matrixH[n - 1, j]; - matrixH[n - 1, j] = (q * z) + (p * matrixH[n, j]); - matrixH[n, j] = (q * matrixH[n, j]) - (p * z); + matrixH[n - 1, j] = (q*z) + (p*matrixH[n, j]); + matrixH[n, j] = (q*matrixH[n, j]) - (p*z); } // Column modification for (var i = 0; i <= n; i++) { z = matrixH[i, n - 1]; - matrixH[i, n - 1] = (q * z) + (p * matrixH[i, n]); - matrixH[i, n] = (q * matrixH[i, n]) - (p * z); + matrixH[i, n - 1] = (q*z) + (p*matrixH[i, n]); + matrixH[i, n] = (q*matrixH[i, n]) - (p*z); } // Accumulate transformations for (var i = 0; i < order; i++) { - z = EigenVectors.At(i, n - 1); - EigenVectors.At(i, n - 1, (q * z) + (p * EigenVectors.At(i, n))); - EigenVectors.At(i, n, (q * EigenVectors.At(i, n)) - (p * z)); + z = eigenVectors.At(i, n - 1); + eigenVectors.At(i, n - 1, (q*z) + (p*eigenVectors.At(i, n))); + eigenVectors.At(i, n, (q*eigenVectors.At(i, n)) - (p*z)); } // Complex pair @@ -698,7 +699,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization if (l < n) { y = matrixH[n - 1, n - 1]; - w = matrixH[n, n - 1] * matrixH[n - 1, n]; + w = matrixH[n, n - 1]*matrixH[n - 1, n]; } // Wilkinson's original ad hoc shift @@ -711,15 +712,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } s = Math.Abs(matrixH[n, n - 1]) + Math.Abs(matrixH[n - 1, n - 2]); - x = y = 0.75 * s; - w = (-0.4375) * s * s; + x = y = 0.75*s; + w = (-0.4375)*s*s; } // MATLAB's new ad hoc shift if (iter == 30) { - s = (y - x) / 2.0; - s = (s * s) + w; + s = (y - x)/2.0; + s = (s*s) + w; if (s > 0) { s = Math.Sqrt(s); @@ -728,7 +729,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization s = -s; } - s = x - (w / (((y - x) / 2.0) + s)); + s = x - (w/(((y - x)/2.0) + s)); for (var i = 0; i <= n; i++) { matrixH[i, i] -= s; @@ -748,20 +749,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization z = matrixH[m, m]; r = x - z; s = y - z; - p = (((r * s) - w) / matrixH[m + 1, m]) + matrixH[m, m + 1]; + p = (((r*s) - w)/matrixH[m + 1, m]) + matrixH[m, m + 1]; q = matrixH[m + 1, m + 1] - z - r - s; r = matrixH[m + 2, m + 1]; s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - p = p / s; - q = q / s; - r = r / s; + p = p/s; + q = q/s; + r = r/s; if (m == l) { break; } - if (Math.Abs(matrixH[m, m - 1]) * (Math.Abs(q) + Math.Abs(r)) < eps * (Math.Abs(p) * (Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1])))) + if (Math.Abs(matrixH[m, m - 1])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1])))) { break; } @@ -791,9 +792,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); if (x != 0.0) { - p = p / x; - q = q / x; - r = r / x; + p = p/x; + q = q/x; + r = r/x; } } @@ -802,7 +803,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization break; } - s = Math.Sqrt((p * p) + (q * q) + (r * r)); + s = Math.Sqrt((p*p) + (q*q) + (r*r)); if (p < 0) { s = -s; @@ -812,7 +813,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { if (k != m) { - matrixH[k, k - 1] = (-s) * x; + matrixH[k, k - 1] = (-s)*x; } else if (l != m) { @@ -820,55 +821,55 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } p = p + s; - x = p / s; - y = q / s; - z = r / s; - q = q / p; - r = r / p; + x = p/s; + y = q/s; + z = r/s; + q = q/p; + r = r/p; // Row modification for (var j = k; j < order; j++) { - p = matrixH[k, j] + (q * matrixH[k + 1, j]); + p = matrixH[k, j] + (q*matrixH[k + 1, j]); if (notlast) { - p = p + (r * matrixH[k + 2, j]); - matrixH[k + 2, j] = matrixH[k + 2, j] - (p * z); + p = p + (r*matrixH[k + 2, j]); + matrixH[k + 2, j] = matrixH[k + 2, j] - (p*z); } - matrixH[k, j] = matrixH[k, j] - (p * x); - matrixH[k + 1, j] = matrixH[k + 1, j] - (p * y); + matrixH[k, j] = matrixH[k, j] - (p*x); + matrixH[k + 1, j] = matrixH[k + 1, j] - (p*y); } // Column modification for (var i = 0; i <= Math.Min(n, k + 3); i++) { - p = (x * matrixH[i, k]) + (y * matrixH[i, k + 1]); + p = (x*matrixH[i, k]) + (y*matrixH[i, k + 1]); if (notlast) { - p = p + (z * matrixH[i, k + 2]); - matrixH[i, k + 2] = matrixH[i, k + 2] - (p * r); + p = p + (z*matrixH[i, k + 2]); + matrixH[i, k + 2] = matrixH[i, k + 2] - (p*r); } matrixH[i, k] = matrixH[i, k] - p; - matrixH[i, k + 1] = matrixH[i, k + 1] - (p * q); + matrixH[i, k + 1] = matrixH[i, k + 1] - (p*q); } // Accumulate transformations for (var i = 0; i < order; i++) { - p = (x * EigenVectors.At(i, k)) + (y * EigenVectors.At(i, k + 1)); + p = (x*eigenVectors.At(i, k)) + (y*eigenVectors.At(i, k + 1)); if (notlast) { - p = p + (z * EigenVectors.At(i, k + 2)); - EigenVectors.At(i, k + 2, EigenVectors.At(i, k + 2) - (p * r)); + p = p + (z*eigenVectors.At(i, k + 2)); + eigenVectors.At(i, k + 2, eigenVectors.At(i, k + 2) - (p*r)); } - EigenVectors.At(i, k, EigenVectors.At(i, k) - p); - EigenVectors.At(i, k + 1, EigenVectors.At(i, k + 1) - (p * q)); + eigenVectors.At(i, k, eigenVectors.At(i, k) - p); + eigenVectors.At(i, k + 1, eigenVectors.At(i, k + 1) - (p*q)); } } // (s != 0) } // k loop @@ -899,7 +900,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization r = 0.0; for (var j = l; j <= n; j++) { - r = r + (matrixH[i, j] * matrixH[j, n]); + r = r + (matrixH[i, j]*matrixH[j, n]); } if (e[i] < 0.0) @@ -914,11 +915,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { if (w != 0.0) { - matrixH[i, n] = (-r) / w; + matrixH[i, n] = (-r)/w; } else { - matrixH[i, n] = (-r) / (eps * norm); + matrixH[i, n] = (-r)/(eps*norm); } // Solve real equations @@ -927,26 +928,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { x = matrixH[i, i + 1]; y = matrixH[i + 1, i]; - q = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]); - t = ((x * s) - (z * r)) / q; + q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); + t = ((x*s) - (z*r))/q; matrixH[i, n] = t; if (Math.Abs(x) > Math.Abs(z)) { - matrixH[i + 1, n] = (-r - (w * t)) / x; + matrixH[i + 1, n] = (-r - (w*t))/x; } else { - matrixH[i + 1, n] = (-s - (y * t)) / z; + matrixH[i + 1, n] = (-s - (y*t))/z; } } // Overflow control t = Math.Abs(matrixH[i, n]); - if ((eps * t) * t > 1) + if ((eps*t)*t > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n] = matrixH[j, n] / t; + matrixH[j, n] = matrixH[j, n]/t; } } } @@ -961,8 +962,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Last vector component imaginary so matrix is triangular if (Math.Abs(matrixH[n, n - 1]) > Math.Abs(matrixH[n - 1, n])) { - matrixH[n - 1, n - 1] = q / matrixH[n, n - 1]; - matrixH[n - 1, n] = (-(matrixH[n, n] - p)) / matrixH[n, n - 1]; + matrixH[n - 1, n - 1] = q/matrixH[n, n - 1]; + matrixH[n - 1, n] = (-(matrixH[n, n] - p))/matrixH[n, n - 1]; } else { @@ -979,8 +980,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization double sa = 0.0; for (var j = l; j <= n; j++) { - ra = ra + (matrixH[i, j] * matrixH[j, n - 1]); - sa = sa + (matrixH[i, j] * matrixH[j, n]); + ra = ra + (matrixH[i, j]*matrixH[j, n - 1]); + sa = sa + (matrixH[i, j]*matrixH[j, n]); } w = matrixH[i, i] - p; @@ -1006,24 +1007,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization x = matrixH[i, i + 1]; y = matrixH[i + 1, i]; - double vr = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]) - (q * q); - double vi = (d[i] - p) * 2.0 * q; + double vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); + double vi = (d[i] - p)*2.0*q; if ((vr == 0.0) && (vi == 0.0)) { - vr = eps * norm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); + vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); } - var res = Cdiv((x * r) - (z * ra) + (q * sa), (x * s) - (z * sa) - (q * ra), vr, vi); + var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); matrixH[i, n - 1] = res.Real; matrixH[i, n] = res.Imaginary; if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) { - matrixH[i + 1, n - 1] = (-ra - (w * matrixH[i, n - 1]) + (q * matrixH[i, n])) / x; - matrixH[i + 1, n] = (-sa - (w * matrixH[i, n]) - (q * matrixH[i, n - 1])) / x; + matrixH[i + 1, n - 1] = (-ra - (w*matrixH[i, n - 1]) + (q*matrixH[i, n]))/x; + matrixH[i + 1, n] = (-sa - (w*matrixH[i, n]) - (q*matrixH[i, n - 1]))/x; } else { - res = Cdiv(-r - (y * matrixH[i, n - 1]), -s - (y * matrixH[i, n]), z, q); + res = Cdiv(-r - (y*matrixH[i, n - 1]), -s - (y*matrixH[i, n]), z, q); matrixH[i + 1, n - 1] = res.Real; matrixH[i + 1, n] = res.Imaginary; } @@ -1031,12 +1032,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Overflow control t = Math.Max(Math.Abs(matrixH[i, n - 1]), Math.Abs(matrixH[i, n])); - if ((eps * t) * t > 1) + if ((eps*t)*t > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n - 1] = matrixH[j, n - 1] / t; - matrixH[j, n] = matrixH[j, n] / t; + matrixH[j, n - 1] = matrixH[j, n - 1]/t; + matrixH[j, n] = matrixH[j, n]/t; } } } @@ -1052,10 +1053,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization z = 0.0; for (var k = 0; k <= j; k++) { - z = z + (EigenVectors.At(i, k) * matrixH[k, j]); + z = z + (eigenVectors.At(i, k)*matrixH[k, j]); } - EigenVectors.At(i, j, z); + eigenVectors.At(i, j, z); } } } @@ -1068,14 +1069,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Real part of Y /// Imaginary part of Y /// Division result as a number. - private static Complex Cdiv(double xreal, double ximag, double yreal, double yimag) + static Complex Cdiv(double xreal, double ximag, double yreal, double yimag) { if (Math.Abs(yimag) < Math.Abs(yreal)) { - return new Complex((xreal + (ximag * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))), (ximag - (xreal * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal)))); + return new Complex((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); } - return new Complex((ximag + (xreal * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))), (-xreal + (ximag * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag)))); + return new Complex((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); } /// @@ -1085,17 +1086,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -1128,7 +1118,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j) * input.At(i, k); + value += EigenVectors.At(i, j)*input.At(i, k); } value /= EigenValues[j].Real; @@ -1142,7 +1132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization double value = 0; for (var i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result.At(j, k, value); @@ -1162,16 +1152,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) @@ -1199,7 +1179,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j) * input[i]; + value += EigenVectors.At(i, j)*input[i]; } value /= EigenValues[j].Real; @@ -1213,7 +1193,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization value = 0; for (int i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result[j] = value; @@ -1225,4 +1205,4 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } } } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 6d40080c..e4b5d074 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -486,7 +486,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Evd Evd() { - return new UserEvd(this); + return UserEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Factorization/Evd.cs index 3b810644..28785ef5 100644 --- a/src/Numerics/LinearAlgebra/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Factorization/Evd.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -57,10 +57,18 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization public abstract class Evd : ISolver where T : struct, IEquatable, IFormattable { + protected Evd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + { + EigenVectors = eigenVectors; + EigenValues = eigenValues; + D = blockDiagonal; + IsSymmetric = isSymmetric; + } + /// /// Gets or sets a value indicating whether matrix is symmetric or not /// - public bool IsSymmetric { get; protected set; } + public bool IsSymmetric { get; private set; } /// /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. @@ -82,17 +90,17 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// /// Gets or sets the eigen values (λ) of matrix in ascending value. /// - public Vector EigenValues { get; protected set; } + public Vector EigenValues { get; private set; } /// /// Gets or sets eigenvectors. /// - public Matrix EigenVectors { get; protected set; } + public Matrix EigenVectors { get; private set; } /// /// Gets or sets the block diagonal eigenvalue matrix. /// - public Matrix D { get; protected set; } + public Matrix D { get; private set; } /// /// Solves a system of linear equations, AX = B, with A SVD factorized. @@ -101,12 +109,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// The left hand side , X. public virtual Matrix Solve(Matrix input) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - var result = EigenVectors.CreateMatrix(EigenVectors.ColumnCount, input.ColumnCount); Solve(input, result); return result; @@ -126,12 +128,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// The left hand side , x. public virtual Vector Solve(Vector input) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - var x = EigenVectors.CreateVector(EigenVectors.ColumnCount); Solve(input, x); return x; diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index b8555c65..99b1c773 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -1062,7 +1062,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Evd Evd() { - return new DenseEvd(this); + return DenseEvd.Create(this); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs index fb76e240..1b281ef4 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { @@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public sealed class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,13 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public DenseEvd(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -79,22 +74,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var order = matrix.RowCount; // Initialize matrices for eigenvalues and eigenvectors - EigenVectors = matrix.CreateMatrix(order, order); - D = matrix.CreateMatrix(order, order); - EigenValues = new LinearAlgebra.Complex.DenseVector(order); + var eigenVectors = new DenseMatrix(order); + var blockDiagonal = new DenseMatrix(order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - IsSymmetric = true; + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i); } } - Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) EigenVectors).Values, - ((LinearAlgebra.Complex.DenseVector) EigenValues).Values, ((DenseMatrix) D).Values); + Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); + + return new DenseEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + DenseEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -1114,17 +1115,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -1191,16 +1181,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs index bb32a123..28b45dd5 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -29,6 +33,7 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; + #if NOSYSNUMERICS using Complex = Numerics.Complex; #else @@ -52,6 +57,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// public abstract class Evd : Evd { + protected Evd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { + } + /// /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. /// @@ -64,7 +74,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { det *= EigenValues[i]; - if (((Numerics.Complex32)EigenValues[i]).AlmostEqual(Numerics.Complex32.Zero)) + if (((Numerics.Complex32) EigenValues[i]).AlmostEqual(Numerics.Complex32.Zero)) { return 0; } @@ -85,7 +95,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var rank = 0; for (var i = 0; i < EigenValues.Count; i++) { - if (((Numerics.Complex32)EigenValues[i]).AlmostEqual(Numerics.Complex32.Zero)) + if (((Numerics.Complex32) EigenValues[i]).AlmostEqual(Numerics.Complex32.Zero)) { continue; } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs index 1f00bc81..3f7c3c15 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { @@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// conditioned, or even singular, so the validity of the equation /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public sealed class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the @@ -63,13 +63,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The matrix to factor. /// If is null. /// If EVD algorithm failed to converge with matrix . - public UserEvd(Matrix matrix) + public static UserEvd Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); @@ -78,57 +73,64 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var order = matrix.RowCount; // Initialize matricies for eigenvalues and eigenvectors - EigenVectors = matrix.CreateMatrix(order, order); - D = matrix.CreateMatrix(order, order); - EigenValues = new LinearAlgebra.Complex.DenseVector(order); - - IsSymmetric = true; + var eigenVectors = matrix.CreateMatrix(order, order); + var blockDiagonal = matrix.CreateMatrix(order, order); + var eigenValues = new LinearAlgebra.Complex.DenseVector(order); + + var isSymmetric = true; - for (var i = 0; IsSymmetric && i < order; i++) + for (var i = 0; isSymmetric && i < order; i++) { - for (var j = 0; IsSymmetric && j < order; j++) + for (var j = 0; isSymmetric && j < order; j++) { - IsSymmetric &= matrix.At(i, j) == matrix.At(j, i); + isSymmetric &= matrix.At(i, j) == matrix.At(j, i); } } var d = new float[order]; var e = new float[order]; - if (IsSymmetric) + if (isSymmetric) { - matrix.CopyTo(EigenVectors); - d = EigenVectors.Row(order - 1).ToArray(); + matrix.CopyTo(eigenVectors); + d = eigenVectors.Row(order - 1).ToArray(); - SymmetricTridiagonalize(d, e, order); - SymmetricDiagonalize(d, e, order); + SymmetricTridiagonalize(eigenVectors, d, e, order); + SymmetricDiagonalize(eigenVectors, d, e, order); } else { var matrixH = matrix.ToArray(); - NonsymmetricReduceToHessenberg(matrixH, order); - NonsymmetricReduceHessenberToRealSchur(matrixH, d, e, order); + NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(eigenVectors, matrixH, d, e, order); } for (var i = 0; i < order; i++) { - D.At(i, i, d[i]); + blockDiagonal.At(i, i, d[i]); if (e[i] > 0) { - D.At(i, i + 1, e[i]); + blockDiagonal.At(i, i + 1, e[i]); } else if (e[i] < 0) { - D.At(i, i - 1, e[i]); + blockDiagonal.At(i, i - 1, e[i]); } } for (var i = 0; i < order; i++) { - EigenValues[i] = new Complex(d[i], e[i]); + eigenValues[i] = new Complex(d[i], e[i]); } + + return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric); + } + + UserEvd(Matrix eigenVectors, Vector eigenValues, Matrix blockDiagonal, bool isSymmetric) + : base(eigenVectors, eigenValues, blockDiagonal, isSymmetric) + { } /// @@ -141,7 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void SymmetricTridiagonalize(float[] d, float[] e, int order) + static void SymmetricTridiagonalize(Matrix eigenVectors, float[] d, float[] e, int order) { // Householder reduction to tridiagonal form. for (var i = order - 1; i > 0; i--) @@ -160,9 +162,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization e[i] = d[i - 1]; for (var j = 0; j < i; j++) { - d[j] = EigenVectors.At(i - 1, j); - EigenVectors.At(i, j, 0.0f); - EigenVectors.At(j, i, 0.0f); + d[j] = eigenVectors.At(i - 1, j); + eigenVectors.At(i, j, 0.0f); + eigenVectors.At(j, i, 0.0f); } } else @@ -171,18 +173,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var k = 0; k < i; k++) { d[k] /= scale; - h += d[k] * d[k]; + h += d[k]*d[k]; } var f = d[i - 1]; - var g = (float)Math.Sqrt(h); + var g = (float) Math.Sqrt(h); if (f > 0) { g = -g; } - e[i] = scale * g; - h = h - (f * g); + e[i] = scale*g; + h = h - (f*g); d[i - 1] = f - g; for (var j = 0; j < i; j++) @@ -194,13 +196,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var j = 0; j < i; j++) { f = d[j]; - EigenVectors.At(j, i, f); - g = e[j] + (EigenVectors.At(j, j) * f); + eigenVectors.At(j, i, f); + g = e[j] + (eigenVectors.At(j, j)*f); for (var k = j + 1; k <= i - 1; k++) { - g += EigenVectors.At(k, j) * d[k]; - e[k] += EigenVectors.At(k, j) * f; + g += eigenVectors.At(k, j)*d[k]; + e[k] += eigenVectors.At(k, j)*f; } e[j] = g; @@ -211,14 +213,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var j = 0; j < i; j++) { e[j] /= h; - f += e[j] * d[j]; + f += e[j]*d[j]; } - var hh = f / (h + h); + var hh = f/(h + h); for (var j = 0; j < i; j++) { - e[j] -= hh * d[j]; + e[j] -= hh*d[j]; } for (var j = 0; j < i; j++) @@ -228,11 +230,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var k = j; k <= i - 1; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - (f * e[k]) - (g * d[k])); + eigenVectors.At(k, j, eigenVectors.At(k, j) - (f*e[k]) - (g*d[k])); } - d[j] = EigenVectors.At(i - 1, j); - EigenVectors.At(i, j, 0.0f); + d[j] = eigenVectors.At(i - 1, j); + eigenVectors.At(i, j, 0.0f); } } @@ -242,14 +244,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Accumulate transformations. for (var i = 0; i < order - 1; i++) { - EigenVectors.At(order - 1, i, EigenVectors.At(i, i)); - EigenVectors.At(i, i, 1.0f); + eigenVectors.At(order - 1, i, eigenVectors.At(i, i)); + eigenVectors.At(i, i, 1.0f); var h = d[i + 1]; if (h != 0.0f) { for (var k = 0; k <= i; k++) { - d[k] = EigenVectors.At(k, i + 1) / h; + d[k] = eigenVectors.At(k, i + 1)/h; } for (var j = 0; j <= i; j++) @@ -257,29 +259,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var g = 0.0f; for (var k = 0; k <= i; k++) { - g += EigenVectors.At(k, i + 1) * EigenVectors.At(k, j); + g += eigenVectors.At(k, i + 1)*eigenVectors.At(k, j); } for (var k = 0; k <= i; k++) { - EigenVectors.At(k, j, EigenVectors.At(k, j) - g * d[k]); + eigenVectors.At(k, j, eigenVectors.At(k, j) - g*d[k]); } } } for (var k = 0; k <= i; k++) { - EigenVectors.At(k, i + 1, 0.0f); + eigenVectors.At(k, i + 1, 0.0f); } } for (var j = 0; j < order; j++) { - d[j] = EigenVectors.At(order - 1, j); - EigenVectors.At(order - 1, j, 0.0f); + d[j] = eigenVectors.At(order - 1, j); + eigenVectors.At(order - 1, j, 0.0f); } - EigenVectors.At(order - 1, order - 1, 1.0f); + eigenVectors.At(order - 1, order - 1, 1.0f); e[0] = 0.0f; } @@ -294,7 +296,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. /// - private void SymmetricDiagonalize(float[] d, float[] e, int order) + static void SymmetricDiagonalize(Matrix eigenVectors, float[] d, float[] e, int order) { const int maxiter = 1000; @@ -315,7 +317,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var m = l; while (m < order) { - if (Math.Abs(e[m]) <= eps * tst1) + if (Math.Abs(e[m]) <= eps*tst1) { break; } @@ -334,15 +336,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Compute implicit shift var g = d[l]; - var p = (d[l + 1] - g) / (2.0f * e[l]); + var p = (d[l + 1] - g)/(2.0f*e[l]); var r = SpecialFunctions.Hypotenuse(p, 1.0f); if (p < 0) { r = -r; } - d[l] = e[l] / (p + r); - d[l + 1] = e[l] * (p + r); + d[l] = e[l]/(p + r); + d[l + 1] = e[l]*(p + r); var dl1 = d[l + 1]; var h = g - d[l]; @@ -366,27 +368,27 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization c3 = c2; c2 = c; s2 = s; - g = c * e[i]; - h = c * p; + g = c*e[i]; + h = c*p; r = SpecialFunctions.Hypotenuse(p, e[i]); - e[i + 1] = s * r; - s = e[i] / r; - c = p / r; - p = (c * d[i]) - (s * g); - d[i + 1] = h + (s * ((c * g) + (s * d[i]))); + e[i + 1] = s*r; + s = e[i]/r; + c = p/r; + p = (c*d[i]) - (s*g); + d[i + 1] = h + (s*((c*g) + (s*d[i]))); // Accumulate transformation. for (var k = 0; k < order; k++) { - h = EigenVectors.At(k, i + 1); - EigenVectors.At(k, i + 1, (s * EigenVectors.At(k, i)) + (c * h)); - EigenVectors.At(k, i, (c * EigenVectors.At(k, i)) - (s * h)); + h = eigenVectors.At(k, i + 1); + eigenVectors.At(k, i + 1, (s*eigenVectors.At(k, i)) + (c*h)); + eigenVectors.At(k, i, (c*eigenVectors.At(k, i)) - (s*h)); } } - p = (-s) * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; + p = (-s)*s2*c3*el1*e[l]/dl1; + e[l] = s*p; + d[l] = c*p; // Check for convergence. If too many iterations have been performed, // throw exception that Convergence Failed @@ -394,8 +396,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { throw new NonConvergenceException(); } - } - while (Math.Abs(e[l]) > eps * tst1); + } while (Math.Abs(e[l]) > eps*tst1); } d[l] = d[l] + f; @@ -422,9 +423,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization d[i] = p; for (var j = 0; j < order; j++) { - p = EigenVectors.At(j, i); - EigenVectors.At(j, i, EigenVectors.At(j, k)); - EigenVectors.At(j, k, p); + p = eigenVectors.At(j, i); + eigenVectors.At(j, i, eigenVectors.At(j, k)); + eigenVectors.At(j, k, p); } } } @@ -439,7 +440,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutines in EISPACK. - private void NonsymmetricReduceToHessenberg(float[,] matrixH, int order) + static void NonsymmetricReduceToHessenberg(Matrix eigenVectors, float[,] matrixH, int order) { var ort = new float[order]; @@ -458,17 +459,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var h = 0.0f; for (var i = order - 1; i >= m; i--) { - ort[i] = matrixH[i, m - 1] / scale; - h += ort[i] * ort[i]; + ort[i] = matrixH[i, m - 1]/scale; + h += ort[i]*ort[i]; } - var g = (float)Math.Sqrt(h); + var g = (float) Math.Sqrt(h); if (ort[m] > 0) { g = -g; } - h = h - (ort[m] * g); + h = h - (ort[m]*g); ort[m] = ort[m] - g; // Apply Householder similarity transformation @@ -478,13 +479,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var f = 0.0f; for (var i = order - 1; i >= m; i--) { - f += ort[i] * matrixH[i, j]; + f += ort[i]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var i = m; i < order; i++) { - matrixH[i, j] -= f * ort[i]; + matrixH[i, j] -= f*ort[i]; } } @@ -493,18 +494,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var f = 0.0f; for (var j = order - 1; j >= m; j--) { - f += ort[j] * matrixH[i, j]; + f += ort[j]*matrixH[i, j]; } - f = f / h; + f = f/h; for (var j = m; j < order; j++) { - matrixH[i, j] -= f * ort[j]; + matrixH[i, j] -= f*ort[j]; } } - ort[m] = scale * ort[m]; - matrixH[m, m - 1] = scale * g; + ort[m] = scale*ort[m]; + matrixH[m, m - 1] = scale*g; } } @@ -513,7 +514,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var j = 0; j < order; j++) { - EigenVectors.At(i, j, i == j ? 1.0f : 0.0f); + eigenVectors.At(i, j, i == j ? 1.0f : 0.0f); } } @@ -531,14 +532,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var g = 0.0f; for (var i = m; i < order; i++) { - g += ort[i] * EigenVectors.At(i, j); + g += ort[i]*eigenVectors.At(i, j); } // Double division avoids possible underflow - g = (g / ort[m]) / matrixH[m, m - 1]; + g = (g/ort[m])/matrixH[m, m - 1]; for (var i = m; i < order; i++) { - EigenVectors.At(i, j, EigenVectors.At(i, j) + g * ort[i]); + eigenVectors.At(i, j, eigenVectors.At(i, j) + g*ort[i]); } } } @@ -556,11 +557,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// by Martin and Wilkinson, Handbook for Auto. Comp., /// Vol.ii-Linear Algebra, and the corresponding /// Fortran subroutine in EISPACK. - private void NonsymmetricReduceHessenberToRealSchur(float[,] matrixH, float[] d, float[] e, int order) + static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVectors, float[,] matrixH, float[] d, float[] e, int order) { // Initialize var n = order - 1; - var eps = (float)Precision.SingleMachinePrecision; + var eps = (float) Precision.SingleMachinePrecision; var exshift = 0.0f; float p = 0, q = 0, r = 0, s = 0, z = 0, w, x, y; @@ -589,7 +590,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization s = norm; } - if (Math.Abs(matrixH[l, l - 1]) < eps * s) + if (Math.Abs(matrixH[l, l - 1]) < eps*s) { break; } @@ -611,10 +612,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } else if (l == n - 1) { - w = matrixH[n, n - 1] * matrixH[n - 1, n]; - p = (matrixH[n - 1, n - 1] - matrixH[n, n]) / 2.0f; - q = (p * p) + w; - z = (float)Math.Sqrt(Math.Abs(q)); + w = matrixH[n, n - 1]*matrixH[n - 1, n]; + p = (matrixH[n - 1, n - 1] - matrixH[n, n])/2.0f; + q = (p*p) + w; + z = (float) Math.Sqrt(Math.Abs(q)); matrixH[n, n] = matrixH[n, n] + exshift; matrixH[n - 1, n - 1] = matrixH[n - 1, n - 1] + exshift; x = matrixH[n, n]; @@ -636,41 +637,41 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization d[n] = d[n - 1]; if (z != 0.0f) { - d[n] = x - (w / z); + d[n] = x - (w/z); } e[n - 1] = 0.0f; e[n] = 0.0f; x = matrixH[n, n - 1]; s = Math.Abs(x) + Math.Abs(z); - p = x / s; - q = z / s; - r = (float)Math.Sqrt((p * p) + (q * q)); - p = p / r; - q = q / r; + p = x/s; + q = z/s; + r = (float) Math.Sqrt((p*p) + (q*q)); + p = p/r; + q = q/r; // Row modification for (var j = n - 1; j < order; j++) { z = matrixH[n - 1, j]; - matrixH[n - 1, j] = (q * z) + (p * matrixH[n, j]); - matrixH[n, j] = (q * matrixH[n, j]) - (p * z); + matrixH[n - 1, j] = (q*z) + (p*matrixH[n, j]); + matrixH[n, j] = (q*matrixH[n, j]) - (p*z); } // Column modification for (var i = 0; i <= n; i++) { z = matrixH[i, n - 1]; - matrixH[i, n - 1] = (q * z) + (p * matrixH[i, n]); - matrixH[i, n] = (q * matrixH[i, n]) - (p * z); + matrixH[i, n - 1] = (q*z) + (p*matrixH[i, n]); + matrixH[i, n] = (q*matrixH[i, n]) - (p*z); } // Accumulate transformations for (var i = 0; i < order; i++) { - z = EigenVectors.At(i, n - 1); - EigenVectors.At(i, n - 1, (q * z) + (p * EigenVectors.At(i, n))); - EigenVectors.At(i, n, (q * EigenVectors.At(i, n)) - (p * z)); + z = eigenVectors.At(i, n - 1); + eigenVectors.At(i, n - 1, (q*z) + (p*eigenVectors.At(i, n))); + eigenVectors.At(i, n, (q*eigenVectors.At(i, n)) - (p*z)); } // Complex pair @@ -697,7 +698,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization if (l < n) { y = matrixH[n - 1, n - 1]; - w = matrixH[n, n - 1] * matrixH[n - 1, n]; + w = matrixH[n, n - 1]*matrixH[n - 1, n]; } // Wilkinson's original ad hoc shift @@ -710,24 +711,24 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } s = Math.Abs(matrixH[n, n - 1]) + Math.Abs(matrixH[n - 1, n - 2]); - x = y = 0.75f * s; - w = (-0.4375f) * s * s; + x = y = 0.75f*s; + w = (-0.4375f)*s*s; } // MATLAB's new ad hoc shift if (iter == 30) { - s = (y - x) / 2.0f; - s = (s * s) + w; + s = (y - x)/2.0f; + s = (s*s) + w; if (s > 0) { - s = (float)Math.Sqrt(s); + s = (float) Math.Sqrt(s); if (y < x) { s = -s; } - s = x - (w / (((y - x) / 2.0f) + s)); + s = x - (w/(((y - x)/2.0f) + s)); for (var i = 0; i <= n; i++) { matrixH[i, i] -= s; @@ -747,20 +748,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization z = matrixH[m, m]; r = x - z; s = y - z; - p = (((r * s) - w) / matrixH[m + 1, m]) + matrixH[m, m + 1]; + p = (((r*s) - w)/matrixH[m + 1, m]) + matrixH[m, m + 1]; q = matrixH[m + 1, m + 1] - z - r - s; r = matrixH[m + 2, m + 1]; s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - p = p / s; - q = q / s; - r = r / s; + p = p/s; + q = q/s; + r = r/s; if (m == l) { break; } - if (Math.Abs(matrixH[m, m - 1]) * (Math.Abs(q) + Math.Abs(r)) < eps * (Math.Abs(p) * (Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1])))) + if (Math.Abs(matrixH[m, m - 1])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1])))) { break; } @@ -790,9 +791,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); if (x != 0.0f) { - p = p / x; - q = q / x; - r = r / x; + p = p/x; + q = q/x; + r = r/x; } } @@ -801,7 +802,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization break; } - s = (float)Math.Sqrt((p * p) + (q * q) + (r * r)); + s = (float) Math.Sqrt((p*p) + (q*q) + (r*r)); if (p < 0) { s = -s; @@ -811,7 +812,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { if (k != m) { - matrixH[k, k - 1] = (-s) * x; + matrixH[k, k - 1] = (-s)*x; } else if (l != m) { @@ -819,55 +820,55 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } p = p + s; - x = p / s; - y = q / s; - z = r / s; - q = q / p; - r = r / p; + x = p/s; + y = q/s; + z = r/s; + q = q/p; + r = r/p; // Row modification for (var j = k; j < order; j++) { - p = matrixH[k, j] + (q * matrixH[k + 1, j]); + p = matrixH[k, j] + (q*matrixH[k + 1, j]); if (notlast) { - p = p + (r * matrixH[k + 2, j]); - matrixH[k + 2, j] = matrixH[k + 2, j] - (p * z); + p = p + (r*matrixH[k + 2, j]); + matrixH[k + 2, j] = matrixH[k + 2, j] - (p*z); } - matrixH[k, j] = matrixH[k, j] - (p * x); - matrixH[k + 1, j] = matrixH[k + 1, j] - (p * y); + matrixH[k, j] = matrixH[k, j] - (p*x); + matrixH[k + 1, j] = matrixH[k + 1, j] - (p*y); } // Column modification for (var i = 0; i <= Math.Min(n, k + 3); i++) { - p = (x * matrixH[i, k]) + (y * matrixH[i, k + 1]); + p = (x*matrixH[i, k]) + (y*matrixH[i, k + 1]); if (notlast) { - p = p + (z * matrixH[i, k + 2]); - matrixH[i, k + 2] = matrixH[i, k + 2] - (p * r); + p = p + (z*matrixH[i, k + 2]); + matrixH[i, k + 2] = matrixH[i, k + 2] - (p*r); } matrixH[i, k] = matrixH[i, k] - p; - matrixH[i, k + 1] = matrixH[i, k + 1] - (p * q); + matrixH[i, k + 1] = matrixH[i, k + 1] - (p*q); } // Accumulate transformations for (var i = 0; i < order; i++) { - p = (x * EigenVectors.At(i, k)) + (y * EigenVectors.At(i, k + 1)); + p = (x*eigenVectors.At(i, k)) + (y*eigenVectors.At(i, k + 1)); if (notlast) { - p = p + (z * EigenVectors.At(i, k + 2)); - EigenVectors.At(i, k + 2, EigenVectors.At(i, k + 2) - (p * r)); + p = p + (z*eigenVectors.At(i, k + 2)); + eigenVectors.At(i, k + 2, eigenVectors.At(i, k + 2) - (p*r)); } - EigenVectors.At(i, k, EigenVectors.At(i, k) - p); - EigenVectors.At(i, k + 1, EigenVectors.At(i, k + 1) - (p * q)); + eigenVectors.At(i, k, eigenVectors.At(i, k) - p); + eigenVectors.At(i, k + 1, eigenVectors.At(i, k + 1) - (p*q)); } } // (s != 0) } // k loop @@ -898,7 +899,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization r = 0.0f; for (var j = l; j <= n; j++) { - r = r + (matrixH[i, j] * matrixH[j, n]); + r = r + (matrixH[i, j]*matrixH[j, n]); } if (e[i] < 0.0f) @@ -913,11 +914,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { if (w != 0.0f) { - matrixH[i, n] = (-r) / w; + matrixH[i, n] = (-r)/w; } else { - matrixH[i, n] = (-r) / (eps * norm); + matrixH[i, n] = (-r)/(eps*norm); } // Solve real equations @@ -926,26 +927,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { x = matrixH[i, i + 1]; y = matrixH[i + 1, i]; - q = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]); - t = ((x * s) - (z * r)) / q; + q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); + t = ((x*s) - (z*r))/q; matrixH[i, n] = t; if (Math.Abs(x) > Math.Abs(z)) { - matrixH[i + 1, n] = (-r - (w * t)) / x; + matrixH[i + 1, n] = (-r - (w*t))/x; } else { - matrixH[i + 1, n] = (-s - (y * t)) / z; + matrixH[i + 1, n] = (-s - (y*t))/z; } } // Overflow control t = Math.Abs(matrixH[i, n]); - if ((eps * t) * t > 1) + if ((eps*t)*t > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n] = matrixH[j, n] / t; + matrixH[j, n] = matrixH[j, n]/t; } } } @@ -960,8 +961,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Last vector component imaginary so matrix is triangular if (Math.Abs(matrixH[n, n - 1]) > Math.Abs(matrixH[n - 1, n])) { - matrixH[n - 1, n - 1] = q / matrixH[n, n - 1]; - matrixH[n - 1, n] = (-(matrixH[n, n] - p)) / matrixH[n, n - 1]; + matrixH[n - 1, n - 1] = q/matrixH[n, n - 1]; + matrixH[n - 1, n] = (-(matrixH[n, n] - p))/matrixH[n, n - 1]; } else { @@ -978,8 +979,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization float sa = 0.0f; for (var j = l; j <= n; j++) { - ra = ra + (matrixH[i, j] * matrixH[j, n - 1]); - sa = sa + (matrixH[i, j] * matrixH[j, n]); + ra = ra + (matrixH[i, j]*matrixH[j, n - 1]); + sa = sa + (matrixH[i, j]*matrixH[j, n]); } w = matrixH[i, i] - p; @@ -1005,24 +1006,24 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization x = matrixH[i, i + 1]; y = matrixH[i + 1, i]; - float vr = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]) - (q * q); - float vi = (d[i] - p) * 2.0f * q; + float vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); + float vi = (d[i] - p)*2.0f*q; if ((vr == 0.0f) && (vi == 0.0f)) { - vr = eps * norm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); + vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); } - var res = Cdiv((x * r) - (z * ra) + (q * sa), (x * s) - (z * sa) - (q * ra), vr, vi); + var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); matrixH[i, n - 1] = res.Real; matrixH[i, n] = res.Imaginary; if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) { - matrixH[i + 1, n - 1] = (-ra - (w * matrixH[i, n - 1]) + (q * matrixH[i, n])) / x; - matrixH[i + 1, n] = (-sa - (w * matrixH[i, n]) - (q * matrixH[i, n - 1])) / x; + matrixH[i + 1, n - 1] = (-ra - (w*matrixH[i, n - 1]) + (q*matrixH[i, n]))/x; + matrixH[i + 1, n] = (-sa - (w*matrixH[i, n]) - (q*matrixH[i, n - 1]))/x; } else { - res = Cdiv(-r - (y * matrixH[i, n - 1]), -s - (y * matrixH[i, n]), z, q); + res = Cdiv(-r - (y*matrixH[i, n - 1]), -s - (y*matrixH[i, n]), z, q); matrixH[i + 1, n - 1] = res.Real; matrixH[i + 1, n] = res.Imaginary; } @@ -1030,12 +1031,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Overflow control t = Math.Max(Math.Abs(matrixH[i, n - 1]), Math.Abs(matrixH[i, n])); - if ((eps * t) * t > 1) + if ((eps*t)*t > 1) { for (var j = i; j <= n; j++) { - matrixH[j, n - 1] = matrixH[j, n - 1] / t; - matrixH[j, n] = matrixH[j, n] / t; + matrixH[j, n - 1] = matrixH[j, n - 1]/t; + matrixH[j, n] = matrixH[j, n]/t; } } } @@ -1051,10 +1052,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization z = 0.0f; for (var k = 0; k <= j; k++) { - z = z + (EigenVectors.At(i, k) * matrixH[k, j]); + z = z + (eigenVectors.At(i, k)*matrixH[k, j]); } - EigenVectors.At(i, j, z); + eigenVectors.At(i, j, z); } } } @@ -1067,14 +1068,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Real part of Y /// Imaginary part of Y /// Division result as a number. - private static Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag) + static Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag) { if (Math.Abs(yimag) < Math.Abs(yreal)) { - return new Complex32((xreal + (ximag * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))), (ximag - (xreal * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal)))); + return new Complex32((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); } - return new Complex32((ximag + (xreal * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))), (-xreal + (ximag * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag)))); + return new Complex32((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); } /// @@ -1084,17 +1085,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // The solution X should have the same number of columns as B if (input.ColumnCount != result.ColumnCount) { @@ -1127,10 +1117,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j) * input.At(i, k); + value += EigenVectors.At(i, j)*input.At(i, k); } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; @@ -1141,7 +1131,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization float value = 0; for (var i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result.At(j, k, value); @@ -1161,16 +1151,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - // Ax=b where A is an m x m matrix // Check that b is a column vector with m entries if (EigenValues.Count != input.Count) @@ -1198,10 +1178,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var i = 0; i < order; i++) { - value += EigenVectors.At(i, j) * input[i]; + value += EigenVectors.At(i, j)*input[i]; } - value /= (float)EigenValues[j].Real; + value /= (float) EigenValues[j].Real; } tmp[j] = value; @@ -1212,7 +1192,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization value = 0; for (int i = 0; i < order; i++) { - value += EigenVectors.At(j, i) * tmp[i]; + value += EigenVectors.At(j, i)*tmp[i]; } result[j] = value; @@ -1224,4 +1204,4 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } } } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 911220c4..109d9419 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -486,7 +486,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Evd Evd() { - return new UserEvd(this); + return UserEvd.Create(this); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs index c7096454..8433fe07 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs @@ -24,28 +24,18 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Complex; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Complex; - using LinearAlgebra.Complex.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for a dense matrix. /// public class EvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNullThrowsArgumentNullException() - { - Assert.Throws(() => new DenseEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs index 083de445..0f51a7da 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs @@ -24,27 +24,17 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Complex.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for an user matrix. /// public class UserEvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs index f0cf8d00..232987bb 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs @@ -24,13 +24,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Complex32; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Complex32; - using LinearAlgebra.Complex32.Factorization; - using NUnit.Framework; using Complex32 = Numerics.Complex32; /// @@ -38,15 +37,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization /// public class EvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNullThrowsArgumentNullException() - { - Assert.Throws(() => new DenseEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs index 4a702093..37a36eba 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs @@ -24,12 +24,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Complex32.Factorization; - using NUnit.Framework; using Complex32 = Numerics.Complex32; /// @@ -37,15 +36,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization /// public class UserEvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs index c2b39503..98357978 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs @@ -24,28 +24,18 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Double; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Double; - using LinearAlgebra.Double.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for a dense matrix. /// public class EvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNullThrowsArgumentNullException() - { - Assert.Throws(() => new DenseEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs index 59993baa..73c8f39a 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs @@ -24,27 +24,17 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Double.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for an user matrix. /// public class UserEvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs index c45f8fbb..827d9b7f 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs @@ -24,28 +24,18 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Single; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Single; - using LinearAlgebra.Single.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for a dense matrix. /// public class EvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNullThrowsArgumentNullException() - { - Assert.Throws(() => new DenseEvd(null)); - } - /// /// Can factorize identity matrix. /// diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs index b284422e..efb2addb 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs @@ -24,27 +24,17 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Single.Factorization; - using NUnit.Framework; /// /// Eigenvalues factorization tests for an user matrix. /// public class UserEvdTests { - /// - /// Constructor null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserEvd(null)); - } - /// /// Can factorize identity matrix. ///