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.
///