diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
index 5b62552e..8d4fdb85 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
@@ -432,9 +432,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
+ /// The work array. On exit, work[0] contains the optimal work size value.
///
/// This is equivalent to the GESVD LAPACK routine.
void SingularValueDecomposition(bool computeVectors, T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] work);
@@ -442,34 +440,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// On entry the B matrix; on exit the X matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- void SvdSolve(T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x);
-
- ///
- /// Solves A*X=B for X using the singular value decomposition of A.
- ///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// On entry the B matrix; on exit the X matrix.
+ /// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- ///
- void SvdSolve(T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x, T[] work);
+ void SvdSolve(T[] a, int rowsA, int columnsA, T[] b, int columnsB, T[] x);
///
/// Solves A*X=B for X using a previously SVD decomposed matrix.
@@ -479,7 +456,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// The s values returned by .
/// The left singular vectors returned by .
/// The right singular vectors returned by .
- /// On entry the B matrix; on exit the X matrix.
+ /// The B matrix
/// The number of columns of B.
/// On exit, the solution matrix.
void SvdSolveFactored(int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x);
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
index bf600e9a..e148e56e 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
@@ -1877,8 +1877,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
- // Actually "work = new Complex[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)];
+ var work = new Complex[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
@@ -1894,9 +1893,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
+ /// The work array. Length should be at least .
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
{
@@ -2571,114 +2568,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- public virtual void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x)
- {
- if (a == null)
- {
- throw new ArgumentNullException("a");
- }
-
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
- if (b == null)
- {
- throw new ArgumentNullException("b");
- }
-
- if (x == null)
- {
- throw new ArgumentNullException("x");
- }
-
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
- if (b.Length != rowsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- if (x.Length != columnsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- // Actually "work = new Complex[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)];
- SvdSolve(a, rowsA, columnsA, s, u, vt, b, columnsB, x, work);
- }
-
- ///
- /// Solves A*X=B for X using the singular value decomposition of A.
- ///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- public virtual void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x, Complex[] work)
+ public virtual void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
if (b == null)
{
throw new ArgumentNullException("b");
@@ -2689,21 +2591,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("x");
}
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
@@ -2714,19 +2601,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
- if (work.Length == 0)
- {
- throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
- }
-
- if (work.Length < rowsA)
- {
- work[0] = rowsA;
- throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
- }
+ var work = new Complex[rowsA];
+ var s = new Complex[Math.Min(rowsA, columnsA)];
+ var u = new Complex[rowsA * rowsA];
+ var vt = new Complex[columnsA * columnsA];
- SingularValueDecomposition(true, a, rowsA, columnsA, s, u, vt, work);
- SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
+ var clone = new Complex[a.Length];
+ Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfComplex);
+ SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
+ SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
///
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
index 406b6bdf..fcfedef1 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
@@ -1877,8 +1877,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
- // Actually "work = new Complex32[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)];
+ var work = new Complex32[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
@@ -1894,9 +1893,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
+ /// The work array. Length should be at least .
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work)
{
@@ -2571,114 +2568,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x)
- {
- if (a == null)
- {
- throw new ArgumentNullException("a");
- }
-
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
- if (b == null)
- {
- throw new ArgumentNullException("b");
- }
-
- if (x == null)
- {
- throw new ArgumentNullException("x");
- }
-
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
- if (b.Length != rowsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- if (x.Length != columnsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- // TODO: Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)];
- SvdSolve(a, rowsA, columnsA, s, u, vt, b, columnsB, x, work);
- }
-
- ///
- /// Solves A*X=B for X using the singular value decomposition of A.
- ///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work)
+ public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
if (b == null)
{
throw new ArgumentNullException("b");
@@ -2689,21 +2591,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("x");
}
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
@@ -2714,18 +2601,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
- if (work.Length == 0)
- {
- throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
- }
-
- if (work.Length < rowsA)
- {
- work[0] = rowsA;
- throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
- }
+ var work = new Complex32[rowsA];
+ var s = new Complex32[Math.Min(rowsA, columnsA)];
+ var u = new Complex32[rowsA * rowsA];
+ var vt = new Complex32[columnsA * columnsA];
- SingularValueDecomposition(true, a, rowsA, columnsA, s, u, vt, work);
+ var clone = new Complex32[a.Length];
+ Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfComplex32);
+ SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
index b12c56b9..46cf3324 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
@@ -1873,8 +1873,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
- // Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))];
+ var work = new double[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
@@ -1890,9 +1889,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
+ /// The work array. Length should be at least .
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work)
{
@@ -2628,114 +2625,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- public virtual void SvdSolve(double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x)
- {
- if (a == null)
- {
- throw new ArgumentNullException("a");
- }
-
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
- if (b == null)
- {
- throw new ArgumentNullException("b");
- }
-
- if (x == null)
- {
- throw new ArgumentNullException("x");
- }
-
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
- if (b.Length != rowsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- if (x.Length != columnsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- // Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))];
- SvdSolve(a, rowsA, columnsA, s, u, vt, b, columnsB, x, work);
- }
-
- ///
- /// Solves A*X=B for X using the singular value decomposition of A.
- ///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- public virtual void SvdSolve(double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x, double[] work)
+ public virtual void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
if (b == null)
{
throw new ArgumentNullException("b");
@@ -2746,21 +2648,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("x");
}
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
@@ -2771,18 +2658,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
- if (work.Length == 0)
- {
- throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
- }
-
- if (work.Length < rowsA)
- {
- work[0] = rowsA;
- throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
- }
+ var work = new double[rowsA];
+ var s = new double[Math.Min(rowsA, columnsA)];
+ var u = new double[rowsA * rowsA];
+ var vt = new double[columnsA * columnsA];
- SingularValueDecomposition(true, a, rowsA, columnsA, s, u, vt, work);
+ var clone = new double[a.Length];
+ Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble);
+ SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
index e1c626a0..7921cbb9 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
@@ -1838,7 +1838,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt)
{
- if (a == null)
+ if (a == null)
{
throw new ArgumentNullException("a");
}
@@ -1873,8 +1873,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
- // Actually "work = new float[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new float[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))];
+ var work = new float[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
@@ -1890,10 +1889,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- /// This is equivalent to the GESVD LAPACK routine.
+ /// The work array. Length should be at least .
public virtual void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work)
{
if (a == null)
@@ -2630,114 +2626,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- public virtual void SvdSolve(float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x)
- {
- if (a == null)
- {
- throw new ArgumentNullException("a");
- }
-
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
- if (b == null)
- {
- throw new ArgumentNullException("b");
- }
-
- if (x == null)
- {
- throw new ArgumentNullException("x");
- }
-
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
- if (b.Length != rowsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- if (x.Length != columnsA * columnsB)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
- }
-
- // Actually "work = new float[aRows]" is acceptable size of work array. I set size proposed in method description
- var work = new float[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))];
- SvdSolve(a, rowsA, columnsA, s, u, vt, b, columnsB, x, work);
- }
-
- ///
- /// Solves A*X=B for X using the singular value decomposition of A.
- ///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- public virtual void SvdSolve(float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x, float[] work)
+ public virtual void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
- if (s == null)
- {
- throw new ArgumentNullException("s");
- }
-
- if (u == null)
- {
- throw new ArgumentNullException("u");
- }
-
- if (vt == null)
- {
- throw new ArgumentNullException("vt");
- }
-
if (b == null)
{
throw new ArgumentNullException("b");
@@ -2748,21 +2649,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("x");
}
- if (u.Length != rowsA * rowsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
- }
-
- if (vt.Length != columnsA * columnsA)
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
- }
-
- if (s.Length != Math.Min(rowsA, columnsA))
- {
- throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
- }
-
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
@@ -2773,18 +2659,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
- if (work.Length == 0)
- {
- throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
- }
-
- if (work.Length < rowsA)
- {
- work[0] = rowsA;
- throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
- }
+ var work = new float[rowsA];
+ var s = new float[Math.Min(rowsA, columnsA)];
+ var u = new float[rowsA * rowsA];
+ var vt = new float[columnsA * columnsA];
- SingularValueDecomposition(true, a, rowsA, columnsA, s, u, vt, work);
+ var clone = new float[a.Length];
+ Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat);
+ SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.tt
index 86a5d533..c87fb0a4 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.tt
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.tt
@@ -7,6 +7,8 @@
<# string one = "Complex.One";#>
<# string prefix = "z";#>
<# string reff = "ref ";#>
+<# string svd_work = "2 * Math.Min(rowsA, columnsA) + Math.Max(rowsA, columnsA)";#>
+<# string data_size = "Constants.SizeOfComplex";#>
<#@ include file="..\native.header.include" #>
<#@ include file="..\native.generic.include" #>
<#@ include file="..\native.vector.include" #>
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.tt
index b956a489..3e99ed8a 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.tt
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.tt
@@ -7,6 +7,8 @@
<# string one = "Complex32.One";#>
<# string prefix = "c";#>
<# string reff = "ref ";#>
+<# string svd_work = "2 * Math.Min(rowsA, columnsA) + Math.Max(rowsA, columnsA)";#>
+<# string data_size = "Constants.SizeOfComplex32";#>
<#@ include file="..\native.header.include" #>
<#@ include file="..\native.generic.include" #>
<#@ include file="..\native.vector.include" #>
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.tt
index 48c1e051..ed72de26 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.tt
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.tt
@@ -7,6 +7,8 @@
<# string one = "1.0";#>
<# string prefix = "d";#>
<# string reff = "";#>
+<# string svd_work = "Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))";#>
+<# string data_size = "Constants.SizeOfDouble";#>
<#@ include file="..\native.header.include" #>
<#@ include file="..\native.generic.include" #>
<#@ include file="..\native.vector.include" #>
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.tt
index 6b003cc9..82cfba7f 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.tt
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.tt
@@ -7,6 +7,8 @@
<# string one = "1.0f";#>
<# string prefix = "s";#>
<# string reff = "";#>
+<# string svd_work = "Math.Max((3 * Math.Min(rowsA, columnsA) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))";#>
+<# string data_size = "Constants.SizeOfFloat";#>
<#@ include file="..\native.header.include" #>
<#@ include file="..\native.generic.include" #>
<#@ include file="..\native.vector.include" #>
diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include
index ea4f4f65..4839eb7f 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include
+++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include
@@ -883,83 +883,161 @@
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt)
{
- throw new NotImplementedException();
- }
+ if (a == null)
+ {
+ throw new ArgumentNullException("a");
+ }
- ///
- /// Computes the singular value decomposition of A.
- ///
- /// Compute the singular U and VT vectors or not.
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// If is true, on exit U contains the left
- /// singular vectors.
- /// If is true, on exit VT contains the transposed
- /// right singular vectors.
- /// The work array. For real matrices, the work array should be at least
- /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
- /// On exit, work[0] contains the optimal work size value.
- /// This is equivalent to the GESVD LAPACK routine.
- [SecuritySafeCritical]
- public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] work)
- {
- throw new NotImplementedException();
+ if (s == null)
+ {
+ throw new ArgumentNullException("s");
+ }
+
+ if (u == null)
+ {
+ throw new ArgumentNullException("u");
+ }
+
+ if (vt == null)
+ {
+ throw new ArgumentNullException("vt");
+ }
+
+ if (u.Length != rowsA * rowsA)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
+ }
+
+ if (vt.Length != columnsA * columnsA)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
+ }
+
+ if (s.Length != Math.Min(rowsA, columnsA))
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
+ }
+
+ var work = new <#=dataType#>[<#=svd_work#>];
+ SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
- /// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
+ /// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
- /// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
- [SecuritySafeCritical]
- public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
+ public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
{
- throw new NotImplementedException();
+ if (a == null)
+ {
+ throw new ArgumentNullException("a");
+ }
+
+ if (b == null)
+ {
+ throw new ArgumentNullException("b");
+ }
+
+ if (x == null)
+ {
+ throw new ArgumentNullException("x");
+ }
+
+ if (b.Length != rowsA * columnsB)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
+ }
+
+ if (x.Length != columnsA * columnsB)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
+ }
+
+ var work = new <#=dataType#>[<#=svd_work#>];
+ var s = new <#=dataType#>[Math.Min(rowsA, columnsA)];
+ var u = new <#=dataType#>[rowsA * rowsA];
+ var vt = new <#=dataType#>[columnsA * columnsA];
+
+ var clone = new <#=dataType#>[a.Length];
+ Buffer.BlockCopy(a, 0, clone, 0, a.Length * <#=data_size#>);
+ SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
+ SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
///
- /// Solves A*X=B for X using the singular value decomposition of A.
+ /// Computes the singular value decomposition of A.
///
+ /// Compute the singular U and VT vectors or not.
/// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The singular values of A in ascending value.
- /// On exit U contains the left singular vectors.
- /// On exit VT contains the transposed right singular vectors.
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
+ /// If is true, on exit U contains the left
+ /// singular vectors.
+ /// If is true, on exit VT contains the transposed
+ /// right singular vectors.
/// The work array. For real matrices, the work array should be at least
/// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N).
/// On exit, work[0] contains the optimal work size value.
+ /// This is equivalent to the GESVD LAPACK routine.
[SecuritySafeCritical]
- public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work)
+ public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] work)
{
- throw new NotImplementedException();
- }
+ if (a == null)
+ {
+ throw new ArgumentNullException("a");
+ }
- ///
- /// Solves A*X=B for X using a previously SVD decomposed matrix.
- ///
- /// The number of rows in the A matrix.
- /// The number of columns in the A matrix.
- /// The s values returned by .
- /// The left singular vectors returned by .
- /// The right singular vectors returned by .
- /// The B matrix.
- /// The number of columns of B.
- /// On exit, the solution matrix.
- [SecuritySafeCritical]
- public override void SvdSolveFactored(int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
- {
- throw new NotImplementedException();
+ if (s == null)
+ {
+ throw new ArgumentNullException("s");
+ }
+
+ if (u == null)
+ {
+ throw new ArgumentNullException("u");
+ }
+
+ if (vt == null)
+ {
+ throw new ArgumentNullException("vt");
+ }
+
+ if (work == null)
+ {
+ throw new ArgumentNullException("work");
+ }
+
+ if (u.Length != rowsA * rowsA)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
+ }
+
+ if (vt.Length != columnsA * columnsA)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
+ }
+
+ if (s.Length != Math.Min(rowsA, columnsA))
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
+ }
+
+ if (work.Length == 0)
+ {
+ throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
+ }
+
+ if (work.Length < <#=svd_work#>)
+ {
+ work[0] = <#=svd_work#>;
+ throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
+ }
+
+ SafeNativeMethods.<#=prefix#>_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length);
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
index 8e2695f3..bcbef2c3 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
+++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
@@ -246,4 +246,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len);
+ [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+ internal static extern int s_svd_factor(bool compute_vectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len);
+
+ [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+ internal static extern int d_svd_factor(bool compute_vectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len);
+
+ [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+ internal static extern int c_svd_factor(bool compute_vectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len);
+
+ [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+ internal static extern int z_svd_factor(bool compute_vectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len);
+
#endregion LAPACK
diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
index e1a82d1d..ac0b42fc 100644
--- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
+++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
@@ -1032,6 +1032,336 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqual(test[1, 1], x[3], 14);
}
+ ///
+ /// Can compute the SVD factorization of a square matrix.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfSquareMatrix()
+ {
+ var matrix = _matrices["Square3x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.RowCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 0], result[2, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 1], result[2, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 2], result[0, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 2], result[1, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 2], result[2, 2], 14);
+ }
+
+ ///
+ /// Can compute the SVD factorization of a tall matrix.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfTallMatrix()
+ {
+ var matrix = _matrices["Tall3x2"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.ColumnCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 0], result[2, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 1], result[2, 1], 14);
+ }
+
+ ///
+ /// Can compute the SVD factorization of a wide matrix.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfWideMatrix()
+ {
+ var matrix = _matrices["Wide2x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.RowCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 2], result[0, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 2], result[1, 2], 14);
+ }
+
+ ///
+ /// Can compute the SVD factorization of a square matrix using
+ /// a work array.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfSquareMatrixWithWorkArray()
+ {
+ var matrix = _matrices["Square3x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.RowCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+ var work = new double[100];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 0], result[2, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 1], result[2, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 2], result[0, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 2], result[1, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 2], result[2, 2], 14);
+ }
+
+ ///
+ /// Can compute the SVD factorization of a tall matrix using
+ /// a work array.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfTallMatrixWithWorkArray()
+ {
+ var matrix = _matrices["Tall3x2"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.ColumnCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+ var work = new double[100];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 0], result[2, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[2, 1], result[2, 1], 14);
+ }
+
+ ///
+ /// Can compute the SVD factorization of a wide matrix using
+ /// a work array.
+ ///
+ [Test]
+ public void CanComputeSVDFactorizationOfWideMatrixWithWorkArray()
+ {
+ var matrix = _matrices["Wide2x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.RowCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+ var work = new double[100];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
+
+ var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
+ for (var index = 0; index < s.Length; index++)
+ {
+ w[index, index] = s[index];
+ }
+
+ var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
+ var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
+ var result = mU * w * mV;
+
+ AssertHelpers.AlmostEqual(matrix[0, 0], result[0, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 0], result[1, 0], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 1], result[0, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 1], result[1, 1], 14);
+ AssertHelpers.AlmostEqual(matrix[0, 2], result[0, 2], 14);
+ AssertHelpers.AlmostEqual(matrix[1, 2], result[1, 2], 14);
+ }
+
+ ///
+ /// Can solve Ax=b using SVD factorization with a square A matrix.
+ ///
+ [Test]
+ public void CanSolveUsingSVDSquareMatrix()
+ {
+ var matrix = _matrices["Square3x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+ var x = new double[matrix.ColumnCount * 2];
+ Provider.SvdSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x);
+
+ NotModified(3, 3, a, matrix);
+
+ var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
+ var mb = matrix * mx;
+
+ AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14);
+ AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14);
+ AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14);
+ AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14);
+ AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14);
+ AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14);
+ }
+
+ ///
+ /// Can solve Ax=b using SVD factorization with a tall A matrix.
+ ///
+ [Test]
+ public void CanSolveUsingSVDTallMatrix()
+ {
+ var matrix = _matrices["Tall3x2"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+ var x = new double[matrix.ColumnCount * 2];
+ Provider.SvdSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x);
+
+ NotModified(3, 2, a, matrix);
+
+ var mb = new DenseMatrix(matrix.RowCount, 2, b);
+ var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb;
+
+ AssertHelpers.AlmostEqual(test[0, 0], x[0], 14);
+ AssertHelpers.AlmostEqual(test[1, 0], x[1], 14);
+ AssertHelpers.AlmostEqual(test[0, 1], x[2], 14);
+ AssertHelpers.AlmostEqual(test[1, 1], x[3], 14);
+ }
+
+ ///
+ /// Can solve Ax=b using SVD factorization with a square A matrix
+ /// using a factored matrix.
+ ///
+ [Test]
+ public void CanSolveUsingSVDSquareMatrixOnFactoredMatrix()
+ {
+ var matrix = _matrices["Square3x3"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.RowCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);
+
+ var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+ var x = new double[matrix.ColumnCount * 2];
+ Provider.SvdSolveFactored(matrix.RowCount, matrix.ColumnCount, s, u, vt, b, 2, x);
+
+ var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
+ var mb = matrix * mx;
+
+ AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14);
+ AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14);
+ AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14);
+ AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14);
+ AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14);
+ AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14);
+ }
+
+ ///
+ /// Can solve Ax=b using SVD factorization with a tall A matrix
+ /// using a factored matrix.
+ ///
+ [Test]
+ public void CanSolveUsingSVDTallMatrixOnFactoredMatrix()
+ {
+ var matrix = _matrices["Tall3x2"];
+ var a = new double[matrix.RowCount * matrix.ColumnCount];
+ Array.Copy(matrix.Data, a, a.Length);
+
+ var s = new double[matrix.ColumnCount];
+ var u = new double[matrix.RowCount * matrix.RowCount];
+ var vt = new double[matrix.ColumnCount * matrix.ColumnCount];
+
+ Provider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt);
+
+ var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+ var x = new double[matrix.ColumnCount * 2];
+ Provider.SvdSolveFactored(matrix.RowCount, matrix.ColumnCount, s, u, vt, b, 2, x);
+
+ var mb = new DenseMatrix(matrix.RowCount, 2, b);
+ var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb;
+
+ AssertHelpers.AlmostEqual(test[0, 0], x[0], 14);
+ AssertHelpers.AlmostEqual(test[1, 0], x[1], 14);
+ AssertHelpers.AlmostEqual(test[0, 1], x[2], 14);
+ AssertHelpers.AlmostEqual(test[1, 1], x[3], 14);
+ }
+
///
/// Checks to see if a matrix and array contain the same values.
///