From f4ea07ce212ab7353731c50793294ae84bf2b8bf Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 8 Sep 2013 19:28:27 +0200 Subject: [PATCH] LA: Simplify SVD decomposition architecture --- .../LinearAlgebra/Factorization/Svd.cs | 7 +- .../LinearAlgebra/Complex/DenseMatrix.cs | 2 +- .../Complex/Factorization/DenseSvd.cs | 58 +--- .../Complex/Factorization/Svd.cs | 11 +- .../Complex/Factorization/UserSvd.cs | 323 ++++++++--------- src/Numerics/LinearAlgebra/Complex/Matrix.cs | 2 +- .../LinearAlgebra/Complex32/DenseMatrix.cs | 2 +- .../Complex32/Factorization/DenseSvd.cs | 58 +--- .../Complex32/Factorization/Svd.cs | 11 +- .../Complex32/Factorization/UserSvd.cs | 326 ++++++++---------- .../LinearAlgebra/Complex32/Matrix.cs | 2 +- .../LinearAlgebra/Double/DenseMatrix.cs | 2 +- .../Double/Factorization/DenseSvd.cs | 58 +--- .../LinearAlgebra/Double/Factorization/Svd.cs | 11 +- .../Double/Factorization/UserSvd.cs | 322 ++++++++--------- src/Numerics/LinearAlgebra/Double/Matrix.cs | 2 +- .../LinearAlgebra/Factorization/Svd.cs | 96 +++--- .../LinearAlgebra/Single/DenseMatrix.cs | 2 +- .../Single/Factorization/DenseSvd.cs | 58 +--- .../LinearAlgebra/Single/Factorization/Svd.cs | 11 +- .../Single/Factorization/UserSvd.cs | 326 ++++++++---------- src/Numerics/LinearAlgebra/Single/Matrix.cs | 2 +- .../Complex/Factorization/SvdTests.cs | 17 +- .../Complex/Factorization/UserSvdTests.cs | 20 +- .../Complex32/Factorization/SvdTests.cs | 17 +- .../Complex32/Factorization/UserSvdTests.cs | 22 +- .../Double/Factorization/SvdTests.cs | 17 +- .../Double/Factorization/UserSvdTests.cs | 21 +- .../Single/Factorization/SvdTests.cs | 17 +- .../Single/Factorization/UserSvdTests.cs | 21 +- 30 files changed, 821 insertions(+), 1023 deletions(-) diff --git a/src/Examples/LinearAlgebra/Factorization/Svd.cs b/src/Examples/LinearAlgebra/Factorization/Svd.cs index cfb950b3..dd8e65f7 100644 --- a/src/Examples/LinearAlgebra/Factorization/Svd.cs +++ b/src/Examples/LinearAlgebra/Factorization/Svd.cs @@ -26,6 +26,7 @@ using System; using System.Globalization; +using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; namespace Examples.LinearAlgebra.FactorizationExamples @@ -97,7 +98,7 @@ namespace Examples.LinearAlgebra.FactorizationExamples // 3. Singular values as diagonal matrix Console.WriteLine(@"3. Singular values as diagonal matrix"); - Console.WriteLine(svd.W().ToString("#0.00\t", formatProvider)); + Console.WriteLine(svd.W.ToString("#0.00\t", formatProvider)); Console.WriteLine(); // 4. Right singular vectors @@ -118,7 +119,7 @@ namespace Examples.LinearAlgebra.FactorizationExamples Console.WriteLine(); // 7. Reconstruct initial matrix: A = U*Σ*VT - var reconstruct = svd.U * svd.W() * svd.VT; + var reconstruct = svd.U * svd.W * svd.VT; Console.WriteLine(@"7. Reconstruct initial matrix: A = U*S*VT"); Console.WriteLine(reconstruct.ToString("#0.00\t", formatProvider)); Console.WriteLine(); @@ -154,7 +155,7 @@ namespace Examples.LinearAlgebra.FactorizationExamples // 13. Singular values as diagonal matrix Console.WriteLine(@"13. Singular values as diagonal matrix"); - Console.WriteLine(svd.W().ToString("#0.00\t", formatProvider)); + Console.WriteLine(svd.W.ToString("#0.00\t", formatProvider)); Console.WriteLine(); // 14. Access to left singular vectors when partial SVD decomposition was performed diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index ce1eea7c..8a687a38 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -1031,7 +1031,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Svd Svd(bool computeVectors) { - return new DenseSvd(this, computeVectors); + return DenseSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs index ed56d62a..60fb6f99 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.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 { @@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public sealed class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,19 +64,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// If SVD algorithm failed to converge with matrix . - public DenseSvd(DenseMatrix matrix, bool computeVectors) + public static DenseSvd Create(DenseMatrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); - S = new DenseVector(nm); - U = new DenseMatrix(matrix.RowCount); - VT = new DenseMatrix(matrix.ColumnCount); - Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values); + var s = new DenseVector(nm); + var u = new DenseMatrix(matrix.RowCount); + var vt = new DenseMatrix(matrix.ColumnCount); + Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix) matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, s.Values, u.Values, vt.Values); + + return new DenseSvd(s, u, vt, computeVectors); + } + + DenseSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -86,18 +87,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -132,7 +122,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, input.ColumnCount, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values); } /// @@ -142,17 +132,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -182,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, 1, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs index fcc2fd86..5597de74 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.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,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.LinearAlgebra.Factorization; -using MathNet.Numerics.Properties; using System; using System.Linq; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { @@ -58,6 +58,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// public abstract class Svd : Svd { + protected Svd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { + } + /// /// Gets the effective numerical matrix rank. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs index 9e6b1fb2..8b4b2002 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.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 { @@ -38,6 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using Numerics; #else using System.Numerics; + #endif /// @@ -54,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public sealed class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -64,20 +65,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// - public UserSvd(Matrix matrix, bool computeVectors) + public static UserSvd Create(Matrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); var matrixCopy = matrix.Clone(); - S = matrixCopy.CreateVector(nm); - U = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); - VT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); + var s = matrixCopy.CreateVector(nm); + var u = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); + var vt = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); const int maxiter = 1000; var e = new Complex[matrixCopy.ColumnCount]; @@ -100,34 +95,34 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization if (l < nct) { // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. - S[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); - if (S[l].Magnitude != 0.0) + s[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); + if (s[l].Magnitude != 0.0) { if (matrixCopy.At(l, l).Magnitude != 0.0) { - S[l] = Csign(S[l], matrixCopy.At(l, l)); + s[l] = Csign(s[l], matrixCopy.At(l, l)); } - CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / S[l]); + CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0/s[l]); matrixCopy.At(l, l, (Complex.One + matrixCopy.At(l, l))); } - S[l] = -S[l]; + s[l] = -s[l]; } for (j = lp1; j < matrixCopy.ColumnCount; j++) { if (l < nct) { - if (S[l].Magnitude != 0.0) + if (s[l].Magnitude != 0.0) { // Apply the transformation. - t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); + t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l)/matrixCopy.At(l, l); if (t != Complex.Zero) { for (var ii = l; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t*matrixCopy.At(ii, l))); } } } @@ -138,12 +133,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization e[j] = matrixCopy.At(l, j).Conjugate(); } - if (ComputeVectors && l < nct) + if (computeVectors && l < nct) { // Place the transformation in u for subsequent back multiplication. for (i = l; i < matrixCopy.RowCount; i++) { - U.At(i, l, matrixCopy.At(i, l)); + u.At(i, l, matrixCopy.At(i, l)); } } @@ -162,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization e[l] = Csign(e[l], e[lp1]); } - CscalVector(e, lp1, 1.0 / e[l]); + CscalVector(e, lp1, 1.0/e[l]); e[lp1] = Complex.One + e[lp1]; } @@ -181,30 +176,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - work[ii] += e[j] * matrixCopy.At(ii, j); + work[ii] += e[j]*matrixCopy.At(ii, j); } } } for (j = lp1; j < matrixCopy.ColumnCount; j++) { - var ww = (-e[j] / e[lp1]).Conjugate(); + var ww = (-e[j]/e[lp1]).Conjugate(); if (ww != Complex.Zero) { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww*work[ii])); } } } } - if (ComputeVectors) + if (computeVectors) { // Place the transformation in v for subsequent back multiplication. for (i = lp1; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, e[i]); + vt.At(i, l, e[i]); } } } @@ -215,12 +210,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var nrtp1 = nrt + 1; if (nct < matrixCopy.ColumnCount) { - S[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); + s[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); } if (matrixCopy.RowCount < m) { - S[m - 1] = Complex.Zero; + s[m - 1] = Complex.Zero; } if (nrtp1 < m) @@ -231,55 +226,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization e[m - 1] = Complex.Zero; // If required, generate u. - if (ComputeVectors) + if (computeVectors) { for (j = nctp1 - 1; j < ncu; j++) { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, j, Complex.Zero); + u.At(i, j, Complex.Zero); } - U.At(j, j, Complex.One); + u.At(j, j, Complex.One); } for (l = nct - 1; l >= 0; l--) { - if (S[l].Magnitude != 0.0) + if (s[l].Magnitude != 0.0) { for (j = l + 1; j < ncu; j++) { - t = -Cdotc(U, matrixCopy.RowCount, l, j, l) / U.At(l, l); + t = -Cdotc(u, matrixCopy.RowCount, l, j, l)/u.At(l, l); if (t != Complex.Zero) { for (var ii = l; ii < matrixCopy.RowCount; ii++) { - U.At(ii, j, U.At(ii, j) + (t * U.At(ii, l))); + u.At(ii, j, u.At(ii, j) + (t*u.At(ii, l))); } } } - CscalColumn(U, matrixCopy.RowCount, l, l, -1.0); - U.At(l, l, Complex.One + U.At(l, l)); + CscalColumn(u, matrixCopy.RowCount, l, l, -1.0); + u.At(l, l, Complex.One + u.At(l, l)); for (i = 0; i < l; i++) { - U.At(i, l, Complex.Zero); + u.At(i, l, Complex.Zero); } } else { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, l, Complex.Zero); + u.At(i, l, Complex.Zero); } - U.At(l, l, Complex.One); + u.At(l, l, Complex.One); } } } // If it is required, generate v. - if (ComputeVectors) + if (computeVectors) { for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) { @@ -290,12 +285,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (j = lp1; j < matrixCopy.ColumnCount; j++) { - t = -Cdotc(VT, matrixCopy.ColumnCount, l, j, lp1) / VT.At(lp1, l); + t = -Cdotc(vt, matrixCopy.ColumnCount, l, j, lp1)/vt.At(lp1, l); if (t != Complex.Zero) { for (var ii = l; ii < matrixCopy.ColumnCount; ii++) { - VT.At(ii, j, VT.At(ii, j) + (t * VT.At(ii, l))); + vt.At(ii, j, vt.At(ii, j) + (t*vt.At(ii, l))); } } } @@ -304,10 +299,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (i = 0; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, Complex.Zero); + vt.At(i, l, Complex.Zero); } - VT.At(l, l, Complex.One); + vt.At(l, l, Complex.One); } } @@ -315,19 +310,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (i = 0; i < m; i++) { Complex r; - if (S[i].Magnitude != 0.0) + if (s[i].Magnitude != 0.0) { - t = S[i].Magnitude; - r = S[i] / t; - S[i] = t; + t = s[i].Magnitude; + r = s[i]/t; + s[i] = t; if (i < m - 1) { - e[i] = e[i] / r; + e[i] = e[i]/r; } - if (ComputeVectors) + if (computeVectors) { - CscalColumn(U, matrixCopy.RowCount, i, 0, r); + CscalColumn(u, matrixCopy.RowCount, i, 0, r); } } @@ -340,12 +335,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization if (e[i].Magnitude != 0.0) { t = e[i].Magnitude; - r = t / e[i]; + r = t/e[i]; e[i] = t; - S[i + 1] = S[i + 1] * r; - if (ComputeVectors) + s[i + 1] = s[i + 1]*r; + if (computeVectors) { - CscalColumn(VT, matrixCopy.ColumnCount, i + 1, 0, r); + CscalColumn(vt, matrixCopy.ColumnCount, i + 1, 0, r); } } } @@ -373,7 +368,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization double test; for (l = m - 2; l >= 0; l--) { - test = S[l].Magnitude + S[l + 1].Magnitude; + test = s[l].Magnitude + s[l + 1].Magnitude; ztest = test + e[l].Magnitude; if (ztest.AlmostEqualInDecimalPlaces(test, 15)) { @@ -403,10 +398,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization test = test + e[ls - 1].Magnitude; } - ztest = test + S[ls].Magnitude; + ztest = test + s[ls].Magnitude; if (ztest.AlmostEqualInDecimalPlaces(test, 15)) { - S[ls] = Complex.Zero; + s[ls] = Complex.Zero; break; } } @@ -435,7 +430,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization double sn; switch (kase) { - // Deflate negligible VectorS[m]. + // Deflate negligible VectorS[m]. case 1: f = e[m - 2].Real; e[m - 2] = Complex.Zero; @@ -443,73 +438,73 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (var kk = l; kk < m - 1; kk++) { k = m - 2 - kk + l; - t1 = S[k].Real; + t1 = s[k].Real; Srotg(ref t1, ref f, out cs, out sn); - S[k] = t1; + s[k] = t1; if (k != l) { - f = -sn * e[k - 1].Real; - e[k - 1] = cs * e[k - 1]; + f = -sn*e[k - 1].Real; + e[k - 1] = cs*e[k - 1]; } - if (ComputeVectors) + if (computeVectors) { - Csrot(VT, matrixCopy.ColumnCount, k, m - 1, cs, sn); + Csrot(vt, matrixCopy.ColumnCount, k, m - 1, cs, sn); } } break; - // Split at negligible VectorS[l]. + // Split at negligible VectorS[l]. case 2: f = e[l - 1].Real; e[l - 1] = Complex.Zero; for (k = l; k < m; k++) { - t1 = S[k].Real; + t1 = s[k].Real; Srotg(ref t1, ref f, out cs, out sn); - S[k] = t1; - f = -sn * e[k].Real; - e[k] = cs * e[k]; - if (ComputeVectors) + s[k] = t1; + f = -sn*e[k].Real; + e[k] = cs*e[k]; + if (computeVectors) { - Csrot(U, matrixCopy.RowCount, k, l - 1, cs, sn); + Csrot(u, matrixCopy.RowCount, k, l - 1, cs, sn); } } break; - // Perform one qr step. + // Perform one qr step. case 3: // Calculate the shift. var scale = 0.0; - scale = Math.Max(scale, S[m - 1].Magnitude); - scale = Math.Max(scale, S[m - 2].Magnitude); + scale = Math.Max(scale, s[m - 1].Magnitude); + scale = Math.Max(scale, s[m - 2].Magnitude); scale = Math.Max(scale, e[m - 2].Magnitude); - scale = Math.Max(scale, S[l].Magnitude); + scale = Math.Max(scale, s[l].Magnitude); scale = Math.Max(scale, e[l].Magnitude); - var sm = S[m - 1].Real / scale; - var smm1 = S[m - 2].Real / scale; - var emm1 = e[m - 2].Real / scale; - var sl = S[l].Real / scale; - var el = e[l].Real / scale; - var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; - var c = (sm * emm1) * (sm * emm1); + var sm = s[m - 1].Real/scale; + var smm1 = s[m - 2].Real/scale; + var emm1 = e[m - 2].Real/scale; + var sl = s[l].Real/scale; + var el = e[l].Real/scale; + var b = (((smm1 + sm)*(smm1 - sm)) + (emm1*emm1))/2.0; + var c = (sm*emm1)*(sm*emm1); var shift = 0.0; if (b != 0.0 || c != 0.0) { - shift = Math.Sqrt((b * b) + c); + shift = Math.Sqrt((b*b) + c); if (b < 0.0) { shift = -shift; } - shift = c / (b + shift); + shift = c/(b + shift); } - f = ((sl + sm) * (sl - sm)) + shift; - var g = sl * el; + f = ((sl + sm)*(sl - sm)) + shift; + var g = sl*el; // Chase zeros. for (k = l; k < m - 1; k++) @@ -520,24 +515,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization e[k - 1] = f; } - f = (cs * S[k].Real) + (sn * e[k].Real); - e[k] = (cs * e[k]) - (sn * S[k]); - g = sn * S[k + 1].Real; - S[k + 1] = cs * S[k + 1]; - if (ComputeVectors) + f = (cs*s[k].Real) + (sn*e[k].Real); + e[k] = (cs*e[k]) - (sn*s[k]); + g = sn*s[k + 1].Real; + s[k + 1] = cs*s[k + 1]; + if (computeVectors) { - Csrot(VT, matrixCopy.ColumnCount, k, k + 1, cs, sn); + Csrot(vt, matrixCopy.ColumnCount, k, k + 1, cs, sn); } Srotg(ref f, ref g, out cs, out sn); - S[k] = f; - f = (cs * e[k].Real) + (sn * S[k + 1].Real); - S[k + 1] = (-sn * e[k]) + (cs * S[k + 1]); - g = sn * e[k + 1].Real; - e[k + 1] = cs * e[k + 1]; - if (ComputeVectors && k < matrixCopy.RowCount) + s[k] = f; + f = (cs*e[k].Real) + (sn*s[k + 1].Real); + s[k + 1] = (-sn*e[k]) + (cs*s[k + 1]); + g = sn*e[k + 1].Real; + e[k + 1] = cs*e[k + 1]; + if (computeVectors && k < matrixCopy.RowCount) { - Csrot(U, matrixCopy.RowCount, k, k + 1, cs, sn); + Csrot(u, matrixCopy.RowCount, k, k + 1, cs, sn); } } @@ -545,37 +540,37 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization iter = iter + 1; break; - // Convergence. + // Convergence. case 4: // Make the singular value positive - if (S[l].Real < 0.0) + if (s[l].Real < 0.0) { - S[l] = -S[l]; - if (ComputeVectors) + s[l] = -s[l]; + if (computeVectors) { - CscalColumn(VT, matrixCopy.ColumnCount, l, 0, -1.0); + CscalColumn(vt, matrixCopy.ColumnCount, l, 0, -1.0); } } // Order the singular value. while (l != mn - 1) { - if (S[l].Real >= S[l + 1].Real) + if (s[l].Real >= s[l + 1].Real) { break; } - t = S[l]; - S[l] = S[l + 1]; - S[l + 1] = t; - if (ComputeVectors && l < matrixCopy.ColumnCount) + t = s[l]; + s[l] = s[l + 1]; + s[l + 1] = t; + if (computeVectors && l < matrixCopy.ColumnCount) { - Swap(VT, matrixCopy.ColumnCount, l, l + 1); + Swap(vt, matrixCopy.ColumnCount, l, l + 1); } - if (ComputeVectors && l < matrixCopy.RowCount) + if (computeVectors && l < matrixCopy.RowCount) { - Swap(U, matrixCopy.RowCount, l, l + 1); + Swap(u, matrixCopy.RowCount, l, l + 1); } l = l + 1; @@ -587,9 +582,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } } - if (ComputeVectors) + if (computeVectors) { - VT = VT.ConjugateTranspose(); + vt = vt.ConjugateTranspose(); } // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses @@ -601,11 +596,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var tmp = matrixCopy.CreateVector(nm); for (i = 0; i < nm; i++) { - tmp[i] = S[i]; + tmp[i] = s[i]; } - S = tmp; + s = tmp; } + + return new UserSvd(s, u, vt, computeVectors); + } + + UserSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -614,9 +616,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Complex value z1 /// Complex value z2 /// Result multiplication of signum function and absolute value - private static Complex Csign(Complex z1, Complex z2) + static Complex Csign(Complex z1, Complex z2) { - return z1.Magnitude * (z2 / z2.Magnitude); + return z1.Magnitude*(z2/z2.Magnitude); } /// @@ -626,7 +628,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The number of rows in /// Column A index to swap /// Column B index to swap - private static void Swap(Matrix a, int rowCount, int columnA, int columnB) + static void Swap(Matrix a, int rowCount, int columnA, int columnB) { for (var i = 0; i < rowCount; i++) { @@ -644,11 +646,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Column to scale /// Row to scale from /// Scale value - private static void CscalColumn(Matrix a, int rowCount, int column, int rowStart, Complex z) + static void CscalColumn(Matrix a, int rowCount, int column, int rowStart, Complex z) { for (var i = rowStart; i < rowCount; i++) { - a.At(i, column, a.At(i, column) * z); + a.At(i, column, a.At(i, column)*z); } } @@ -658,11 +660,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Source vector /// Row to scale from /// Scale value - private static void CscalVector(Complex[] a, int start, Complex z) + static void CscalVector(Complex[] a, int start, Complex z) { for (var i = start; i < a.Length; i++) { - a[i] = a[i] * z; + a[i] = a[i]*z; } } @@ -675,7 +677,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Contains the parameter c associated with the Givens rotation /// Contains the parameter s associated with the Givens rotation /// This is equivalent to the DROTG LAPACK routine. - private static void Srotg(ref double da, ref double db, out double c, out double s) + static void Srotg(ref double da, ref double db, out double c, out double s) { double r, z; @@ -697,16 +699,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } else { - var sda = da / scale; - var sdb = db / scale; - r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); + var sda = da/scale; + var sdb = db/scale; + r = scale*Math.Sqrt((sda*sda) + (sdb*sdb)); if (roe < 0.0) { r = -r; } - c = da / r; - s = db / r; + c = da/r; + s = db/r; z = 1.0; if (absda > absdb) { @@ -715,7 +717,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization if (absdb >= absda && c != 0.0) { - z = 1.0 / c; + z = 1.0/c; } } @@ -731,12 +733,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Column index /// Start row index /// Norm2 (Euclidean norm) of the column - private static double Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) + static double Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) { var s = 0.0; for (var i = rowStart; i < rowCount; i++) { - s += a.At(i, column).Magnitude * a.At(i, column).Magnitude; + s += a.At(i, column).Magnitude*a.At(i, column).Magnitude; } return Math.Sqrt(s); @@ -748,12 +750,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Source vector /// Start index /// Norm2 (Euclidean norm) of the vector - private static double Cnrm2Vector(Complex[] a, int rowStart) + static double Cnrm2Vector(Complex[] a, int rowStart) { var s = 0.0; for (var i = rowStart; i < a.Length; i++) { - s += a[i].Magnitude * a[i].Magnitude; + s += a[i].Magnitude*a[i].Magnitude; } return Math.Sqrt(s); @@ -768,12 +770,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Index of column B /// Starting row index /// Dot product value - private static Complex Cdotc(Matrix a, int rowCount, int columnA, int columnB, int rowStart) + static Complex Cdotc(Matrix a, int rowCount, int columnA, int columnB, int rowStart) { var z = Complex.Zero; for (var i = rowStart; i < rowCount; i++) { - z += a.At(i, columnA).Conjugate() * a.At(i, columnB); + z += a.At(i, columnA).Conjugate()*a.At(i, columnB); } return z; @@ -789,17 +791,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Index of column B /// scalar cos value /// scalar sin value - private static void Csrot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) + static void Csrot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) { for (var i = 0; i < rowCount; i++) { - var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); - var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); + var z = (c*a.At(i, columnA)) + (s*a.At(i, columnB)); + var tmp = (c*a.At(i, columnB)) - (s*a.At(i, columnA)); a.At(i, columnB, tmp); a.At(i, columnA, z); } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -807,18 +809,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -855,7 +846,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j).Conjugate() * input.At(i, k); + value += U.At(i, j).Conjugate()*input.At(i, k); } value /= S[j]; @@ -869,7 +860,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var value = Complex.Zero; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j).Conjugate() * tmp[i]; + value += VT.At(i, j).Conjugate()*tmp[i]; } result.At(j, k, value); @@ -884,17 +875,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -921,7 +902,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j).Conjugate() * input[i]; + value += U.At(i, j).Conjugate()*input[i]; } value /= S[j]; @@ -935,7 +916,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var value = Complex.Zero; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j).Conjugate() * tmp[i]; + value += VT.At(i, j).Conjugate()*tmp[i]; } result[j] = value; diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index 2fde3cc5..a95f76bc 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -480,7 +480,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Svd Svd(bool computeVectors) { - return new UserSvd(this, computeVectors); + return UserSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index 1d6f10d1..4f247fee 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Svd Svd(bool computeVectors) { - return new DenseSvd(this, computeVectors); + return DenseSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs index dcf22c84..918b74c8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.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 { @@ -49,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public sealed class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -59,19 +59,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// If SVD algorithm failed to converge with matrix . - public DenseSvd(DenseMatrix matrix, bool computeVectors) + public static DenseSvd Create(DenseMatrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); - S = new DenseVector(nm); - U = new DenseMatrix(matrix.RowCount); - VT = new DenseMatrix(matrix.ColumnCount); - Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values); + var s = new DenseVector(nm); + var u = new DenseMatrix(matrix.RowCount); + var vt = new DenseMatrix(matrix.ColumnCount); + Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix) matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, s.Values, u.Values, vt.Values); + + return new DenseSvd(s, u, vt, computeVectors); + } + + DenseSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -81,18 +82,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -127,7 +117,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, input.ColumnCount, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values); } /// @@ -137,17 +127,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -177,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, 1, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs index 8a125b75..c6f5316f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.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,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.LinearAlgebra.Factorization; -using MathNet.Numerics.Properties; using System; using System.Linq; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { @@ -53,6 +53,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// public abstract class Svd : Svd { + protected Svd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { + } + /// /// Gets the effective numerical matrix rank. /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs index 7cce7edc..ce998c45 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.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 { @@ -49,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public sealed class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -59,20 +59,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// - public UserSvd(Matrix matrix, bool computeVectors) + public static UserSvd Create(Matrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); var matrixCopy = matrix.Clone(); - S = matrixCopy.CreateVector(nm); - U = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); - VT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); + var s = matrixCopy.CreateVector(nm); + var u = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); + var vt = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); const int maxiter = 1000; var e = new Complex32[matrixCopy.ColumnCount]; @@ -95,34 +89,34 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization if (l < nct) { // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. - S[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); - if (S[l].Magnitude != 0.0f) + s[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); + if (s[l].Magnitude != 0.0f) { if (matrixCopy.At(l, l).Magnitude != 0.0f) { - S[l] = Csign(S[l], matrixCopy.At(l, l)); + s[l] = Csign(s[l], matrixCopy.At(l, l)); } - CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f / S[l]); + CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f/s[l]); matrixCopy.At(l, l, (Complex32.One + matrixCopy.At(l, l))); } - S[l] = -S[l]; + s[l] = -s[l]; } for (j = lp1; j < matrixCopy.ColumnCount; j++) { if (l < nct) { - if (S[l].Magnitude != 0.0f) + if (s[l].Magnitude != 0.0f) { // Apply the transformation. - t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); + t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l)/matrixCopy.At(l, l); if (t != Complex32.Zero) { for (var ii = l; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t*matrixCopy.At(ii, l))); } } } @@ -133,12 +127,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization e[j] = matrixCopy.At(l, j).Conjugate(); } - if (ComputeVectors && l < nct) + if (computeVectors && l < nct) { // Place the transformation in u for subsequent back multiplication. for (i = l; i < matrixCopy.RowCount; i++) { - U.At(i, l, matrixCopy.At(i, l)); + u.At(i, l, matrixCopy.At(i, l)); } } @@ -157,7 +151,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization e[l] = Csign(e[l], e[lp1]); } - CscalVector(e, lp1, 1.0f / e[l]); + CscalVector(e, lp1, 1.0f/e[l]); e[lp1] = Complex32.One + e[lp1]; } @@ -176,30 +170,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - work[ii] += e[j] * matrixCopy.At(ii, j); + work[ii] += e[j]*matrixCopy.At(ii, j); } } } for (j = lp1; j < matrixCopy.ColumnCount; j++) { - var ww = (-e[j] / e[lp1]).Conjugate(); + var ww = (-e[j]/e[lp1]).Conjugate(); if (ww != Complex32.Zero) { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww*work[ii])); } } } } - if (ComputeVectors) + if (computeVectors) { // Place the transformation in v for subsequent back multiplication. for (i = lp1; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, e[i]); + vt.At(i, l, e[i]); } } } @@ -210,12 +204,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var nrtp1 = nrt + 1; if (nct < matrixCopy.ColumnCount) { - S[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); + s[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); } if (matrixCopy.RowCount < m) { - S[m - 1] = Complex32.Zero; + s[m - 1] = Complex32.Zero; } if (nrtp1 < m) @@ -226,55 +220,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization e[m - 1] = Complex32.Zero; // If required, generate u. - if (ComputeVectors) + if (computeVectors) { for (j = nctp1 - 1; j < ncu; j++) { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, j, Complex32.Zero); + u.At(i, j, Complex32.Zero); } - U.At(j, j, Complex32.One); + u.At(j, j, Complex32.One); } for (l = nct - 1; l >= 0; l--) { - if (S[l].Magnitude != 0.0f) + if (s[l].Magnitude != 0.0f) { for (j = l + 1; j < ncu; j++) { - t = -Cdotc(U, matrixCopy.RowCount, l, j, l) / U.At(l, l); + t = -Cdotc(u, matrixCopy.RowCount, l, j, l)/u.At(l, l); if (t != Complex32.Zero) { for (var ii = l; ii < matrixCopy.RowCount; ii++) { - U.At(ii, j, U.At(ii, j) + (t * U.At(ii, l))); + u.At(ii, j, u.At(ii, j) + (t*u.At(ii, l))); } } } - CscalColumn(U, matrixCopy.RowCount, l, l, -1.0f); - U.At(l, l, Complex32.One + U.At(l, l)); + CscalColumn(u, matrixCopy.RowCount, l, l, -1.0f); + u.At(l, l, Complex32.One + u.At(l, l)); for (i = 0; i < l; i++) { - U.At(i, l, Complex32.Zero); + u.At(i, l, Complex32.Zero); } } else { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, l, Complex32.Zero); + u.At(i, l, Complex32.Zero); } - U.At(l, l, Complex32.One); + u.At(l, l, Complex32.One); } } } // If it is required, generate v. - if (ComputeVectors) + if (computeVectors) { for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) { @@ -285,12 +279,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (j = lp1; j < matrixCopy.ColumnCount; j++) { - t = -Cdotc(VT, matrixCopy.ColumnCount, l, j, lp1) / VT.At(lp1, l); + t = -Cdotc(vt, matrixCopy.ColumnCount, l, j, lp1)/vt.At(lp1, l); if (t != Complex32.Zero) { for (var ii = l; ii < matrixCopy.ColumnCount; ii++) { - VT.At(ii, j, VT.At(ii, j) + (t * VT.At(ii, l))); + vt.At(ii, j, vt.At(ii, j) + (t*vt.At(ii, l))); } } } @@ -299,10 +293,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (i = 0; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, Complex32.Zero); + vt.At(i, l, Complex32.Zero); } - VT.At(l, l, Complex32.One); + vt.At(l, l, Complex32.One); } } @@ -310,19 +304,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (i = 0; i < m; i++) { Complex32 r; - if (S[i].Magnitude != 0.0f) + if (s[i].Magnitude != 0.0f) { - t = S[i].Magnitude; - r = S[i] / t; - S[i] = t; + t = s[i].Magnitude; + r = s[i]/t; + s[i] = t; if (i < m - 1) { - e[i] = e[i] / r; + e[i] = e[i]/r; } - if (ComputeVectors) + if (computeVectors) { - CscalColumn(U, matrixCopy.RowCount, i, 0, r); + CscalColumn(u, matrixCopy.RowCount, i, 0, r); } } @@ -335,12 +329,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization if (e[i].Magnitude != 0.0f) { t = e[i].Magnitude; - r = t / e[i]; + r = t/e[i]; e[i] = t; - S[i + 1] = S[i + 1] * r; - if (ComputeVectors) + s[i + 1] = s[i + 1]*r; + if (computeVectors) { - CscalColumn(VT, matrixCopy.ColumnCount, i + 1, 0, r); + CscalColumn(vt, matrixCopy.ColumnCount, i + 1, 0, r); } } } @@ -368,7 +362,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization float test; for (l = m - 2; l >= 0; l--) { - test = S[l].Magnitude + S[l + 1].Magnitude; + test = s[l].Magnitude + s[l + 1].Magnitude; ztest = test + e[l].Magnitude; if (ztest.AlmostEqualInDecimalPlaces(test, 7)) { @@ -398,10 +392,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization test = test + e[ls - 1].Magnitude; } - ztest = test + S[ls].Magnitude; + ztest = test + s[ls].Magnitude; if (ztest.AlmostEqualInDecimalPlaces(test, 7)) { - S[ls] = Complex32.Zero; + s[ls] = Complex32.Zero; break; } } @@ -430,7 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization float cs; switch (kase) { - // Deflate negligible VectorS[m]. + // Deflate negligible VectorS[m]. case 1: f = e[m - 2].Real; e[m - 2] = Complex32.Zero; @@ -438,73 +432,73 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (var kk = l; kk < m - 1; kk++) { k = m - 2 - kk + l; - t1 = S[k].Real; + t1 = s[k].Real; Srotg(ref t1, ref f, out cs, out sn); - S[k] = t1; + s[k] = t1; if (k != l) { - f = -sn * e[k - 1].Real; - e[k - 1] = cs * e[k - 1]; + f = -sn*e[k - 1].Real; + e[k - 1] = cs*e[k - 1]; } - if (ComputeVectors) + if (computeVectors) { - Csrot(VT, matrixCopy.ColumnCount, k, m - 1, cs, sn); + Csrot(vt, matrixCopy.ColumnCount, k, m - 1, cs, sn); } } break; - // Split at negligible VectorS[l]. + // Split at negligible VectorS[l]. case 2: f = e[l - 1].Real; e[l - 1] = Complex32.Zero; for (k = l; k < m; k++) { - t1 = S[k].Real; + t1 = s[k].Real; Srotg(ref t1, ref f, out cs, out sn); - S[k] = t1; - f = -sn * e[k].Real; - e[k] = cs * e[k]; - if (ComputeVectors) + s[k] = t1; + f = -sn*e[k].Real; + e[k] = cs*e[k]; + if (computeVectors) { - Csrot(U, matrixCopy.RowCount, k, l - 1, cs, sn); + Csrot(u, matrixCopy.RowCount, k, l - 1, cs, sn); } } break; - // Perform one qr step. + // Perform one qr step. case 3: // Calculate the shift. var scale = 0.0f; - scale = Math.Max(scale, S[m - 1].Magnitude); - scale = Math.Max(scale, S[m - 2].Magnitude); + scale = Math.Max(scale, s[m - 1].Magnitude); + scale = Math.Max(scale, s[m - 2].Magnitude); scale = Math.Max(scale, e[m - 2].Magnitude); - scale = Math.Max(scale, S[l].Magnitude); + scale = Math.Max(scale, s[l].Magnitude); scale = Math.Max(scale, e[l].Magnitude); - var sm = S[m - 1].Real / scale; - var smm1 = S[m - 2].Real / scale; - var emm1 = e[m - 2].Real / scale; - var sl = S[l].Real / scale; - var el = e[l].Real / scale; - var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0f; - var c = (sm * emm1) * (sm * emm1); + var sm = s[m - 1].Real/scale; + var smm1 = s[m - 2].Real/scale; + var emm1 = e[m - 2].Real/scale; + var sl = s[l].Real/scale; + var el = e[l].Real/scale; + var b = (((smm1 + sm)*(smm1 - sm)) + (emm1*emm1))/2.0f; + var c = (sm*emm1)*(sm*emm1); var shift = 0.0f; if (b != 0.0f || c != 0.0f) { - shift = (float)Math.Sqrt((b * b) + c); + shift = (float) Math.Sqrt((b*b) + c); if (b < 0.0f) { shift = -shift; } - shift = c / (b + shift); + shift = c/(b + shift); } - f = ((sl + sm) * (sl - sm)) + shift; - var g = sl * el; + f = ((sl + sm)*(sl - sm)) + shift; + var g = sl*el; // Chase zeros. for (k = l; k < m - 1; k++) @@ -515,24 +509,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization e[k - 1] = f; } - f = (cs * S[k].Real) + (sn * e[k].Real); - e[k] = (cs * e[k]) - (sn * S[k]); - g = sn * S[k + 1].Real; - S[k + 1] = cs * S[k + 1]; - if (ComputeVectors) + f = (cs*s[k].Real) + (sn*e[k].Real); + e[k] = (cs*e[k]) - (sn*s[k]); + g = sn*s[k + 1].Real; + s[k + 1] = cs*s[k + 1]; + if (computeVectors) { - Csrot(VT, matrixCopy.ColumnCount, k, k + 1, cs, sn); + Csrot(vt, matrixCopy.ColumnCount, k, k + 1, cs, sn); } Srotg(ref f, ref g, out cs, out sn); - S[k] = f; - f = (cs * e[k].Real) + (sn * S[k + 1].Real); - S[k + 1] = (-sn * e[k]) + (cs * S[k + 1]); - g = sn * e[k + 1].Real; - e[k + 1] = cs * e[k + 1]; - if (ComputeVectors && k < matrixCopy.RowCount) + s[k] = f; + f = (cs*e[k].Real) + (sn*s[k + 1].Real); + s[k + 1] = (-sn*e[k]) + (cs*s[k + 1]); + g = sn*e[k + 1].Real; + e[k + 1] = cs*e[k + 1]; + if (computeVectors && k < matrixCopy.RowCount) { - Csrot(U, matrixCopy.RowCount, k, k + 1, cs, sn); + Csrot(u, matrixCopy.RowCount, k, k + 1, cs, sn); } } @@ -540,37 +534,37 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization iter = iter + 1; break; - // Convergence. + // Convergence. case 4: // Make the singular value positive - if (S[l].Real < 0.0f) + if (s[l].Real < 0.0f) { - S[l] = -S[l]; - if (ComputeVectors) + s[l] = -s[l]; + if (computeVectors) { - CscalColumn(VT, matrixCopy.ColumnCount, l, 0, -1.0f); + CscalColumn(vt, matrixCopy.ColumnCount, l, 0, -1.0f); } } // Order the singular value. while (l != mn - 1) { - if (S[l].Real >= S[l + 1].Real) + if (s[l].Real >= s[l + 1].Real) { break; } - t = S[l]; - S[l] = S[l + 1]; - S[l + 1] = t; - if (ComputeVectors && l < matrixCopy.ColumnCount) + t = s[l]; + s[l] = s[l + 1]; + s[l + 1] = t; + if (computeVectors && l < matrixCopy.ColumnCount) { - Swap(VT, matrixCopy.ColumnCount, l, l + 1); + Swap(vt, matrixCopy.ColumnCount, l, l + 1); } - if (ComputeVectors && l < matrixCopy.RowCount) + if (computeVectors && l < matrixCopy.RowCount) { - Swap(U, matrixCopy.RowCount, l, l + 1); + Swap(u, matrixCopy.RowCount, l, l + 1); } l = l + 1; @@ -582,9 +576,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } } - if (ComputeVectors) + if (computeVectors) { - VT = VT.ConjugateTranspose(); + vt = vt.ConjugateTranspose(); } // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses @@ -596,11 +590,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var tmp = matrixCopy.CreateVector(nm); for (i = 0; i < nm; i++) { - tmp[i] = S[i]; + tmp[i] = s[i]; } - S = tmp; + s = tmp; } + + return new UserSvd(s, u, vt, computeVectors); + } + + UserSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -609,9 +610,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Complex32 value z1 /// Complex32 value z2 /// Result multiplication of signum function and absolute value - private static Complex32 Csign(Complex32 z1, Complex32 z2) + static Complex32 Csign(Complex32 z1, Complex32 z2) { - return z1.Magnitude * (z2 / z2.Magnitude); + return z1.Magnitude*(z2/z2.Magnitude); } /// @@ -621,7 +622,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The number of rows in /// Column A index to swap /// Column B index to swap - private static void Swap(Matrix a, int rowCount, int columnA, int columnB) + static void Swap(Matrix a, int rowCount, int columnA, int columnB) { for (var i = 0; i < rowCount; i++) { @@ -639,11 +640,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Column to scale /// Row to scale from /// Scale value - private static void CscalColumn(Matrix a, int rowCount, int column, int rowStart, Complex32 z) + static void CscalColumn(Matrix a, int rowCount, int column, int rowStart, Complex32 z) { for (var i = rowStart; i < rowCount; i++) { - a.At(i, column, a.At(i, column) * z); + a.At(i, column, a.At(i, column)*z); } } @@ -653,11 +654,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Source vector /// Row to scale from /// Scale value - private static void CscalVector(Complex32[] a, int start, Complex32 z) + static void CscalVector(Complex32[] a, int start, Complex32 z) { for (var i = start; i < a.Length; i++) { - a[i] = a[i] * z; + a[i] = a[i]*z; } } @@ -670,7 +671,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Contains the parameter c associated with the Givens rotation /// Contains the parameter s associated with the Givens rotation /// This is equivalent to the DROTG LAPACK routine. - private static void Srotg(ref float da, ref float db, out float c, out float s) + static void Srotg(ref float da, ref float db, out float c, out float s) { float r, z; @@ -692,16 +693,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } else { - var sda = da / scale; - var sdb = db / scale; - r = scale * (float)Math.Sqrt((sda * sda) + (sdb * sdb)); + var sda = da/scale; + var sdb = db/scale; + r = scale*(float) Math.Sqrt((sda*sda) + (sdb*sdb)); if (roe < 0.0f) { r = -r; } - c = da / r; - s = db / r; + c = da/r; + s = db/r; z = 1.0f; if (absda > absdb) { @@ -710,7 +711,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization if (absdb >= absda && c != 0.0f) { - z = 1.0f / c; + z = 1.0f/c; } } @@ -726,15 +727,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Column index /// Start row index /// Norm2 (Euclidean norm) of the column - private static float Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) + static float Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) { var s = 0.0f; for (var i = rowStart; i < rowCount; i++) { - s += a.At(i, column).Magnitude * a.At(i, column).Magnitude; + s += a.At(i, column).Magnitude*a.At(i, column).Magnitude; } - return (float)Math.Sqrt(s); + return (float) Math.Sqrt(s); } /// @@ -743,15 +744,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Source vector /// Start index /// Norm2 (Euclidean norm) of the vector - private static float Cnrm2Vector(Complex32[] a, int rowStart) + static float Cnrm2Vector(Complex32[] a, int rowStart) { var s = 0.0f; for (var i = rowStart; i < a.Length; i++) { - s += a[i].Magnitude * a[i].Magnitude; + s += a[i].Magnitude*a[i].Magnitude; } - return (float)Math.Sqrt(s); + return (float) Math.Sqrt(s); } /// @@ -763,12 +764,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Index of column B /// Starting row index /// Dot product value - private static Complex32 Cdotc(Matrix a, int rowCount, int columnA, int columnB, int rowStart) + static Complex32 Cdotc(Matrix a, int rowCount, int columnA, int columnB, int rowStart) { var z = Complex32.Zero; for (var i = rowStart; i < rowCount; i++) { - z += a.At(i, columnA).Conjugate() * a.At(i, columnB); + z += a.At(i, columnA).Conjugate()*a.At(i, columnB); } return z; @@ -784,17 +785,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Index of column B /// scalar cos value /// scalar sin value - private static void Csrot(Matrix a, int rowCount, int columnA, int columnB, float c, float s) + static void Csrot(Matrix a, int rowCount, int columnA, int columnB, float c, float s) { for (var i = 0; i < rowCount; i++) { - var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); - var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); + var z = (c*a.At(i, columnA)) + (s*a.At(i, columnB)); + var tmp = (c*a.At(i, columnB)) - (s*a.At(i, columnA)); a.At(i, columnB, tmp); a.At(i, columnA, z); } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -802,18 +803,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -850,7 +840,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j).Conjugate() * input.At(i, k); + value += U.At(i, j).Conjugate()*input.At(i, k); } value /= S[j]; @@ -864,7 +854,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var value = Complex32.Zero; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j).Conjugate() * tmp[i]; + value += VT.At(i, j).Conjugate()*tmp[i]; } result.At(j, k, value); @@ -879,17 +869,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -916,7 +896,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j).Conjugate() * input[i]; + value += U.At(i, j).Conjugate()*input[i]; } value /= S[j]; @@ -930,7 +910,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var value = Complex32.Zero; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j).Conjugate() * tmp[i]; + value += VT.At(i, j).Conjugate()*tmp[i]; } result[j] = value; diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index 156964b6..189f9aa3 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -475,7 +475,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Svd Svd(bool computeVectors) { - return new UserSvd(this, computeVectors); + return UserSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index 87d3480e..de183c98 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -1057,7 +1057,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Svd Svd(bool computeVectors) { - return new DenseSvd(this, computeVectors); + return DenseSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs index fa98016c..510a1767 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.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 { @@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public sealed class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,19 +57,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// If SVD algorithm failed to converge with matrix . - public DenseSvd(DenseMatrix matrix, bool computeVectors) + public static DenseSvd Create(DenseMatrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); - S = new DenseVector(nm); - U = new DenseMatrix(matrix.RowCount); - VT = new DenseMatrix(matrix.ColumnCount); - Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values); + var s = new DenseVector(nm); + var u = new DenseMatrix(matrix.RowCount); + var vt = new DenseMatrix(matrix.ColumnCount); + Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix) matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, s.Values, u.Values, vt.Values); + + return new DenseSvd(s, u, vt, computeVectors); + } + + DenseSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -79,18 +80,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -125,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, input.ColumnCount, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values); } /// @@ -135,17 +125,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -175,7 +155,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, 1, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs index bbb99648..95a103f8 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.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,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.LinearAlgebra.Factorization; -using MathNet.Numerics.Properties; using System; using System.Linq; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { @@ -51,6 +51,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// public abstract class Svd : Svd { + protected Svd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { + } + /// /// Gets the effective numerical matrix rank. /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs index 06cae758..feb9c62d 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.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 { @@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public sealed class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,20 +57,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// - public UserSvd(Matrix matrix, bool computeVectors) + public static UserSvd Create(Matrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); var matrixCopy = matrix.Clone(); - S = matrixCopy.CreateVector(nm); - U = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); - VT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); + var s = matrixCopy.CreateVector(nm); + var u = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); + var vt = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); const int maxiter = 1000; var e = new double[matrixCopy.ColumnCount]; @@ -94,32 +88,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); - S[l] = xnorm; - if (S[l] != 0.0) + s[l] = xnorm; + if (s[l] != 0.0) { if (matrixCopy.At(l, l) != 0.0) { - S[l] = Dsign(S[l], matrixCopy.At(l, l)); + s[l] = Dsign(s[l], matrixCopy.At(l, l)); } - DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / S[l]); + DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0/s[l]); matrixCopy.At(l, l, (1.0 + matrixCopy.At(l, l))); } - S[l] = -S[l]; + s[l] = -s[l]; } for (j = lp1; j < matrixCopy.ColumnCount; j++) { if (l < nct) { - if (S[l] != 0.0) + if (s[l] != 0.0) { // Apply the transformation. - t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); + t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l)/matrixCopy.At(l, l); for (var ii = l; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t*matrixCopy.At(ii, l))); } } } @@ -129,12 +123,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization e[j] = matrixCopy.At(l, j); } - if (ComputeVectors && l < nct) + if (computeVectors && l < nct) { // Place the transformation in u for subsequent back multiplication. for (i = l; i < matrixCopy.RowCount; i++) { - U.At(i, l, matrixCopy.At(i, l)); + u.At(i, l, matrixCopy.At(i, l)); } } @@ -153,7 +147,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization e[l] = Dsign(e[l], e[lp1]); } - DscalVector(e, lp1, 1.0 / e[l]); + DscalVector(e, lp1, 1.0/e[l]); e[lp1] = 1.0 + e[lp1]; } @@ -170,26 +164,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - work[ii] += e[j] * matrixCopy.At(ii, j); + work[ii] += e[j]*matrixCopy.At(ii, j); } } for (j = lp1; j < matrixCopy.ColumnCount; j++) { - var ww = -e[j] / e[lp1]; + var ww = -e[j]/e[lp1]; for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww*work[ii])); } } } - if (ComputeVectors) + if (computeVectors) { // Place the transformation in v for subsequent back multiplication. for (i = lp1; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, e[i]); + vt.At(i, l, e[i]); } } } @@ -200,12 +194,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var nrtp1 = nrt + 1; if (nct < matrixCopy.ColumnCount) { - S[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); + s[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); } if (matrixCopy.RowCount < m) { - S[m - 1] = 0.0; + s[m - 1] = 0.0; } if (nrtp1 < m) @@ -216,52 +210,52 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization e[m - 1] = 0.0; // If required, generate u. - if (ComputeVectors) + if (computeVectors) { for (j = nctp1 - 1; j < ncu; j++) { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, j, 0.0); + u.At(i, j, 0.0); } - U.At(j, j, 1.0); + u.At(j, j, 1.0); } for (l = nct - 1; l >= 0; l--) { - if (S[l] != 0.0) + if (s[l] != 0.0) { for (j = l + 1; j < ncu; j++) { - t = -Ddot(U, matrixCopy.RowCount, l, j, l) / U.At(l, l); + t = -Ddot(u, matrixCopy.RowCount, l, j, l)/u.At(l, l); for (var ii = l; ii < matrixCopy.RowCount; ii++) { - U.At(ii, j, U.At(ii, j) + (t * U.At(ii, l))); + u.At(ii, j, u.At(ii, j) + (t*u.At(ii, l))); } } - DscalColumn(U, matrixCopy.RowCount, l, l, -1.0); - U.At(l, l, 1.0 + U.At(l, l)); + DscalColumn(u, matrixCopy.RowCount, l, l, -1.0); + u.At(l, l, 1.0 + u.At(l, l)); for (i = 0; i < l; i++) { - U.At(i, l, 0.0); + u.At(i, l, 0.0); } } else { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, l, 0.0); + u.At(i, l, 0.0); } - U.At(l, l, 1.0); + u.At(l, l, 1.0); } } } // If it is required, generate v. - if (ComputeVectors) + if (computeVectors) { for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) { @@ -272,10 +266,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (j = lp1; j < matrixCopy.ColumnCount; j++) { - t = -Ddot(VT, matrixCopy.ColumnCount, l, j, lp1) / VT.At(lp1, l); + t = -Ddot(vt, matrixCopy.ColumnCount, l, j, lp1)/vt.At(lp1, l); for (var ii = l; ii < matrixCopy.ColumnCount; ii++) { - VT.At(ii, j, VT.At(ii, j) + (t * VT.At(ii, l))); + vt.At(ii, j, vt.At(ii, j) + (t*vt.At(ii, l))); } } } @@ -283,10 +277,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (i = 0; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, 0.0); + vt.At(i, l, 0.0); } - VT.At(l, l, 1.0); + vt.At(l, l, 1.0); } } @@ -294,19 +288,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (i = 0; i < m; i++) { double r; - if (S[i] != 0.0) + if (s[i] != 0.0) { - t = S[i]; - r = S[i] / t; - S[i] = t; + t = s[i]; + r = s[i]/t; + s[i] = t; if (i < m - 1) { - e[i] = e[i] / r; + e[i] = e[i]/r; } - if (ComputeVectors) + if (computeVectors) { - DscalColumn(U, matrixCopy.RowCount, i, 0, r); + DscalColumn(u, matrixCopy.RowCount, i, 0, r); } } @@ -319,12 +313,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization if (e[i] != 0.0) { t = e[i]; - r = t / e[i]; + r = t/e[i]; e[i] = t; - S[i + 1] = S[i + 1] * r; - if (ComputeVectors) + s[i + 1] = s[i + 1]*r; + if (computeVectors) { - DscalColumn(VT, matrixCopy.ColumnCount, i + 1, 0, r); + DscalColumn(vt, matrixCopy.ColumnCount, i + 1, 0, r); } } } @@ -352,7 +346,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization double test; for (l = m - 2; l >= 0; l--) { - test = Math.Abs(S[l]) + Math.Abs(S[l + 1]); + test = Math.Abs(s[l]) + Math.Abs(s[l + 1]); ztest = test + Math.Abs(e[l]); if (ztest.AlmostEqualInDecimalPlaces(test, 15)) { @@ -382,10 +376,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization test = test + Math.Abs(e[ls - 1]); } - ztest = test + Math.Abs(S[ls]); + ztest = test + Math.Abs(s[ls]); if (ztest.AlmostEqualInDecimalPlaces(test, 15)) { - S[ls] = 0.0; + s[ls] = 0.0; break; } } @@ -414,7 +408,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization double cs; switch (kase) { - // Deflate negligible VectorS[m]. + // Deflate negligible VectorS[m]. case 1: f = e[m - 2]; e[m - 2] = 0.0; @@ -422,72 +416,72 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var kk = l; kk < m - 1; kk++) { k = m - 2 - kk + l; - t1 = S[k]; + t1 = s[k]; Drotg(ref t1, ref f, out cs, out sn); - S[k] = t1; + s[k] = t1; if (k != l) { - f = -sn * e[k - 1]; - e[k - 1] = cs * e[k - 1]; + f = -sn*e[k - 1]; + e[k - 1] = cs*e[k - 1]; } - if (ComputeVectors) + if (computeVectors) { - Drot(VT, matrixCopy.ColumnCount, k, m - 1, cs, sn); + Drot(vt, matrixCopy.ColumnCount, k, m - 1, cs, sn); } } break; - // Split at negligible VectorS[l]. + // Split at negligible VectorS[l]. case 2: f = e[l - 1]; e[l - 1] = 0.0; for (k = l; k < m; k++) { - t1 = S[k]; + t1 = s[k]; Drotg(ref t1, ref f, out cs, out sn); - S[k] = t1; - f = -sn * e[k]; - e[k] = cs * e[k]; - if (ComputeVectors) + s[k] = t1; + f = -sn*e[k]; + e[k] = cs*e[k]; + if (computeVectors) { - Drot(U, matrixCopy.RowCount, k, l - 1, cs, sn); + Drot(u, matrixCopy.RowCount, k, l - 1, cs, sn); } } break; - // Perform one qr step. + // Perform one qr step. case 3: // Calculate the shift. var scale = 0.0; - scale = Math.Max(scale, Math.Abs(S[m - 1])); - scale = Math.Max(scale, Math.Abs(S[m - 2])); + scale = Math.Max(scale, Math.Abs(s[m - 1])); + scale = Math.Max(scale, Math.Abs(s[m - 2])); scale = Math.Max(scale, Math.Abs(e[m - 2])); - scale = Math.Max(scale, Math.Abs(S[l])); + scale = Math.Max(scale, Math.Abs(s[l])); scale = Math.Max(scale, Math.Abs(e[l])); - var sm = S[m - 1] / scale; - var smm1 = S[m - 2] / scale; - var emm1 = e[m - 2] / scale; - var sl = S[l] / scale; - var el = e[l] / scale; - var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; - var c = (sm * emm1) * (sm * emm1); + var sm = s[m - 1]/scale; + var smm1 = s[m - 2]/scale; + var emm1 = e[m - 2]/scale; + var sl = s[l]/scale; + var el = e[l]/scale; + var b = (((smm1 + sm)*(smm1 - sm)) + (emm1*emm1))/2.0; + var c = (sm*emm1)*(sm*emm1); var shift = 0.0; if (b != 0.0 || c != 0.0) { - shift = Math.Sqrt((b * b) + c); + shift = Math.Sqrt((b*b) + c); if (b < 0.0) { shift = -shift; } - shift = c / (b + shift); + shift = c/(b + shift); } - f = ((sl + sm) * (sl - sm)) + shift; - var g = sl * el; + f = ((sl + sm)*(sl - sm)) + shift; + var g = sl*el; // Chase zeros. for (k = l; k < m - 1; k++) @@ -498,24 +492,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization e[k - 1] = f; } - f = (cs * S[k]) + (sn * e[k]); - e[k] = (cs * e[k]) - (sn * S[k]); - g = sn * S[k + 1]; - S[k + 1] = cs * S[k + 1]; - if (ComputeVectors) + f = (cs*s[k]) + (sn*e[k]); + e[k] = (cs*e[k]) - (sn*s[k]); + g = sn*s[k + 1]; + s[k + 1] = cs*s[k + 1]; + if (computeVectors) { - Drot(VT, matrixCopy.ColumnCount, k, k + 1, cs, sn); + Drot(vt, matrixCopy.ColumnCount, k, k + 1, cs, sn); } Drotg(ref f, ref g, out cs, out sn); - S[k] = f; - f = (cs * e[k]) + (sn * S[k + 1]); - S[k + 1] = (-sn * e[k]) + (cs * S[k + 1]); - g = sn * e[k + 1]; - e[k + 1] = cs * e[k + 1]; - if (ComputeVectors && k < matrixCopy.RowCount) + s[k] = f; + f = (cs*e[k]) + (sn*s[k + 1]); + s[k + 1] = (-sn*e[k]) + (cs*s[k + 1]); + g = sn*e[k + 1]; + e[k + 1] = cs*e[k + 1]; + if (computeVectors && k < matrixCopy.RowCount) { - Drot(U, matrixCopy.RowCount, k, k + 1, cs, sn); + Drot(u, matrixCopy.RowCount, k, k + 1, cs, sn); } } @@ -523,37 +517,37 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization iter = iter + 1; break; - // Convergence. + // Convergence. case 4: // Make the singular value positive - if (S[l] < 0.0) + if (s[l] < 0.0) { - S[l] = -S[l]; - if (ComputeVectors) + s[l] = -s[l]; + if (computeVectors) { - DscalColumn(VT, matrixCopy.ColumnCount, l, 0, -1.0); + DscalColumn(vt, matrixCopy.ColumnCount, l, 0, -1.0); } } // Order the singular value. while (l != mn - 1) { - if (S[l] >= S[l + 1]) + if (s[l] >= s[l + 1]) { break; } - t = S[l]; - S[l] = S[l + 1]; - S[l + 1] = t; - if (ComputeVectors && l < matrixCopy.ColumnCount) + t = s[l]; + s[l] = s[l + 1]; + s[l + 1] = t; + if (computeVectors && l < matrixCopy.ColumnCount) { - Dswap(VT, matrixCopy.ColumnCount, l, l + 1); + Dswap(vt, matrixCopy.ColumnCount, l, l + 1); } - if (ComputeVectors && l < matrixCopy.RowCount) + if (computeVectors && l < matrixCopy.RowCount) { - Dswap(U, matrixCopy.RowCount, l, l + 1); + Dswap(u, matrixCopy.RowCount, l, l + 1); } l = l + 1; @@ -565,9 +559,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } } - if (ComputeVectors) + if (computeVectors) { - VT = VT.Transpose(); + vt = vt.Transpose(); } // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses @@ -579,11 +573,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var tmp = matrixCopy.CreateVector(nm); for (i = 0; i < nm; i++) { - tmp[i] = S[i]; + tmp[i] = s[i]; } - S = tmp; + s = tmp; } + + return new UserSvd(s, u, vt, computeVectors); + } + + UserSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -592,9 +593,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Double value z1 /// Double value z2 /// Result multiplication of signum function and absolute value - private static double Dsign(double z1, double z2) + static double Dsign(double z1, double z2) { - return Math.Abs(z1) * (z2 / Math.Abs(z2)); + return Math.Abs(z1)*(z2/Math.Abs(z2)); } /// @@ -604,7 +605,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The number of rows in /// Column A index to swap /// Column B index to swap - private static void Dswap(Matrix a, int rowCount, int columnA, int columnB) + static void Dswap(Matrix a, int rowCount, int columnA, int columnB) { for (var i = 0; i < rowCount; i++) { @@ -622,11 +623,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Column to scale /// Row to scale from /// Scale value - private static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, double z) + static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, double z) { for (var i = rowStart; i < rowCount; i++) { - a.At(i, column, a.At(i, column) * z); + a.At(i, column, a.At(i, column)*z); } } @@ -636,11 +637,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Source vector /// Row to scale from /// Scale value - private static void DscalVector(double[] a, int start, double z) + static void DscalVector(double[] a, int start, double z) { for (var i = start; i < a.Length; i++) { - a[i] = a[i] * z; + a[i] = a[i]*z; } } @@ -653,7 +654,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Contains the parameter c associated with the Givens rotation /// Contains the parameter s associated with the Givens rotation /// This is equivalent to the DROTG LAPACK routine. - private static void Drotg(ref double da, ref double db, out double c, out double s) + static void Drotg(ref double da, ref double db, out double c, out double s) { double r, z; @@ -675,16 +676,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } else { - var sda = da / scale; - var sdb = db / scale; - r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); + var sda = da/scale; + var sdb = db/scale; + r = scale*Math.Sqrt((sda*sda) + (sdb*sdb)); if (roe < 0.0) { r = -r; } - c = da / r; - s = db / r; + c = da/r; + s = db/r; z = 1.0; if (absda > absdb) { @@ -693,7 +694,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization if (absdb >= absda && c != 0.0) { - z = 1.0 / c; + z = 1.0/c; } } @@ -709,12 +710,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Column index /// Start row index /// Norm2 (Euclidean norm) of the column - private static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) + static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) { double s = 0; for (var i = rowStart; i < rowCount; i++) { - s += a.At(i, column) * a.At(i, column); + s += a.At(i, column)*a.At(i, column); } return Math.Sqrt(s); @@ -726,12 +727,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Source vector /// Start index /// Norm2 (Euclidean norm) of the vector - private static double Dnrm2Vector(double[] a, int rowStart) + static double Dnrm2Vector(double[] a, int rowStart) { double s = 0; for (var i = rowStart; i < a.Length; i++) { - s += a[i] * a[i]; + s += a[i]*a[i]; } return Math.Sqrt(s); @@ -746,12 +747,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Index of column B /// Starting row index /// Dot product value - private static double Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) + static double Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) { var z = 0.0; for (var i = rowStart; i < rowCount; i++) { - z += a.At(i, columnB) * a.At(i, columnA); + z += a.At(i, columnB)*a.At(i, columnA); } return z; @@ -767,17 +768,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Index of column B /// Scalar "c" value /// Scalar "s" value - private static void Drot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) + static void Drot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) { for (var i = 0; i < rowCount; i++) { - var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); - var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); + var z = (c*a.At(i, columnA)) + (s*a.At(i, columnB)); + var tmp = (c*a.At(i, columnB)) - (s*a.At(i, columnA)); a.At(i, columnB, tmp); a.At(i, columnA, z); } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -785,18 +786,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -833,7 +823,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j) * input.At(i, k); + value += U.At(i, j)*input.At(i, k); } value /= S[j]; @@ -847,7 +837,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization double value = 0; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j) * tmp[i]; + value += VT.At(i, j)*tmp[i]; } result.At(j, k, value); @@ -862,17 +852,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -900,7 +880,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j) * input[i]; + value += U.At(i, j)*input[i]; } value /= S[j]; @@ -914,7 +894,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization value = 0; for (int i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j) * tmp[i]; + value += VT.At(i, j)*tmp[i]; } result[j] = value; diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 5ad949dd..336850d2 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -481,7 +481,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Svd Svd(bool computeVectors) { - return new UserSvd(this, computeVectors); + return UserSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Factorization/Svd.cs index b8b7cca5..227c0b36 100644 --- a/src/Numerics/LinearAlgebra/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Factorization/Svd.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 @@ -49,27 +49,66 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// /// Supported data types are double, single, , and . public abstract class Svd : ISolver - where T : struct, IEquatable, IFormattable + where T : struct, IEquatable, IFormattable { + readonly Lazy> _lazyW; + + /// Indicating whether U and VT matrices have been computed during SVD factorization. + protected readonly bool VectorsComputed; + + protected Svd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + { + S = s; + U = u; + VT = vt; + + VectorsComputed = vectorsComputed; + + _lazyW = new Lazy>(ComputeW); + } + + Matrix ComputeW() + { + var rows = U.RowCount; + var columns = VT.ColumnCount; + var result = U.CreateMatrix(rows, columns); + for (var i = 0; i < rows; i++) + { + for (var j = 0; j < columns; j++) + { + if (i == j) + { + result.At(i, i, S[i]); + } + } + } + + return result; + } + /// - /// Gets or sets a value indicating whether to compute U and VT matrices during SVD factorization or not + /// Gets the singular values (Σ) of matrix in ascending value. /// - protected bool ComputeVectors { get; set; } + public Vector S { get; private set; } /// - /// Gets or sets the singular values (Σ) of matrix in ascending value. + /// Gets the left singular vectors (U - m-by-m unitary matrix) /// - public Vector S { get; protected set; } + public Matrix U { get; private set; } /// - /// Gets or sets left singular vectors (U - m-by-m unitary matrix) + /// Gets the transpose right singular vectors (transpose of V, an n-by-n unitary matrix) /// - public Matrix U { get; protected set; } + public Matrix VT { get; private set; } /// - /// Gets or sets transpose right singular vectors (transpose of V, an n-by-n unitary matrix + /// Returns the singular values as a diagonal . /// - public Matrix VT { get; protected set; } + /// The singular values as a diagonal . + public Matrix W + { + get { return _lazyW.Value; } + } /// /// Gets the effective numerical matrix rank. @@ -94,27 +133,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// public abstract T Determinant { get; } - /// Returns the singular values as a diagonal . - /// The singular values as a diagonal . - public Matrix W() - { - var rows = U.RowCount; - var columns = VT.ColumnCount; - var result = U.CreateMatrix(rows, columns); - for (var i = 0; i < rows; i++) - { - for (var j = 0; j < columns; j++) - { - if (i == j) - { - result.At(i, i, S[i]); - } - } - } - - return result; - } - /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -122,13 +140,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -152,13 +164,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index 88d2cb11..42865d30 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -1057,7 +1057,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Svd Svd(bool computeVectors) { - return new DenseSvd(this, computeVectors); + return DenseSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs index 60efd1f9..3fc7330b 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.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 { @@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public sealed class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,19 +57,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// If SVD algorithm failed to converge with matrix . - public DenseSvd(DenseMatrix matrix, bool computeVectors) + public static DenseSvd Create(DenseMatrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); - S = new DenseVector(nm); - U = new DenseMatrix(matrix.RowCount); - VT = new DenseMatrix(matrix.ColumnCount); - Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values); + var s = new DenseVector(nm); + var u = new DenseMatrix(matrix.RowCount); + var vt = new DenseMatrix(matrix.ColumnCount); + Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix) matrix.Clone()).Values, matrix.RowCount, matrix.ColumnCount, s.Values, u.Values, vt.Values); + + return new DenseSvd(s, u, vt, computeVectors); + } + + DenseSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -79,18 +80,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -125,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, input.ColumnCount, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values); } /// @@ -135,17 +125,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -175,7 +155,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector)S).Values, ((DenseMatrix)U).Values, ((DenseMatrix)VT).Values, dinput.Values, 1, dresult.Values); + Control.LinearAlgebraProvider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs index 49dfc58d..ff699323 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Svd.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,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.LinearAlgebra.Factorization; -using MathNet.Numerics.Properties; using System; using System.Linq; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { @@ -51,6 +51,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// public abstract class Svd : Svd { + protected Svd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { + } + /// /// Gets the effective numerical matrix rank. /// diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs index 2743144b..43c96429 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.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 { @@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public sealed class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,20 +57,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Compute the singular U and VT vectors or not. /// If is null. /// - public UserSvd(Matrix matrix, bool computeVectors) + public static UserSvd Create(Matrix matrix, bool computeVectors) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); var matrixCopy = matrix.Clone(); - S = matrixCopy.CreateVector(nm); - U = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); - VT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); + var s = matrixCopy.CreateVector(nm); + var u = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); + var vt = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); const int maxiter = 1000; var e = new float[matrixCopy.ColumnCount]; @@ -94,32 +88,32 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); - S[l] = xnorm; - if (S[l] != 0.0) + s[l] = xnorm; + if (s[l] != 0.0) { if (matrixCopy.At(l, l) != 0.0) { - S[l] = Dsign(S[l], matrixCopy.At(l, l)); + s[l] = Dsign(s[l], matrixCopy.At(l, l)); } - DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f / S[l]); + DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f/s[l]); matrixCopy.At(l, l, (1.0f + matrixCopy.At(l, l))); } - S[l] = -S[l]; + s[l] = -s[l]; } for (j = lp1; j < matrixCopy.ColumnCount; j++) { if (l < nct) { - if (S[l] != 0.0) + if (s[l] != 0.0) { // Apply the transformation. - t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); + t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l)/matrixCopy.At(l, l); for (var ii = l; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t*matrixCopy.At(ii, l))); } } } @@ -129,12 +123,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization e[j] = matrixCopy.At(l, j); } - if (ComputeVectors && l < nct) + if (computeVectors && l < nct) { // Place the transformation in u for subsequent back multiplication. for (i = l; i < matrixCopy.RowCount; i++) { - U.At(i, l, matrixCopy.At(i, l)); + u.At(i, l, matrixCopy.At(i, l)); } } @@ -153,7 +147,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization e[l] = Dsign(e[l], e[lp1]); } - DscalVector(e, lp1, 1.0f / e[l]); + DscalVector(e, lp1, 1.0f/e[l]); e[lp1] = 1.0f + e[lp1]; } @@ -170,26 +164,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - work[ii] += e[j] * matrixCopy.At(ii, j); + work[ii] += e[j]*matrixCopy.At(ii, j); } } for (j = lp1; j < matrixCopy.ColumnCount; j++) { - var ww = -e[j] / e[lp1]; + var ww = -e[j]/e[lp1]; for (var ii = lp1; ii < matrixCopy.RowCount; ii++) { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww*work[ii])); } } } - if (ComputeVectors) + if (computeVectors) { // Place the transformation in v for subsequent back multiplication. for (i = lp1; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, e[i]); + vt.At(i, l, e[i]); } } } @@ -200,12 +194,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var nrtp1 = nrt + 1; if (nct < matrixCopy.ColumnCount) { - S[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); + s[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); } if (matrixCopy.RowCount < m) { - S[m - 1] = 0.0f; + s[m - 1] = 0.0f; } if (nrtp1 < m) @@ -216,52 +210,52 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization e[m - 1] = 0.0f; // If required, generate u. - if (ComputeVectors) + if (computeVectors) { for (j = nctp1 - 1; j < ncu; j++) { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, j, 0.0f); + u.At(i, j, 0.0f); } - U.At(j, j, 1.0f); + u.At(j, j, 1.0f); } for (l = nct - 1; l >= 0; l--) { - if (S[l] != 0.0) + if (s[l] != 0.0) { for (j = l + 1; j < ncu; j++) { - t = -Ddot(U, matrixCopy.RowCount, l, j, l) / U.At(l, l); + t = -Ddot(u, matrixCopy.RowCount, l, j, l)/u.At(l, l); for (var ii = l; ii < matrixCopy.RowCount; ii++) { - U.At(ii, j, U.At(ii, j) + (t * U.At(ii, l))); + u.At(ii, j, u.At(ii, j) + (t*u.At(ii, l))); } } - DscalColumn(U, matrixCopy.RowCount, l, l, -1.0f); - U.At(l, l, 1.0f + U.At(l, l)); + DscalColumn(u, matrixCopy.RowCount, l, l, -1.0f); + u.At(l, l, 1.0f + u.At(l, l)); for (i = 0; i < l; i++) { - U.At(i, l, 0.0f); + u.At(i, l, 0.0f); } } else { for (i = 0; i < matrixCopy.RowCount; i++) { - U.At(i, l, 0.0f); + u.At(i, l, 0.0f); } - U.At(l, l, 1.0f); + u.At(l, l, 1.0f); } } } // If it is required, generate v. - if (ComputeVectors) + if (computeVectors) { for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) { @@ -272,10 +266,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (j = lp1; j < matrixCopy.ColumnCount; j++) { - t = -Ddot(VT, matrixCopy.ColumnCount, l, j, lp1) / VT.At(lp1, l); + t = -Ddot(vt, matrixCopy.ColumnCount, l, j, lp1)/vt.At(lp1, l); for (var ii = l; ii < matrixCopy.ColumnCount; ii++) { - VT.At(ii, j, VT.At(ii, j) + (t * VT.At(ii, l))); + vt.At(ii, j, vt.At(ii, j) + (t*vt.At(ii, l))); } } } @@ -283,10 +277,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (i = 0; i < matrixCopy.ColumnCount; i++) { - VT.At(i, l, 0.0f); + vt.At(i, l, 0.0f); } - VT.At(l, l, 1.0f); + vt.At(l, l, 1.0f); } } @@ -294,19 +288,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (i = 0; i < m; i++) { float r; - if (S[i] != 0.0) + if (s[i] != 0.0) { - t = S[i]; - r = S[i] / t; - S[i] = t; + t = s[i]; + r = s[i]/t; + s[i] = t; if (i < m - 1) { - e[i] = e[i] / r; + e[i] = e[i]/r; } - if (ComputeVectors) + if (computeVectors) { - DscalColumn(U, matrixCopy.RowCount, i, 0, r); + DscalColumn(u, matrixCopy.RowCount, i, 0, r); } } @@ -319,12 +313,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization if (e[i] != 0.0) { t = e[i]; - r = t / e[i]; + r = t/e[i]; e[i] = t; - S[i + 1] = S[i + 1] * r; - if (ComputeVectors) + s[i + 1] = s[i + 1]*r; + if (computeVectors) { - DscalColumn(VT, matrixCopy.ColumnCount, i + 1, 0, r); + DscalColumn(vt, matrixCopy.ColumnCount, i + 1, 0, r); } } } @@ -352,7 +346,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization float test; for (l = m - 2; l >= 0; l--) { - test = Math.Abs(S[l]) + Math.Abs(S[l + 1]); + test = Math.Abs(s[l]) + Math.Abs(s[l + 1]); ztest = test + Math.Abs(e[l]); if (ztest.AlmostEqualInDecimalPlaces(test, 7)) { @@ -382,10 +376,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization test = test + Math.Abs(e[ls - 1]); } - ztest = test + Math.Abs(S[ls]); + ztest = test + Math.Abs(s[ls]); if (ztest.AlmostEqualInDecimalPlaces(test, 7)) { - S[ls] = 0.0f; + s[ls] = 0.0f; break; } } @@ -414,7 +408,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization float cs; switch (kase) { - // Deflate negligible VectorS[m]. + // Deflate negligible VectorS[m]. case 1: f = e[m - 2]; e[m - 2] = 0.0f; @@ -422,72 +416,72 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var kk = l; kk < m - 1; kk++) { k = m - 2 - kk + l; - t1 = S[k]; + t1 = s[k]; Drotg(ref t1, ref f, out cs, out sn); - S[k] = t1; + s[k] = t1; if (k != l) { - f = -sn * e[k - 1]; - e[k - 1] = cs * e[k - 1]; + f = -sn*e[k - 1]; + e[k - 1] = cs*e[k - 1]; } - if (ComputeVectors) + if (computeVectors) { - Drot(VT, matrixCopy.ColumnCount, k, m - 1, cs, sn); + Drot(vt, matrixCopy.ColumnCount, k, m - 1, cs, sn); } } break; - // Split at negligible VectorS[l]. + // Split at negligible VectorS[l]. case 2: f = e[l - 1]; e[l - 1] = 0.0f; for (k = l; k < m; k++) { - t1 = S[k]; + t1 = s[k]; Drotg(ref t1, ref f, out cs, out sn); - S[k] = t1; - f = -sn * e[k]; - e[k] = cs * e[k]; - if (ComputeVectors) + s[k] = t1; + f = -sn*e[k]; + e[k] = cs*e[k]; + if (computeVectors) { - Drot(U, matrixCopy.RowCount, k, l - 1, cs, sn); + Drot(u, matrixCopy.RowCount, k, l - 1, cs, sn); } } break; - // Perform one qr step. + // Perform one qr step. case 3: // Calculate the shift. var scale = 0.0f; - scale = Math.Max(scale, Math.Abs(S[m - 1])); - scale = Math.Max(scale, Math.Abs(S[m - 2])); + scale = Math.Max(scale, Math.Abs(s[m - 1])); + scale = Math.Max(scale, Math.Abs(s[m - 2])); scale = Math.Max(scale, Math.Abs(e[m - 2])); - scale = Math.Max(scale, Math.Abs(S[l])); + scale = Math.Max(scale, Math.Abs(s[l])); scale = Math.Max(scale, Math.Abs(e[l])); - var sm = S[m - 1] / scale; - var smm1 = S[m - 2] / scale; - var emm1 = e[m - 2] / scale; - var sl = S[l] / scale; - var el = e[l] / scale; - var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0f; - var c = (sm * emm1) * (sm * emm1); + var sm = s[m - 1]/scale; + var smm1 = s[m - 2]/scale; + var emm1 = e[m - 2]/scale; + var sl = s[l]/scale; + var el = e[l]/scale; + var b = (((smm1 + sm)*(smm1 - sm)) + (emm1*emm1))/2.0f; + var c = (sm*emm1)*(sm*emm1); var shift = 0.0f; if (b != 0.0 || c != 0.0) { - shift = (float)Math.Sqrt((b * b) + c); + shift = (float) Math.Sqrt((b*b) + c); if (b < 0.0) { shift = -shift; } - shift = c / (b + shift); + shift = c/(b + shift); } - f = ((sl + sm) * (sl - sm)) + shift; - var g = sl * el; + f = ((sl + sm)*(sl - sm)) + shift; + var g = sl*el; // Chase zeros. for (k = l; k < m - 1; k++) @@ -498,24 +492,24 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization e[k - 1] = f; } - f = (cs * S[k]) + (sn * e[k]); - e[k] = (cs * e[k]) - (sn * S[k]); - g = sn * S[k + 1]; - S[k + 1] = cs * S[k + 1]; - if (ComputeVectors) + f = (cs*s[k]) + (sn*e[k]); + e[k] = (cs*e[k]) - (sn*s[k]); + g = sn*s[k + 1]; + s[k + 1] = cs*s[k + 1]; + if (computeVectors) { - Drot(VT, matrixCopy.ColumnCount, k, k + 1, cs, sn); + Drot(vt, matrixCopy.ColumnCount, k, k + 1, cs, sn); } Drotg(ref f, ref g, out cs, out sn); - S[k] = f; - f = (cs * e[k]) + (sn * S[k + 1]); - S[k + 1] = (-sn * e[k]) + (cs * S[k + 1]); - g = sn * e[k + 1]; - e[k + 1] = cs * e[k + 1]; - if (ComputeVectors && k < matrixCopy.RowCount) + s[k] = f; + f = (cs*e[k]) + (sn*s[k + 1]); + s[k + 1] = (-sn*e[k]) + (cs*s[k + 1]); + g = sn*e[k + 1]; + e[k + 1] = cs*e[k + 1]; + if (computeVectors && k < matrixCopy.RowCount) { - Drot(U, matrixCopy.RowCount, k, k + 1, cs, sn); + Drot(u, matrixCopy.RowCount, k, k + 1, cs, sn); } } @@ -523,37 +517,37 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization iter = iter + 1; break; - // Convergence. + // Convergence. case 4: // Make the singular value positive - if (S[l] < 0.0) + if (s[l] < 0.0) { - S[l] = -S[l]; - if (ComputeVectors) + s[l] = -s[l]; + if (computeVectors) { - DscalColumn(VT, matrixCopy.ColumnCount, l, 0, -1.0f); + DscalColumn(vt, matrixCopy.ColumnCount, l, 0, -1.0f); } } // Order the singular value. while (l != mn - 1) { - if (S[l] >= S[l + 1]) + if (s[l] >= s[l + 1]) { break; } - t = S[l]; - S[l] = S[l + 1]; - S[l + 1] = t; - if (ComputeVectors && l < matrixCopy.ColumnCount) + t = s[l]; + s[l] = s[l + 1]; + s[l + 1] = t; + if (computeVectors && l < matrixCopy.ColumnCount) { - Dswap(VT, matrixCopy.ColumnCount, l, l + 1); + Dswap(vt, matrixCopy.ColumnCount, l, l + 1); } - if (ComputeVectors && l < matrixCopy.RowCount) + if (computeVectors && l < matrixCopy.RowCount) { - Dswap(U, matrixCopy.RowCount, l, l + 1); + Dswap(u, matrixCopy.RowCount, l, l + 1); } l = l + 1; @@ -565,9 +559,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } } - if (ComputeVectors) + if (computeVectors) { - VT = VT.Transpose(); + vt = vt.Transpose(); } // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses @@ -579,11 +573,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var tmp = matrixCopy.CreateVector(nm); for (i = 0; i < nm; i++) { - tmp[i] = S[i]; + tmp[i] = s[i]; } - S = tmp; + s = tmp; } + + return new UserSvd(s, u, vt, computeVectors); + } + + UserSvd(Vector s, Matrix u, Matrix vt, bool vectorsComputed) + : base(s, u, vt, vectorsComputed) + { } /// @@ -592,9 +593,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Double value z1 /// Double value z2 /// Result multiplication of signum function and absolute value - private static float Dsign(float z1, float z2) + static float Dsign(float z1, float z2) { - return Math.Abs(z1) * (z2 / Math.Abs(z2)); + return Math.Abs(z1)*(z2/Math.Abs(z2)); } /// @@ -604,7 +605,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The number of rows in /// Column A index to swap /// Column B index to swap - private static void Dswap(Matrix a, int rowCount, int columnA, int columnB) + static void Dswap(Matrix a, int rowCount, int columnA, int columnB) { for (var i = 0; i < rowCount; i++) { @@ -622,11 +623,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Column to scale /// Row to scale from /// Scale value - private static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, float z) + static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, float z) { for (var i = rowStart; i < rowCount; i++) { - a.At(i, column, a.At(i, column) * z); + a.At(i, column, a.At(i, column)*z); } } @@ -636,11 +637,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Source vector /// Row to scale from /// Scale value - private static void DscalVector(float[] a, int start, float z) + static void DscalVector(float[] a, int start, float z) { for (var i = start; i < a.Length; i++) { - a[i] = a[i] * z; + a[i] = a[i]*z; } } @@ -653,7 +654,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Contains the parameter c associated with the Givens rotation /// Contains the parameter s associated with the Givens rotation /// This is equivalent to the DROTG LAPACK routine. - private static void Drotg(ref float da, ref float db, out float c, out float s) + static void Drotg(ref float da, ref float db, out float c, out float s) { float r, z; @@ -675,16 +676,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } else { - var sda = da / scale; - var sdb = db / scale; - r = scale * (float)Math.Sqrt((sda * sda) + (sdb * sdb)); + var sda = da/scale; + var sdb = db/scale; + r = scale*(float) Math.Sqrt((sda*sda) + (sdb*sdb)); if (roe < 0.0) { r = -r; } - c = da / r; - s = db / r; + c = da/r; + s = db/r; z = 1.0f; if (absda > absdb) { @@ -693,7 +694,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization if (absdb >= absda && c != 0.0) { - z = 1.0f / c; + z = 1.0f/c; } } @@ -709,15 +710,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Column index /// Start row index /// Norm2 (Euclidean norm) of the column - private static float Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) + static float Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) { float s = 0; for (var i = rowStart; i < rowCount; i++) { - s += a.At(i, column) * a.At(i, column); + s += a.At(i, column)*a.At(i, column); } - return (float)Math.Sqrt(s); + return (float) Math.Sqrt(s); } /// @@ -726,15 +727,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Source vector /// Start index /// Norm2 (Euclidean norm) of the vector - private static float Dnrm2Vector(float[] a, int rowStart) + static float Dnrm2Vector(float[] a, int rowStart) { float s = 0; for (var i = rowStart; i < a.Length; i++) { - s += a[i] * a[i]; + s += a[i]*a[i]; } - return (float)Math.Sqrt(s); + return (float) Math.Sqrt(s); } /// @@ -746,12 +747,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Index of column B /// Starting row index /// Dot product value - private static float Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) + static float Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) { var z = 0.0f; for (var i = rowStart; i < rowCount; i++) { - z += a.At(i, columnB) * a.At(i, columnA); + z += a.At(i, columnB)*a.At(i, columnA); } return z; @@ -767,17 +768,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Index of column B /// Scalar "c" value /// Scalar "s" value - private static void Drot(Matrix a, int rowCount, int columnA, int columnB, float c, float s) + static void Drot(Matrix a, int rowCount, int columnA, int columnB, float c, float s) { for (var i = 0; i < rowCount; i++) { - var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); - var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); + var z = (c*a.At(i, columnA)) + (s*a.At(i, columnB)); + var tmp = (c*a.At(i, columnB)) - (s*a.At(i, columnA)); a.At(i, columnB, tmp); a.At(i, columnA, z); } } - + /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// @@ -785,18 +786,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -833,7 +823,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j) * input.At(i, k); + value += U.At(i, j)*input.At(i, k); } value /= S[j]; @@ -847,7 +837,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization float value = 0; for (var i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j) * tmp[i]; + value += VT.At(i, j)*tmp[i]; } result.At(j, k, value); @@ -862,17 +852,7 @@ 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"); - } - - if (!ComputeVectors) + if (!VectorsComputed) { throw new InvalidOperationException(Resources.SingularVectorsNotComputed); } @@ -900,7 +880,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { for (var i = 0; i < U.RowCount; i++) { - value += U.At(i, j) * input[i]; + value += U.At(i, j)*input[i]; } value /= S[j]; @@ -914,7 +894,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization value = 0; for (int i = 0; i < VT.ColumnCount; i++) { - value += VT.At(i, j) * tmp[i]; + value += VT.At(i, j)*tmp[i]; } result[j] = value; diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 8ba5eb8e..759605f2 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -481,7 +481,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Svd Svd(bool computeVectors) { - return new UserSvd(this, computeVectors); + return UserSvd.Create(this, computeVectors); } public override Evd Evd() diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/SvdTests.cs index 84e700a6..1405c1fc 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/SvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/SvdTests.cs @@ -24,10 +24,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Complex; -using MathNet.Numerics.LinearAlgebra.Complex.Factorization; using NUnit.Framework; -using System; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization { @@ -38,15 +38,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization /// public class SvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new DenseSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -60,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -97,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserSvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserSvdTests.cs index 7569af0b..29635d10 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserSvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserSvdTests.cs @@ -24,27 +24,19 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization { - using System; using System.Numerics; - using LinearAlgebra.Complex.Factorization; - using NUnit.Framework; /// /// Svd factorization tests for a user matrix. /// public class UserSvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -58,7 +50,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -95,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/SvdTests.cs index d27f432f..114a39dc 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/SvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/SvdTests.cs @@ -24,10 +24,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Complex32; -using MathNet.Numerics.LinearAlgebra.Complex32.Factorization; using NUnit.Framework; -using System; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { @@ -38,15 +38,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization /// public class SvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new DenseSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -60,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -97,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserSvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserSvdTests.cs index dace2e43..2206b465 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserSvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserSvdTests.cs @@ -24,27 +24,19 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { - using System; - using LinearAlgebra.Complex32.Factorization; - using NUnit.Framework; - using Complex32 = Numerics.Complex32; + using Numerics; /// /// Svd factorization tests for a user matrix. /// public class UserSvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -58,7 +50,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -95,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs index 872f1004..ea52f329 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs @@ -24,10 +24,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; -using MathNet.Numerics.LinearAlgebra.Double.Factorization; using NUnit.Framework; -using System; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { @@ -36,15 +36,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization /// public class SvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new DenseSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -58,7 +49,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -95,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs index b7354389..7f8b7adc 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs @@ -24,26 +24,17 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { - using System; - using LinearAlgebra.Double.Factorization; - using NUnit.Framework; - /// /// Svd factorization tests for a user matrix. /// public class UserSvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -57,7 +48,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -94,7 +85,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/SvdTests.cs index 172885fd..2d32f7d0 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/SvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/SvdTests.cs @@ -24,10 +24,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Single; -using MathNet.Numerics.LinearAlgebra.Single.Factorization; using NUnit.Framework; -using System; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization { @@ -36,15 +36,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization /// public class SvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new DenseSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -58,7 +49,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -95,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserSvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserSvdTests.cs index 4c13cf54..3eb3dba9 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserSvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserSvdTests.cs @@ -24,26 +24,17 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; + namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization { - using System; - using LinearAlgebra.Single.Factorization; - using NUnit.Framework; - /// /// Svd factorization tests for a user matrix. /// public class UserSvdTests { - /// - /// Constructor with null throws ArgumentNullException. - /// - [Test] - public void ConstructorNull() - { - Assert.Throws(() => new UserSvd(null, true)); - } - /// /// Can factorize identity matrix. /// @@ -57,7 +48,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var factorSvd = matrixI.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; Assert.AreEqual(matrixI.RowCount, u.RowCount); Assert.AreEqual(matrixI.RowCount, u.ColumnCount); @@ -94,7 +85,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var factorSvd = matrixA.Svd(true); var u = factorSvd.U; var vt = factorSvd.VT; - var w = factorSvd.W(); + var w = factorSvd.W; // Make sure the U has the right dimensions. Assert.AreEqual(row, u.RowCount);