Browse Source

native: added SVD and modified the linear algebra provider as needed

la-knuth
Marcus Cuda 16 years ago
parent
commit
34bc517324
  1. 33
      src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
  2. 141
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  3. 139
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  4. 139
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  5. 142
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  6. 2
      src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.tt
  7. 2
      src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.tt
  8. 2
      src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.tt
  9. 2
      src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.tt
  10. 184
      src/Numerics/Algorithms/LinearAlgebra/native.generic.include
  11. 12
      src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
  12. 330
      src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

33
src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs

@ -432,9 +432,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.
/// <param name="work">The work array. On exit, work[0] contains the optimal work size value.
/// </param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
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
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value. </param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
void SvdSolve(T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x);
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value. </param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.
/// </param>
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);
/// <summary>
/// Solves A*X=B for X using a previously SVD decomposed matrix.
@ -479,7 +456,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="s">The s values returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[])"/>.</param>
/// <param name="u">The left singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int, T[],T[],T[])"/>.</param>
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[],T[])"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="b">The B matrix</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
void SvdSolveFactored(int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x);

141
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.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
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
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
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);
}
/// <summary>

139
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.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
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
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
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);
}

139
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.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
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
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
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);
}

142
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@ -1838,7 +1838,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
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.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
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
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
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);
}

2
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" #>

2
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" #>

2
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" #>

2
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" #>

184
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");
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">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.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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");
}
/// <summary>
/// Solves A*X=B for X using a previously SVD decomposed matrix.
/// </summary>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The s values returned by <see cref="SingularValueDecomposition(bool,<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[],<#=dataType#>[])"/>.</param>
/// <param name="u">The left singular vectors returned by <see cref="SingularValueDecomposition(bool,<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[],<#=dataType#>[])"/>.</param>
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[],<#=dataType#>[])"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[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);
}

12
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

330
src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

@ -1032,6 +1032,336 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqual(test[1, 1], x[3], 14);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix.
/// </summary>
[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);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix.
/// </summary>
[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);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix.
/// </summary>
[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);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix using
/// a work array.
/// </summary>
[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);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix using
/// a work array.
/// </summary>
[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);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix using
/// a work array.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a tall A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix
/// using a factored matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a tall A matrix
/// using a factored matrix.
/// </summary>
[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);
}
/// <summary>
/// Checks to see if a matrix and array contain the same values.
/// </summary>

Loading…
Cancel
Save