diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs index 41c52f77..a2331ba1 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs @@ -347,8 +347,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(double[] a) + public void LUInverse(double[] a, int order) { throw new NotImplementedException(); } @@ -357,9 +358,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(double[] a, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -368,11 +370,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(double[] a, double[] work) + public void LUInverse(double[] a, int order, double[] work) { throw new NotImplementedException(); } @@ -381,12 +384,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(double[] a, int[] ipiv, double[] work) + public void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) { throw new NotImplementedException(); } @@ -396,9 +400,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, double[] a, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { throw new NotImplementedException(); } @@ -408,10 +413,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { throw new NotImplementedException(); } @@ -422,9 +428,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, int order, double[] b) { throw new NotImplementedException(); } @@ -435,10 +442,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { throw new NotImplementedException(); } @@ -958,7 +966,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } - /// /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -977,8 +984,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(float[] a) + public void LUInverse(float[] a, int order) { throw new NotImplementedException(); } @@ -987,9 +995,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(float[] a, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -998,11 +1007,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(float[] a, float[] work) + public void LUInverse(float[] a, int order, float[] work) { throw new NotImplementedException(); } @@ -1011,12 +1021,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(float[] a, int[] ipiv, float[] work) + public void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) { throw new NotImplementedException(); } @@ -1026,9 +1037,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, float[] a, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } @@ -1038,10 +1050,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } @@ -1052,9 +1065,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } @@ -1065,10 +1079,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } @@ -1605,8 +1620,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex[] a) + public void LUInverse(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1615,9 +1631,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex[] a, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -1626,11 +1643,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex[] a, Complex[] work) + public void LUInverse(Complex[] a, int order, Complex[] work) { throw new NotImplementedException(); } @@ -1639,12 +1657,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) { throw new NotImplementedException(); } @@ -1654,9 +1673,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } @@ -1666,10 +1686,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } @@ -1680,9 +1701,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } @@ -1693,10 +1715,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } @@ -2233,8 +2256,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex32[] a) + public void LUInverse(Complex32[] a, int order) { throw new NotImplementedException(); } @@ -2243,9 +2267,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex32[] a, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -2254,11 +2279,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex32[] a, Complex32[] work) + public void LUInverse(Complex32[] a, int order, Complex32[] work) { throw new NotImplementedException(); } @@ -2267,12 +2293,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex32[] a, int[] ipiv, Complex32[] work) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) { throw new NotImplementedException(); } @@ -2282,9 +2309,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } @@ -2294,10 +2322,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } @@ -2308,9 +2337,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } @@ -2321,10 +2351,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 3beb89c7..bf1a7532 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -226,56 +226,62 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - void LUInverse(T[] a); + void LUInverse(T[] a, int order); /// /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - void LUInverseFactored(T[] a, int[] ipiv); + void LUInverseFactored(T[] a, int order, int[] ipiv); /// /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - void LUInverse(T[] a, T[] work); + void LUInverse(T[] a, int order, T[] work); /// /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - void LUInverseFactored(T[] a, int[] ipiv, T[] work); + void LUInverseFactored(T[] a, int order, int[] ipiv, T[] work); /// /// Solves A*X=B for X using LU factorization. /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - void LUSolve(int columnsOfB, T[] a, T[] b); + void LUSolve(int columnsOfB, T[] a, int order, T[] b); /// /// Solves A*X=B for X using a previously factored A matrix. /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - void LUSolveFactored(int columnsOfB, T[] a, int ipiv, T[] b); + void LUSolveFactored(int columnsOfB, T[] a, int order, int[] ipiv, T[] b); /// /// Solves A*X=B for X using LU factorization. @@ -283,9 +289,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - void LUSolve(Transpose transposeA, int columnsOfB, T[] a, T[] b); + void LUSolve(Transpose transposeA, int columnsOfB, T[] a, int order, T[] b); /// /// Solves A*X=B for X using a previously factored A matrix. @@ -293,10 +300,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - void LUSolveFactored(Transpose transposeA, int columnsOfB, T[] a, int ipiv, T[] b); + void LUSolveFactored(Transpose transposeA, int columnsOfB, T[] a, int order, int[] ipiv, T[] b); /// /// Computes the Cholesky factorization of A. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 6fcf52a5..54c8a990 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -724,6 +724,26 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRF LAPACK routine. public void LUFactor(double[] data, int order, int[] ipiv) { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + // Initialize the pivot matrix to the identity permutation. for (var i = 0; i < order; i++) { @@ -794,44 +814,322 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - public void LUInverse(double[] a) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(double[] a, int order) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUInverseFactored(a, order, ipiv); } - public void LUInverseFactored(double[] a, int[] ipiv) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(double[] a, int order, int[] ipiv) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + var inverse = new double[a.Length]; + for (var i = 0; i < order; i++) + { + inverse[i + order * i] = 1.0; + } + + LUSolveFactored(order, a, order, ipiv, inverse); + Buffer.BlockCopy(inverse, 0, a, 0, a.Length * Constants.SizeOfDouble); } - public void LUInverse(double[] a, double[] work) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(double[] a, int order, double[] work) { - throw new NotImplementedException(); + LUInverse(a, order); } - public void LUInverseFactored(double[] a, int[] ipiv, double[] work) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) { - throw new NotImplementedException(); + LUInverseFactored(a, order, ipiv); } - public void LUSolve(int columnsOfB, double[] a, double[] b) + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUSolveFactored(columnsOfB, a, order, ipiv, b); } - public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b) + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// The B matrix. + /// This is equivalent to the GETRS LAPACK routine. + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + // Compute the column vector P*B + for (var i = 0; i < ipiv.Length; i++) + { + if (ipiv[i] == i) + { + continue; + } + + var p = ipiv[i]; + for (var j = 0; j < columnsOfB; j++) + { + var indexk = j*order; + var indexkp = indexk + p; + var indexkj = indexk + i; + var temp = b[indexkp]; + b[indexkp] = b[indexkj]; + b[indexkj] = temp; + } + } + + // Solve L*Y = P*B + for (var k = 0; k < order; k++) + { + var korder = k * order; + for (var i = k + 1; i < order; i++) + { + for (var j = 0; j < columnsOfB; j++) + { + var index = j * order; + b[i + index] -= b[k + index] * a[i + korder]; + } + } + } + + // Solve U*X = Y; + for (var k = order - 1; k >= 0; k--) + { + var korder = k + k * order; + for (var j = 0; j < columnsOfB; j++) + { + b[k + j * order] /= a[korder]; + } + + korder = k * order; + for (var i = 0; i < k; i++) + { + for (var j = 0; j < columnsOfB; j++) + { + var index = j * order; + b[i + index] -= b[k + index] * a[i + korder]; + } + } + } } - public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b) + /// + /// Solves A*X=B for X using LU factorization. + /// + /// How to transpose the matrix. + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, int order, double[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); } - public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b) + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// How to transpose the matrix. + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// The B matrix. + /// This is equivalent to the GETRS LAPACK routine. + public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if ((transposeA == Transpose.Transpose) || (transposeA == Transpose.ConjugateTranspose)) + { + var aT = new double[a.Length]; + for (var i = 0; i < order; i++) + { + for (var j = 0; j < order; j++) + { + aT[(j * order) + i] = a[(i * order) + j]; + } + } + LUSolveFactored(columnsOfB, aT, order, ipiv, b); + } + else + { + LUSolveFactored(columnsOfB, a, order, ipiv, b); + } } /// @@ -3724,42 +4022,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void LUInverse(float[] a) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(float[] a, int order) { throw new NotImplementedException(); } - public void LUInverseFactored(float[] a, int[] ipiv) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(float[] a, int order, int[] ipiv) { throw new NotImplementedException(); } - public void LUInverse(float[] a, float[] work) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(float[] a, int order, float[] work) { throw new NotImplementedException(); } - public void LUInverseFactored(float[] a, int[] ipiv, float[] work) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) { throw new NotImplementedException(); } - public void LUSolve(int columnsOfB, float[] a, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } - public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } @@ -4677,42 +5007,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void LUInverse(Complex[] a) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(Complex[] a, int order) { throw new NotImplementedException(); } - public void LUInverseFactored(Complex[] a, int[] ipiv) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { throw new NotImplementedException(); } - public void LUInverse(Complex[] a, Complex[] work) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(Complex[] a, int order, Complex[] work) { throw new NotImplementedException(); } - public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) { throw new NotImplementedException(); } - public void LUSolve(int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } - public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } @@ -5595,42 +5957,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void LUInverse(Complex32[] a) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(Complex32[] a, int order) { throw new NotImplementedException(); } - public void LUInverseFactored(Complex32[] a, int[] ipiv) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { throw new NotImplementedException(); } - public void LUInverse(Complex32[] a, Complex32[] work) + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + public void LUInverse(Complex32[] a, int order, Complex32[] work) { throw new NotImplementedException(); } - public void LUInverseFactored(Complex32[] a, int[] ipiv, Complex32[] work) + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRI LAPACK routine. + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) { throw new NotImplementedException(); } - public void LUSolve(int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } - public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 43d26741..dbe3826d 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -346,8 +346,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(double[] a) + public void LUInverse(double[] a, int order) { throw new NotImplementedException(); } @@ -356,9 +357,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(double[] a, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -367,11 +369,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(double[] a, double[] work) + public void LUInverse(double[] a, int order, double[] work) { throw new NotImplementedException(); } @@ -380,12 +383,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(double[] a, int[] ipiv, double[] work) + public void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) { throw new NotImplementedException(); } @@ -395,9 +399,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, double[] a, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { throw new NotImplementedException(); } @@ -407,10 +412,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { throw new NotImplementedException(); } @@ -421,9 +427,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, int order, double[] b) { throw new NotImplementedException(); } @@ -434,10 +441,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { throw new NotImplementedException(); } @@ -957,7 +965,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } - /// /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -976,8 +983,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(float[] a) + public void LUInverse(float[] a, int order) { throw new NotImplementedException(); } @@ -986,9 +994,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(float[] a, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -997,11 +1006,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(float[] a, float[] work) + public void LUInverse(float[] a, int order, float[] work) { throw new NotImplementedException(); } @@ -1010,12 +1020,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(float[] a, int[] ipiv, float[] work) + public void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) { throw new NotImplementedException(); } @@ -1025,9 +1036,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, float[] a, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } @@ -1037,10 +1049,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } @@ -1051,9 +1064,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, int order, float[] b) { throw new NotImplementedException(); } @@ -1064,10 +1078,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { throw new NotImplementedException(); } @@ -1604,8 +1619,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex[] a) + public void LUInverse(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1614,9 +1630,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex[] a, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -1625,11 +1642,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex[] a, Complex[] work) + public void LUInverse(Complex[] a, int order, Complex[] work) { throw new NotImplementedException(); } @@ -1638,12 +1656,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) { throw new NotImplementedException(); } @@ -1653,9 +1672,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } @@ -1665,10 +1685,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } @@ -1679,9 +1700,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) { throw new NotImplementedException(); } @@ -1692,10 +1714,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { throw new NotImplementedException(); } @@ -2232,8 +2255,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex32[] a) + public void LUInverse(Complex32[] a, int order) { throw new NotImplementedException(); } @@ -2242,9 +2266,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex32[] a, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { throw new NotImplementedException(); } @@ -2253,11 +2278,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of matrix using LU factorization. /// /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. - public void LUInverse(Complex32[] a, Complex32[] work) + public void LUInverse(Complex32[] a, int order, Complex32[] work) { throw new NotImplementedException(); } @@ -2266,12 +2292,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the inverse of a previously factored matrix. /// /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . /// The pivot indices of . /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. - public void LUInverseFactored(Complex32[] a, int[] ipiv, Complex32[] work) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) { throw new NotImplementedException(); } @@ -2281,9 +2308,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } @@ -2293,10 +2321,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } @@ -2307,9 +2336,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The square matrix A. + /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, Complex32[] b) + public void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, int order, Complex32[] b) { throw new NotImplementedException(); } @@ -2320,10 +2350,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// How to transpose the matrix. /// The number of columns of B. /// The factored A matrix. + /// The order of the square matrix . /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int ipiv, Complex32[] b) + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { throw new NotImplementedException(); } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs index 3283fb06..8524e4ef 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -122,8 +122,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; - throw new NotImplementedException(); - //Control.LinearAlgebraProvider.LUSolveFactored(dfactors.Data, dfactors.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); + Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); } /// @@ -171,9 +170,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // LU solve by overwriting result. - var dfactors = Factors as DenseMatrix; - throw new NotImplementedException(); - //Control.LinearAlgebraProvider.LUSolveFactored(dfactors.Data, dfactors.RowCount, dresult.Data, dresult.Count, 1); + var dfactors = (DenseMatrix)Factors; + Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); + } + + /// + /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition. + /// + /// The inverse of this matrix. + public override Matrix Inverse() + { + var result = (DenseMatrix)Factors.Clone(); + Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots); + return result; } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs index 465330fd..2b0c35e5 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -186,5 +186,11 @@ 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); + + /// + /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition. + /// + /// The inverse of this matrix. + public abstract Matrix Inverse(); } } diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs index 357fb0a8..1ee9c6b2 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs @@ -28,6 +28,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { using System; using Distributions; + using Factorization; using Properties; using Threading; @@ -831,7 +832,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// effective numerical rank, obtained from SVD public virtual int Rank() { - throw new NotImplementedException(); + var svd = this.Svd(false); + return svd.Rank; } /// Calculates the condition number of this matrix. @@ -839,14 +841,34 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The condition number is calculated using singular value decomposition. public virtual double ConditionNumber() { - throw new NotImplementedException(); + var svd = this.Svd(false); + return svd.ConditionNumber; } /// Computes the determinant of this matrix. /// The determinant of this matrix. public virtual double Determinant() { - throw new NotImplementedException(); + if (RowCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var lu = this.LU(); + return lu.Determinant; + } + + /// Computes the inverse of this matrix. + /// The inverse of this matrix. + public virtual Matrix Inverse() + { + if (RowCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var lu = this.LU(); + return lu.Inverse(); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index f5888801..e279548d 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -1088,12 +1088,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (columnLength > subMatrix.ColumnCount) { - throw new ArgumentOutOfRangeException("columnLength", "columnLength can be at most the number of columns in subMatrix."); + throw new ArgumentOutOfRangeException("columnLength", @"columnLength can be at most the number of columns in subMatrix."); } if (rowLength > subMatrix.RowCount) { - throw new ArgumentOutOfRangeException("rowLength", "rowLength can be at most the number of rows in subMatrix."); + throw new ArgumentOutOfRangeException("rowLength", @"rowLength can be at most the number of rows in subMatrix."); } var colMax = columnIndex + columnLength; diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs index 66b6a183..2fb150ca 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { - using System.Collections.Generic; using MbUnit.Framework; using LinearAlgebra.Double; using LinearAlgebra.Double.Factorization; @@ -43,44 +42,30 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [Row(100)] public void CanFactorizeIdentity(int order) { - var I = DenseMatrix.Identity(order); - var lu = I.LU(); + var matrixI = DenseMatrix.Identity(order); + var factorLU = matrixI.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++) + var matrixL = factorLU.L; + Assert.AreEqual(matrixI.RowCount, matrixL.RowCount); + Assert.AreEqual(matrixI.ColumnCount, matrixL.ColumnCount); + for (var i = 0; i < matrixL.RowCount; i++) { - for (var j = 0; j < L.ColumnCount; j++) + for (var j = 0; j < matrixL.ColumnCount; j++) { - if (i == j) - { - Assert.AreEqual(1.0, L[i, j]); - } - else - { - Assert.AreEqual(0.0, L[i, j]); - } + Assert.AreEqual(i == j ? 1.0 : 0.0, matrixL[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++) + var matrixU = factorLU.U; + Assert.AreEqual(matrixI.RowCount, matrixU.RowCount); + Assert.AreEqual(matrixI.ColumnCount, matrixU.ColumnCount); + for (var i = 0; i < matrixU.RowCount; i++) { - for (var j = 0; j < U.ColumnCount; j++) + for (var j = 0; j < matrixU.ColumnCount; j++) { - if (i == j) - { - Assert.AreEqual(1.0, U[i, j]); - } - else - { - Assert.AreEqual(0.0, U[i, j]); - } + Assert.AreEqual(i == j ? 1.0 : 0.0, matrixU[i, j]); } } } @@ -92,7 +77,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization public void LUFailsWithNonSquareMatrix(int row, int col) { var I = new DenseMatrix(row, col); - var lu = I.LU(); + I.LU(); } [Test] @@ -116,47 +101,264 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanFactorizeRandomMatrix(int order) { - var X = MatrixLoader.GenerateRandomMatrix(order, order); - var lu = X.LU(); - var L = lu.L; - var U = lu.U; + var matrixX = MatrixLoader.GenerateRandomMatrix(order, order); + var factorLU = matrixX.LU(); + var matrixL = factorLU.L; + var matrixU = factorLU.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); + Assert.AreEqual(order, matrixL.RowCount); + Assert.AreEqual(order, matrixL.ColumnCount); + Assert.AreEqual(order, matrixU.RowCount); + Assert.AreEqual(order, matrixU.ColumnCount); // Make sure the L factor is lower triangular. - for (int i = 0; i < L.RowCount; i++) + for (var i = 0; i < matrixL.RowCount; i++) { - Assert.AreEqual(1.0, L[i, i]); - for (int j = i+1; j < L.ColumnCount; j++) + Assert.AreEqual(1.0, matrixL[i, i]); + for (var j = i+1; j < matrixL.ColumnCount; j++) { - Assert.AreEqual(0.0, L[i, j]); + Assert.AreEqual(0.0, matrixL[i, j]); } } // Make sure the U factor is upper triangular. - for (int i = 0; i < L.RowCount; i++) + for (var i = 0; i < matrixL.RowCount; i++) + { + for (var j = 0; j < i; j++) + { + Assert.AreEqual(0.0, matrixU[i, j]); + } + } + + // Make sure the LU factor times it's transpose is the original matrix. + var matrixXfromLU = matrixL * matrixU; + var permutationInverse = factorLU.P.Inverse(); + matrixXfromLU.PermuteRows(permutationInverse); + for (var i = 0; i < matrixXfromLU.RowCount; i++) + { + for (var j = 0; j < matrixXfromLU.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixX[i, j], matrixXfromLU[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 matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var vectorb = MatrixLoader.GenerateRandomVector(order); + var resultx = factorLU.Solve(vectorb); + + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < i; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreEqual(0.0, U[i, j]); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixX = factorLU.Solve(matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); } } - // Make sure the cholesky factor times it's transpose is the original matrix. - var XfromLU = L * U; - var Pinv = lu.P.Inverse(); - XfromLU.PermuteRows(Pinv); - for (int i = 0; i < XfromLU.RowCount; i++) + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < XfromLU.ColumnCount; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreApproximatelyEqual(X[i, j], XfromLU[i, j], 1.0e-11); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); } } } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomVectorWhenResultVectorGiven(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorbCopy = vectorb.Clone(); + var resultx = new DenseVector(order); + factorLU.Solve(vectorb, resultx); + + Assert.AreEqual(vectorb.Count, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure b didn't change. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreEqual(vectorbCopy[i], vectorb[i]); + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new DenseMatrix(order, order); + factorLU.Solve(matrixB, matrixX); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); + } + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure B didn't change. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); + } + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanInverse(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixAInverse = factorLU.Inverse(); + + // The inverse dimension is equal A + Assert.AreEqual(matrixAInverse.RowCount, matrixAInverse.RowCount); + Assert.AreEqual(matrixAInverse.ColumnCount, matrixAInverse.ColumnCount); + + var matrixIdentity = matrixA * matrixAInverse; + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Check if multiplication of A and AI produced identity matrix. + for (var i = 0; i < matrixIdentity.RowCount; i++) + { + Assert.AreApproximatelyEqual(matrixIdentity[i, i], 1.0, 1.0e-11); + } + } } }