Browse Source

Small additions to LU and initial unit tests.

la-knuth
Jurgen Van Gael 17 years ago
committed by Christoph Ruegg
parent
commit
9755692cf8
  1. 62
      src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs
  2. 14
      src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
  3. 111
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
  4. 62
      src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs
  5. 57
      src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
  6. 12
      src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs
  7. 21
      src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
  8. 3
      src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs
  9. 332
      src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs
  10. 1
      src/UnitTests/UnitTests.csproj

62
src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs

@ -24,7 +24,7 @@
/* This file is automatically generated - do not modify it.
Change NativeLinearAlgebraProvider.include instead.
Last generated on: 16/05/2010 10:39:25
Last generated on: 18/05/2010 22:54:23
*/
namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
{
@ -319,16 +319,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
SafeNativeMethods.d_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(double[] a, int[] ipiv)
public void LUFactor(double[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -919,16 +920,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(float[] a, int[] ipiv)
public void LUFactor(float[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -1519,16 +1521,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
SafeNativeMethods.z_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex[] a, int[] ipiv)
public void LUFactor(Complex[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -2119,16 +2122,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
SafeNativeMethods.c_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex32[] a, int[] ipiv)
public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}

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

@ -208,17 +208,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="c">The c matrix.</param>
void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, T alpha, T[] a,
int aRows, int aColumns, T[] b, int bRows, int bColumns, T beta, T[] c);
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An <paramref name="aOrder"/> by <paramref name="aOrder"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="A"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="A"/>.</param>
/// <param name="aOrder">The order of the square matrix <paramref name="A"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="aOrder"/>.</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
void LUFactor(T[] a, int aOrder, int[] ipiv);
void LUFactor(T[] data, int order, int[] ipiv);
/// <summary>
/// Computes the inverse of matrix using LU factorization.

111
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@ -691,15 +691,79 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An <paramref name="aOrder"/> by <paramref name="aOrder"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="A"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="A"/>.</param>
/// <param name="aOrder">The order of the square matrix <paramref name="A"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="aOrder"/>.</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(double[] a, int aOrder, int[] ipiv)
public void LUFactor(double[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
// Initialize the pivot matrix to the identity permutation.
for (int i = 0; i < order; i++)
{
ipiv[i] = i;
}
double[] LUcolj = new double[order];
// Outer loop.
for (int j = 0; j < order; j++)
{
int indexj = j * order;
int indexjj = indexj + j;
// Make a copy of the j-th column to localize references.
for (int i = 0; i < order; i++)
{
LUcolj[i] = data[indexj + i];
}
// Apply previous transformations.
for (int i = 0; i < order; i++)
{
// Most of the time is spent in the following dot product.
int kmax = System.Math.Min(i, j);
double s = 0.0;
for (int k = 0; k < kmax; k++)
{
s += data[k * order + i] * LUcolj[k];
}
data[indexj + i] = LUcolj[i] -= s;
}
// Find pivot and exchange if necessary.
int p = j;
for (int i = j + 1; i < order; i++)
{
if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p]))
{
p = i;
}
}
if (p != j)
{
for (int k = 0; k < order; k++)
{
int indexk = k * order;
int indexkp = indexk + p;
int indexkj = indexk + j;
double temp = data[indexkp];
data[indexkp] = data[indexkj];
data[indexkj] = temp;
}
ipiv[j] = p;
}
// Compute multipliers.
if (j < order & data[indexjj] != 0.0)
{
for (int i = j + 1; i < order; i++)
{
data[indexj + i] /= data[indexjj];
}
}
}
}
public void LUInverse(double[] a)
@ -1562,7 +1626,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
public void LUFactor(float[] a, int[] ipiv)
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(float[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -2379,7 +2452,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
public void LUFactor(Complex[] a, int[] ipiv)
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -3161,7 +3243,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
}
public void LUFactor(Complex32[] a, int[] ipiv)
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}

62
src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs

@ -24,7 +24,7 @@
/* This file is automatically generated - do not modify it.
Change NativeLinearAlgebraProvider.include instead.
Last generated on: 16/05/2010 10:39:49
Last generated on: 18/05/2010 22:54:25
*/
namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
{
@ -318,16 +318,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
SafeNativeMethods.d_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(double[] a, int[] ipiv)
public void LUFactor(double[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -918,16 +919,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(float[] a, int[] ipiv)
public void LUFactor(float[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -1518,16 +1520,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
SafeNativeMethods.z_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex[] a, int[] ipiv)
public void LUFactor(Complex[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -2118,16 +2121,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
SafeNativeMethods.c_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex32[] a, int[] ipiv)
public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}

57
src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include

@ -337,13 +337,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An <paramref name="aOrder"/> by <paramref name="aOrder"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="A"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="A"/>.</param>
/// <param name="aOrder">The order of the square matrix <paramref name="A"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="aOrder"/>.</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(double[] a, int aOrder, int[] ipiv)
public void LUFactor(double[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -946,16 +946,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(float[] a, int[] ipiv)
public void LUFactor(float[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -1558,16 +1559,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
SafeNativeMethods.z_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex[] a, int[] ipiv)
public void LUFactor(Complex[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
@ -2170,16 +2172,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
SafeNativeMethods.c_matrix_multiply(transposeA, transposeB, m, n, k, ref alpha, a, b, ref beta, c);
}
/// <summary>
/// Computes the LU factorization of A.
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="a">An m by n matrix. The matrix is overwritten with the
/// the LU factorization On exit.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size
/// of the array must be min(m,n).</param>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
public void LUFactor(Complex32[] a, int[] ipiv)
public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}

12
src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs

@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
public static class ExtensionMethods
{
/// <summary>
/// Computes the cholesky decomposition for a matrix.
/// Computes the Cholesky decomposition for a matrix.
/// </summary>
/// <param name="matrix">The matrix to factor.</param>
/// <returns>The Cholesky decomposition object.</returns>
@ -47,5 +47,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{
return Factorization.Cholesky.Create(matrix);
}
/// <summary>
/// Computes the LU decomposition for a matrix.
/// </summary>
/// <param name="matrix">The matrix to factor.</param>
/// <returns>The LU decomposition object.</returns>
public static LU LU(this Matrix matrix)
{
return Factorization.LU.Create(matrix);
}
}
}

21
src/Numerics/LinearAlgebra/Double/Factorization/LU.cs

@ -167,5 +167,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
/// <param name="input">The right hand side vector, <b>b</b>.</param>
/// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
public abstract void Solve(Vector input, Vector result);
/// <summary>
/// Pivot a matrix according to this LU decomposition.
/// </summary>
/// <param name="data">The matrix to pivot.</param>
public void Pivot(Matrix data)
{
for (int i = 0; i < mPivots.Length; i++)
{
if (mPivots[i] != i)
{
int p = mPivots[i];
for (int j = 0; j < data.ColumnCount; j++)
{
double temp = data.At(p, j);
data.At(p, j, data.At(i, j));
data.At(i, j, temp);
}
}
}
}
}
}

3
src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs

@ -46,6 +46,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
var I = DenseMatrix.Identity(order);
var C = I.Cholesky();
Assert.AreEqual(I.RowCount, C.Factor.RowCount);
Assert.AreEqual(I.ColumnCount, C.Factor.ColumnCount);
for (var i = 0; i < C.Factor.RowCount; i++)
{
for (var j = 0; j < C.Factor.ColumnCount; j++)

332
src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs

@ -0,0 +1,332 @@
// <copyright file="LUTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2010 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
{
using System.Collections.Generic;
using MbUnit.Framework;
using LinearAlgebra.Double;
using LinearAlgebra.Double.Factorization;
public class LUTests
{
[Test]
[Row(1)]
[Row(10)]
[Row(100)]
public void CanFactorizeIdentity(int order)
{
var I = DenseMatrix.Identity(order);
var lu = I.LU();
// Check lower triangular part.
var L = lu.L;
Assert.AreEqual(I.RowCount, L.RowCount);
Assert.AreEqual(I.ColumnCount, L.ColumnCount);
for (var i = 0; i < L.RowCount; i++)
{
for (var j = 0; j < L.ColumnCount; j++)
{
if (i == j)
{
Assert.AreEqual(1.0, L[i, j]);
}
else
{
Assert.AreEqual(0.0, L[i, j]);
}
}
}
// Check upper triangular part.
var U = lu.U;
Assert.AreEqual(I.RowCount, U.RowCount);
Assert.AreEqual(I.ColumnCount, U.ColumnCount);
for (var i = 0; i < U.RowCount; i++)
{
for (var j = 0; j < U.ColumnCount; j++)
{
if (i == j)
{
Assert.AreEqual(1.0, U[i, j]);
}
else
{
Assert.AreEqual(0.0, U[i, j]);
}
}
}
}
[Test]
[Row(3,5)]
[Row(5,3)]
[ExpectedArgumentException]
public void LUFailsWithNonSquareMatrix(int row, int col)
{
var I = new DenseMatrix(row, col);
var lu = I.LU();
}
[Test]
[Row(1)]
[Row(10)]
[Row(100)]
public void IdentityDeterminantIsOne(int order)
{
var I = DenseMatrix.Identity(order);
var lu = I.LU();
Assert.AreEqual(1.0, lu.Determinant);
}
[Test]
[Row(1)]
[Row(2)]
[Row(5)]
[Row(10)]
[Row(50)]
[Row(100)]
[MultipleAsserts]
public void CanFactorizeRandomMatrix(int order)
{
var X = MatrixLoader.GenerateRandomMatrix(order, order);
var lu = X.LU();
var L = lu.L;
var U = lu.U;
// Make sure the factors have the right dimensions.
Assert.AreEqual(order, L.RowCount);
Assert.AreEqual(order, L.ColumnCount);
Assert.AreEqual(order, U.RowCount);
Assert.AreEqual(order, U.ColumnCount);
// Make sure the L factor is lower triangular.
for (int i = 0; i < L.RowCount; i++)
{
Assert.AreEqual(1.0, L[i, i]);
for (int j = i+1; j < L.ColumnCount; j++)
{
Assert.AreEqual(0.0, L[i, j]);
}
}
// Make sure the U factor is upper triangular.
for (int i = 0; i < L.RowCount; i++)
{
for (int j = 0; j < i; j++)
{
Assert.AreEqual(0.0, U[i, j]);
}
}
// Make sure the cholesky factor times it's transpose is the original matrix.
var XfromLU = L * U;
lu.Pivot(XfromLU);
for (int i = 0; i < XfromLU.RowCount; i++)
{
for (int j = 0; j < XfromLU.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(X[i, j], XfromLU[i, j], 1.0e-11);
}
}
}
/*[Test]
[Row(1)]
[Row(2)]
[Row(5)]
[Row(10)]
[Row(50)]
[Row(100)]
[MultipleAsserts]
public void CanSolveForRandomVector(int order)
{
var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
var ACopy = A.Clone();
var chol = A.Cholesky();
var b = MatrixLoader.GenerateRandomVector(order);
var x = chol.Solve(b);
Assert.AreEqual(b.Count, x.Count);
var bReconstruct = A * x;
// Check the reconstruction.
for (int i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
}
// Make sure A didn't change.
for (int i = 0; i < A.RowCount; i++)
{
for (int j = 0; j < A.ColumnCount; j++)
{
Assert.AreEqual(ACopy[i, j], A[i, j]);
}
}
}
[Test]
[Row(1,1)]
[Row(2,4)]
[Row(5,8)]
[Row(10,3)]
[Row(50,10)]
[Row(100,100)]
[MultipleAsserts]
public void CanSolveForRandomMatrix(int row, int col)
{
var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row);
var ACopy = A.Clone();
var chol = A.Cholesky();
var B = MatrixLoader.GenerateRandomMatrix(row, col);
var X = chol.Solve(B);
Assert.AreEqual(B.RowCount, X.RowCount);
Assert.AreEqual(B.ColumnCount, X.ColumnCount);
var BReconstruct = A * X;
// Check the reconstruction.
for (int i = 0; i < B.RowCount; i++)
{
for (int j = 0; j < B.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
}
}
// Make sure A didn't change.
for (int i = 0; i < A.RowCount; i++)
{
for (int j = 0; j < A.ColumnCount; j++)
{
Assert.AreEqual(ACopy[i, j], A[i, j]);
}
}
}
[Test]
[Row(1)]
[Row(2)]
[Row(5)]
[Row(10)]
[Row(50)]
[Row(100)]
[MultipleAsserts]
public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
{
var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
var ACopy = A.Clone();
var chol = A.Cholesky();
var b = MatrixLoader.GenerateRandomVector(order);
var bCopy = b.Clone();
var x = new DenseVector(order);
chol.Solve(b, x);
Assert.AreEqual(b.Count, x.Count);
var bReconstruct = A * x;
// Check the reconstruction.
for (int i = 0; i < order; i++)
{
Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
}
// Make sure A didn't change.
for (int i = 0; i < A.RowCount; i++)
{
for (int j = 0; j < A.ColumnCount; j++)
{
Assert.AreEqual(ACopy[i, j], A[i, j]);
}
}
// Make sure b didn't change.
for (int i = 0; i < order; i++)
{
Assert.AreEqual(bCopy[i], b[i]);
}
}
[Test]
[Row(1, 1)]
[Row(2, 4)]
[Row(5, 8)]
[Row(10, 3)]
[Row(50, 10)]
[Row(100, 100)]
[MultipleAsserts]
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int col)
{
var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row);
var ACopy = A.Clone();
var chol = A.Cholesky();
var B = MatrixLoader.GenerateRandomMatrix(row, col);
var BCopy = B.Clone();
var X = new DenseMatrix(row, col);
chol.Solve(B, X);
Assert.AreEqual(B.RowCount, X.RowCount);
Assert.AreEqual(B.ColumnCount, X.ColumnCount);
var BReconstruct = A * X;
// Check the reconstruction.
for (int i = 0; i < B.RowCount; i++)
{
for (int j = 0; j < B.ColumnCount; j++)
{
Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
}
}
// Make sure A didn't change.
for (int i = 0; i < A.RowCount; i++)
{
for (int j = 0; j < A.ColumnCount; j++)
{
Assert.AreEqual(ACopy[i, j], A[i, j]);
}
}
// Make sure B didn't change.
for (int i = 0; i < B.RowCount; i++)
{
for (int j = 0; j < B.ColumnCount; j++)
{
Assert.AreEqual(BCopy[i, j], B[i, j]);
}
}
}*/
}
}

1
src/UnitTests/UnitTests.csproj

@ -117,6 +117,7 @@
<Compile Include="InterpolationTests\InterpolationInfrastructureContract.cs" />
<Compile Include="InterpolationTests\InterpolationFunctionalTest.cs" />
<Compile Include="LinearAlgebraTests\Double\Factorization\CholeskyTests.cs" />
<Compile Include="LinearAlgebraTests\Double\Factorization\LUTests.cs" />
<Compile Include="LinearAlgebraTests\Double\MatrixLoader.cs" />
<Compile Include="LinearAlgebraTests\Double\MatrixTests.Arithmetic.cs" />
<Compile Include="LinearAlgebraTests\Double\LinearAlgebraProviderTests.cs" />

Loading…
Cancel
Save