diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs
index def67602..b7d363f0 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs
+++ b/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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- public void LUFactor(Complex32[] a, int[] ipiv)
+ public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
index d42ac4cf..02a536f4 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
@@ -208,17 +208,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// The c matrix.
void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, T alpha, T[] a,
int aRows, int aColumns, T[] b, int bRows, int bColumns, T beta, T[] c);
-
+
///
/// Computes the LUP factorization of A. P*A = L*U.
///
- /// An by matrix. The matrix is overwritten with the
- /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
- /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
- /// The order of the square matrix .
- /// On exit, it contains the pivot indices. The size of the array must be .
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- void LUFactor(T[] a, int aOrder, int[] ipiv);
+ void LUFactor(T[] data, int order, int[] ipiv);
///
/// Computes the inverse of matrix using LU factorization.
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
index ae47b7df..d859b219 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
@@ -691,15 +691,79 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
///
/// Computes the LUP factorization of A. P*A = L*U.
///
- /// An by matrix. The matrix is overwritten with the
- /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
- /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
- /// The order of the square matrix .
- /// On exit, it contains the pivot indices. The size of the array must be .
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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)
+ ///
+ /// Computes the LUP factorization of A. P*A = L*U.
+ ///
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
+ /// This is equivalent to the GETRF LAPACK routine.
+ 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)
+ ///
+ /// Computes the LUP factorization of A. P*A = L*U.
+ ///
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
+ /// This is equivalent to the GETRF LAPACK routine.
+ 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)
+ ///
+ /// Computes the LUP factorization of A. P*A = L*U.
+ ///
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
+ /// This is equivalent to the GETRF LAPACK routine.
+ public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs
index 7154353e..998fd41b 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs
+++ b/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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- public void LUFactor(Complex32[] a, int[] ipiv)
+ public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
index 9ee36977..e1fbfddb 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
+++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
@@ -337,13 +337,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
///
/// Computes the LUP factorization of A. P*A = L*U.
///
- /// An by matrix. The matrix is overwritten with the
- /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
- /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
- /// The order of the square matrix .
- /// On exit, it contains the pivot indices. The size of the array must be .
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- 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);
}
-
+
///
- /// Computes the LU factorization of A.
+ /// Computes the LUP factorization of A. P*A = L*U.
///
- /// An m by n matrix. The matrix is overwritten with the
- /// the LU factorization On exit.
- /// On exit, it contains the pivot indices. The size
- /// of the array must be min(m,n).
+ /// An by matrix. The matrix is overwritten with the
+ /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
+ /// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
+ /// The order of the square matrix .
+ /// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
- public void LUFactor(Complex32[] a, int[] ipiv)
+ public void LUFactor(Complex32[] data, int order, int[] ipiv)
{
throw new NotImplementedException();
}
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs
index fe56050e..d6fecc66 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs
@@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
public static class ExtensionMethods
{
///
- /// Computes the cholesky decomposition for a matrix.
+ /// Computes the Cholesky decomposition for a matrix.
///
/// The matrix to factor.
/// The Cholesky decomposition object.
@@ -47,5 +47,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{
return Factorization.Cholesky.Create(matrix);
}
+
+ ///
+ /// Computes the LU decomposition for a matrix.
+ ///
+ /// The matrix to factor.
+ /// The LU decomposition object.
+ public static LU LU(this Matrix matrix)
+ {
+ return Factorization.LU.Create(matrix);
+ }
}
}
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
index de7a25c2..f9a6b3ba 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
@@ -167,5 +167,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
/// The right hand side vector, b.
/// The left hand side , x.
public abstract void Solve(Vector input, Vector result);
+
+ ///
+ /// Pivot a matrix according to this LU decomposition.
+ ///
+ /// The matrix to pivot.
+ 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);
+ }
+ }
+ }
+ }
}
}
diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs
index d78c7cb4..99cc3bc1 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs
+++ b/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++)
diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs
new file mode 100644
index 00000000..e97c7c04
--- /dev/null
+++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs
@@ -0,0 +1,332 @@
+//
+// 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.
+//
+
+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]);
+ }
+ }
+ }*/
+ }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 5a65d53a..ab2aef88 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -117,6 +117,7 @@
+