From 9755692cf8e072b704e1fcb980ebac5e1a381e3d Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Tue, 18 May 2010 22:56:03 +0100 Subject: [PATCH] Small additions to LU and initial unit tests. --- .../Atlas/AtlasLinearAlgebraProvider.cs | 62 ++-- .../ILinearAlgebraProviderOfT.cs | 14 +- .../ManagedLinearAlgebraProvider.cs | 111 +++++- .../Mkl/MklLinearAlgebraProvider.cs | 62 ++-- .../NativeAlgebraProvider.include | 57 +-- .../Double/Factorization/ExtensionMethods.cs | 12 +- .../LinearAlgebra/Double/Factorization/LU.cs | 21 ++ .../Double/Factorization/CholeskyTests.cs | 3 + .../Double/Factorization/LUTests.cs | 332 ++++++++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 10 files changed, 572 insertions(+), 103 deletions(-) create mode 100644 src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs 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 @@ +