From de4eb2a7a53a2ac00cbb295efe113523dad79078 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Sun, 27 Feb 2011 22:35:34 +0800 Subject: [PATCH] native: added Cholesky solve and removed redundant number of B parameter from its interface --- src/NativeWrappers/MKL/lapack.cpp | 16 +-- .../ILinearAlgebraProviderOfT.cs | 24 +++-- .../ManagedLinearAlgebraProvider.Complex.cs | 36 ++++--- .../ManagedLinearAlgebraProvider.Complex32.cs | 38 +++++--- .../ManagedLinearAlgebraProvider.Double.cs | 38 +++++--- .../ManagedLinearAlgebraProvider.Single.cs | 38 +++++--- .../LinearAlgebra/native.generic.include | 70 ++++++++++--- .../LinearAlgebra/safe.native.common.include | 25 +++++ src/Numerics/Constants.cs | 9 ++ .../Complex/Factorization/DenseCholesky.cs | 4 +- .../Complex32/Factorization/DenseCholesky.cs | 4 +- .../Double/Factorization/DenseCholesky.cs | 4 +- .../Single/Factorization/DenseCholesky.cs | 4 +- .../Double/LinearAlgebraProviderTests.cs | 97 ++++++++++++++----- 14 files changed, 286 insertions(+), 121 deletions(-) diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 080332c4..33ab37e3 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -417,14 +417,14 @@ extern "C" { std::memcpy(clone, a, n*n*sizeof(float)); char uplo = 'L'; int info = 0; - SPOTRF(&uplo, &n, a, &n, &info); + SPOTRF(&uplo, &n, clone, &n, &info); if (info != 0){ delete[] clone; return info; } - SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + SPOTRS(&uplo, &n, &nrhs, clone, &n, b, &n, &info); return info; } @@ -434,14 +434,14 @@ extern "C" { std::memcpy(clone, a, n*n*sizeof(double)); char uplo = 'L'; int info = 0; - DPOTRF(&uplo, &n, a, &n, &info); + DPOTRF(&uplo, &n, clone, &n, &info); if (info != 0){ delete[] clone; return info; } - DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + DPOTRS(&uplo, &n, &nrhs, clone, &n, b, &n, &info); return info; } @@ -451,14 +451,14 @@ extern "C" { std::memcpy(clone, a, n*n*sizeof(MKL_Complex8)); char uplo = 'L'; int info = 0; - CPOTRF(&uplo, &n, a, &n, &info); + CPOTRF(&uplo, &n, clone, &n, &info); if (info != 0){ delete[] clone; return info; } - CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + CPOTRS(&uplo, &n, &nrhs, clone, &n, b, &n, &info); return info; } @@ -468,14 +468,14 @@ extern "C" { std::memcpy(clone, a, n*n*sizeof(MKL_Complex16)); char uplo = 'L'; int info = 0; - ZPOTRF(&uplo, &n, a, &n, &info); + ZPOTRF(&uplo, &n, clone, &n, &info); if (info != 0){ delete[] clone; return info; } - ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + ZPOTRS(&uplo, &n, &nrhs, clone, &n, b, &n, &info); return info; } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 08bea682..e0a9b4c6 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -268,7 +268,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. void LUSolve(int columnsOfB, T[] a, int order, T[] b); @@ -279,7 +279,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. void LUSolveFactored(int columnsOfB, T[] a, int order, int[] ipiv, T[] b); @@ -297,22 +297,20 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - void CholeskySolve(T[] a, int orderA, T[] b, int rowsB, int columnsB); + void CholeskySolve(T[] a, int orderA, T[] b, int columnsB); /// /// Solves A*X=B for X using a previously factored A matrix. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - void CholeskySolveFactored(T[] a, int orderA, T[] b, int rowsB, int columnsB); + void CholeskySolveFactored(T[] a, int orderA, T[] b, int columnsB); /// /// Computes the QR factorization of A. @@ -364,7 +362,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. /// The work array. The array must have a length of at least N, @@ -379,7 +377,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] b, int columnsB, T[] x); @@ -427,7 +425,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The singular values of A in ascending value. /// On exit U contains the left singular vectors. /// On exit VT contains the transposed right singular vectors. - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. void SvdSolve(T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x); @@ -441,7 +439,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The singular values of A in ascending value. /// On exit U contains the left singular vectors. /// On exit VT contains the transposed right singular vectors. - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. /// The work array. For real matrices, the work array should be at least @@ -458,7 +456,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The s values returned by . /// The left singular vectors returned by . /// The right singular vectors returned by . - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. void SvdSolveFactored(int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x); diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 43027bba..06e497a9 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { @@ -1050,6 +1050,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + var ipiv = new int[order]; var clone = new Complex[a.Length]; Array.Copy(a, 0, clone, 0, a.Length); @@ -1064,7 +1069,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { @@ -1098,6 +1103,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + // Compute the column vector P*B for (var i = 0; i < ipiv.Length; i++) { @@ -1244,11 +1254,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(Complex[] a, int orderA, Complex[] b, int rowsB, int columnsB) + public virtual void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -1260,9 +1269,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) @@ -1270,8 +1279,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CholeskyFactor(a, orderA); - CholeskySolveFactored(a, orderA, b, rowsB, columnsB); + var clone = new Complex[a.Length]; + Array.Copy(a, 0, clone, 0, a.Length); + CholeskyFactor(clone, orderA); + CholeskySolveFactored(clone, orderA, b, columnsB); } /// @@ -1280,10 +1291,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int rowsB, int columnsB) + public virtual void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -1295,9 +1305,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index e284ba47..8a8f6b99 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { @@ -1050,6 +1050,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + var ipiv = new int[order]; var clone = new Complex32[a.Length]; Array.Copy(a, 0, clone, 0, a.Length); @@ -1064,7 +1069,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { @@ -1098,6 +1103,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + // Compute the column vector P*B for (var i = 0; i < ipiv.Length; i++) { @@ -1244,11 +1254,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int rowsB, int columnsB) + public virtual void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -1260,9 +1269,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) @@ -1270,8 +1279,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CholeskyFactor(a, orderA); - CholeskySolveFactored(a, orderA, b, rowsB, columnsB); + var clone = new Complex32[a.Length]; + Array.Copy(a, 0, clone, 0, a.Length); + CholeskyFactor(clone, orderA); + CholeskySolveFactored(clone, orderA, b, columnsB); } /// @@ -1279,11 +1290,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int rowsB, int columnsB) + public virtual void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -1295,9 +1305,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index ea906a31..10ec9bd3 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, double[] a, int order, double[] b) { @@ -1044,6 +1044,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + var ipiv = new int[order]; var clone = new double[a.Length]; Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); @@ -1058,7 +1063,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { @@ -1092,6 +1097,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + // Compute the column vector P*B for (var i = 0; i < ipiv.Length; i++) { @@ -1238,11 +1248,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(double[] a, int orderA, double[] b, int rowsB, int columnsB) + public virtual void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -1254,9 +1263,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) @@ -1264,8 +1273,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CholeskyFactor(a, orderA); - CholeskySolveFactored(a, orderA, b, rowsB, columnsB); + var clone = new double[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + CholeskyFactor(clone, orderA); + CholeskySolveFactored(clone, orderA, b, columnsB); } /// @@ -1273,11 +1284,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. Has to be different than . /// The number of rows and columns in A. - /// The B matrix. Has to be different than . - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(double[] a, int orderA, double[] b, int rowsB, int columnsB) + public virtual void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -1289,9 +1299,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 8a125cc7..d9f99f7d 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -1021,7 +1021,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, float[] a, int order, float[] b) { @@ -1045,6 +1045,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + var ipiv = new int[order]; var clone = new float[a.Length]; Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); @@ -1059,7 +1064,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// On input the B matrix; on output the X matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { @@ -1092,6 +1097,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } // Compute the column vector P*B for (var i = 0; i < ipiv.Length; i++) @@ -1239,11 +1249,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(float[] a, int orderA, float[] b, int rowsB, int columnsB) + public virtual void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -1255,9 +1264,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) @@ -1265,8 +1274,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CholeskyFactor(a, orderA); - CholeskySolveFactored(a, orderA, b, rowsB, columnsB); + var clone = new float[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + CholeskyFactor(clone, orderA); + CholeskySolveFactored(clone, orderA, b, columnsB); } /// @@ -1274,11 +1285,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(float[] a, int orderA, float[] b, int rowsB, int columnsB) + public virtual void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -1290,9 +1300,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("b"); } - if (orderA != rowsB) + if (b.Length != orderA * columnsB) { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } if (ReferenceEquals(a, b)) diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index b3f61c7c..7da5fb3f 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include @@ -338,7 +338,7 @@ /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] public override void LUSolve(int columnsOfB, <#=dataType#>[] a, int order, <#=dataType#>[] b) @@ -357,6 +357,11 @@ { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } SafeNativeMethods.<#=prefix#>_lu_solve(order, columnsOfB, a, b); } @@ -368,7 +373,7 @@ /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] public override void LUSolveFactored(int columnsOfB, <#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] b) @@ -398,6 +403,11 @@ throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + SafeNativeMethods.<#=prefix#>_lu_solve_factored(order, columnsOfB, a, ipiv, b); } @@ -434,15 +444,34 @@ /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB) + public override void CholeskySolve(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int columnsB) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + SafeNativeMethods.<#=prefix#>_cholesky_solve(orderA, columnsB, a, b); } /// @@ -450,14 +479,33 @@ /// /// The square, positive definite matrix A. /// The number of rows and columns in A. - /// The B matrix. - /// The number of rows in the B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB) + public override void CholeskySolveFactored(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int columnsB) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + SafeNativeMethods.<#=prefix#>_cholesky_solve_factored(orderA, columnsB, a, b); } /// @@ -504,7 +552,7 @@ /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. - /// The B matrix. + /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. [SecuritySafeCritical] diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index 8d71cda7..adae155b 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include @@ -185,4 +185,29 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, float[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, double[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, Complex32[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, Complex[] b); + #endregion LAPACK diff --git a/src/Numerics/Constants.cs b/src/Numerics/Constants.cs index 513d20ea..e13b427c 100644 --- a/src/Numerics/Constants.cs +++ b/src/Numerics/Constants.cs @@ -177,6 +177,15 @@ namespace MathNet.Numerics /// public const int SizeOfFloat = sizeof(float); + /// + /// The size of a Complex in bytes. + /// + public const int SizeOfComplex = 2 * SizeOfDouble; + + /// + /// The size of a Complex in bytes. + /// + public const int SizeOfComplex32 = 2 * SizeOfFloat; #endregion #region UNIVERSAL CONSTANTS diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs index 8cf98d70..13d71b61 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs @@ -124,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.ColumnCount); } /// @@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, 1); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs index 204b1dc8..14ce07ff 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs @@ -124,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.ColumnCount); } /// @@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, 1); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index ef80191b..0562a2b2 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -122,7 +122,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.ColumnCount); } /// @@ -171,7 +171,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, 1); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs index e885d196..4ef94134 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs @@ -122,7 +122,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.ColumnCount); } /// @@ -171,7 +171,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; - Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); + Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, 1); } } } diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index bc70778b..7aa38802 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -393,32 +393,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } - /// - /// Can compute the Cholesky factorization. - /// - [Test] - public void CanComputeCholeskyFactor() - { - var matrix = new double[] { 1, 1, 1, 1, 1, 5, 5, 5, 1, 5, 14, 14, 1, 5, 14, 15 }; - Provider.CholeskyFactor(matrix, 4); - Assert.AreEqual(matrix[0], 1); - Assert.AreEqual(matrix[1], 1); - Assert.AreEqual(matrix[2], 1); - Assert.AreEqual(matrix[3], 1); - Assert.AreEqual(matrix[4], 0); - Assert.AreEqual(matrix[5], 2); - Assert.AreEqual(matrix[6], 2); - Assert.AreEqual(matrix[7], 2); - Assert.AreEqual(matrix[8], 0); - Assert.AreEqual(matrix[9], 0); - Assert.AreEqual(matrix[10], 3); - Assert.AreEqual(matrix[11], 3); - Assert.AreEqual(matrix[12], 0); - Assert.AreEqual(matrix[13], 0); - Assert.AreEqual(matrix[14], 0); - Assert.AreEqual(matrix[15], 1); - } - /// /// Can compute the LU factor of a matrix. /// @@ -598,6 +572,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14); } + /// + /// Can compute the Cholesky factorization. + /// + [Test] + public void CanComputeCholeskyFactor() + { + var matrix = new double[] { 1, 1, 1, 1, 1, 5, 5, 5, 1, 5, 14, 14, 1, 5, 14, 15 }; + Provider.CholeskyFactor(matrix, 4); + Assert.AreEqual(matrix[0], 1); + Assert.AreEqual(matrix[1], 1); + Assert.AreEqual(matrix[2], 1); + Assert.AreEqual(matrix[3], 1); + Assert.AreEqual(matrix[4], 0); + Assert.AreEqual(matrix[5], 2); + Assert.AreEqual(matrix[6], 2); + Assert.AreEqual(matrix[7], 2); + Assert.AreEqual(matrix[8], 0); + Assert.AreEqual(matrix[9], 0); + Assert.AreEqual(matrix[10], 3); + Assert.AreEqual(matrix[11], 3); + Assert.AreEqual(matrix[12], 0); + Assert.AreEqual(matrix[13], 0); + Assert.AreEqual(matrix[14], 0); + Assert.AreEqual(matrix[15], 1); + } + + /// + /// Can solve Ax=b using Cholesky factorization. + /// + [Test] + public void CanSolveUsingCholesky() + { + var matrix = new DenseMatrix(3, 3, new double[] { 1, 1, 1, 1, 2, 3, 1, 3, 6 }); + var a = new double[] { 1, 1, 1, 1, 2, 3, 1, 3, 6 }; + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + Provider.CholeskySolve(a, 3, b, 2); + + AssertHelpers.AlmostEqual(b[0], 0, 14); + AssertHelpers.AlmostEqual(b[1], 1, 14); + AssertHelpers.AlmostEqual(b[2], 0, 14); + AssertHelpers.AlmostEqual(b[3], 3, 14); + AssertHelpers.AlmostEqual(b[4], 1, 14); + AssertHelpers.AlmostEqual(b[5], 0, 14); + + NotModified(3, 3, a, matrix); + } + + /// + /// Can solve Ax=b using LU factorization using a factored matrix. + /// + [Test] + public void CanSolveUsingCholeskyOnFactoredMatrix() + { + var a = new double[] { 1, 1, 1, 1, 2, 3, 1, 3, 6 }; + + Provider.CholeskyFactor(a, 3); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + Provider.CholeskySolveFactored(a, 3, b, 2); + + AssertHelpers.AlmostEqual(b[0], 0, 14); + AssertHelpers.AlmostEqual(b[1], 1, 14); + AssertHelpers.AlmostEqual(b[2], 0, 14); + AssertHelpers.AlmostEqual(b[3], 3, 14); + AssertHelpers.AlmostEqual(b[4], 1, 14); + AssertHelpers.AlmostEqual(b[5], 0, 14); + } + + + /// /// Checks to see if a matrix and array contain the same values. ///