tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

native: added Cholesky solve and removed redundant number of B parameter from its interface

pull/36/head
Marcus Cuda 16 years ago
parent
d25e59247e
commit
de4eb2a7a5
14 changed files with 286 additions and 121 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 16
      src/NativeWrappers/MKL/lapack.cpp
  2. 24
      src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
  3. 36
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  4. 38
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  5. 38
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  6. 38
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  7. 70
      src/Numerics/Algorithms/LinearAlgebra/native.generic.include
  8. 25
      src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
  9. 9
      src/Numerics/Constants.cs
  10. 4
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs
  11. 4
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs
  12. 4
      src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs
  13. 4
      src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs
  14. 97
      src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

16
src/NativeWrappers/MKL/lapack.cpp
View File

@ -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;
}

24
src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
View File

@ -268,7 +268,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
void LUSolve(int columnsOfB, T[] a, int order, T[] b);
@ -279,7 +279,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
void LUSolveFactored(int columnsOfB, T[] a, int order, int[] ipiv, T[] b);
@ -297,22 +297,20 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
void CholeskySolve(T[] a, int orderA, T[] b, int rowsB, int columnsB);
void CholeskySolve(T[] a, int orderA, T[] b, int columnsB);
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
void CholeskySolveFactored(T[] a, int orderA, T[] b, int rowsB, int columnsB);
void CholeskySolveFactored(T[] a, int orderA, T[] b, int columnsB);
/// <summary>
/// Computes the QR factorization of A.
@ -364,7 +362,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
@ -379,7 +377,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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
/// <param name="s">The singular values of A in ascending value. </param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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
/// <param name="s">The singular values of A in ascending value. </param>
/// <param name="u">On exit U contains the left singular vectors.</param>
/// <param name="vt">On exit VT contains the transposed right singular vectors.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. For real matrices, the work array should be at least
@ -458,7 +456,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="s">The s values returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[])"/>.</param>
/// <param name="u">The left singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int, T[],T[],T[])"/>.</param>
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[],T[])"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
void SvdSolveFactored(int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] b, int columnsB, T[] x);

36
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
View File

@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
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
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
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
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
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);
}
/// <summary>
@ -1280,10 +1291,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
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))

38
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
View File

@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
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
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
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
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
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);
}
/// <summary>
@ -1279,11 +1290,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
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))

38
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
View File

@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
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
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
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
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
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);
}
/// <summary>
@ -1273,11 +1284,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// </summary>
/// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix. Has to be different than <paramref name="a"/>.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
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))

38
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
View File

@ -1021,7 +1021,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
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
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On input the B matrix; on output the X matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
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
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
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);
}
/// <summary>
@ -1274,11 +1285,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
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))

70
src/Numerics/Algorithms/LinearAlgebra/native.generic.include
View File

@ -338,7 +338,7 @@
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[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 @@
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[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 @@
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[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);
}
/// <summary>
@ -450,14 +479,33 @@
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">The B matrix.</param>
/// <param name="rowsB">The number of rows in the B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[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);
}
/// <summary>
@ -504,7 +552,7 @@
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="b">The B matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]

25
src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
View File

@ -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

9
src/Numerics/Constants.cs
View File

@ -177,6 +177,15 @@ namespace MathNet.Numerics
/// </summary>
public const int SizeOfFloat = sizeof(float);
/// <summary>
/// The size of a Complex in bytes.
/// </summary>
public const int SizeOfComplex = 2 * SizeOfDouble;
/// <summary>
/// The size of a Complex in bytes.
/// </summary>
public const int SizeOfComplex32 = 2 * SizeOfFloat;
#endregion
#region UNIVERSAL CONSTANTS

4
src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs
View File

@ -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);
}
/// <summary>
@ -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);
}
}
}

4
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs
View File

@ -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);
}
/// <summary>
@ -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);
}
}
}

4
src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs
View File

@ -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);
}
/// <summary>
@ -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);
}
}
}

4
src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs
View File

@ -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);
}
/// <summary>
@ -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);
}
}
}

97
src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
View File

@ -393,32 +393,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
}
}
/// <summary>
/// Can compute the <c>Cholesky</c> factorization.
/// </summary>
[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);
}
/// <summary>
/// Can compute the LU factor of a matrix.
/// </summary>
@ -598,6 +572,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14);
}
/// <summary>
/// Can compute the <c>Cholesky</c> factorization.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using Cholesky factorization.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using LU factorization using a factored matrix.
/// </summary>
[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);
}
/// <summary>
/// Checks to see if a matrix and array contain the same values.
/// </summary>

Write Preview
Loading…
Cancel
Save