From df73c6e828caa9922c149fd8060f419b4b328d5e Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Thu, 9 Dec 2010 21:48:52 +0800 Subject: [PATCH] native: added lu bug:BulirschStoerRationalInterpolation was always returing NaN --- build/mathnet-numerics.shfbproj | 21 +- .../Double/LinearAlgebraProviderTests.cs | 39 +- src/NativeWrappers/MKL/lapack.cpp | 557 +++++++++++++++--- .../MKLWrapper32Tests.csproj | 5 +- .../LinearAlgebra/native.generic.include | 27 +- .../LinearAlgebra/safe.native.common.include | 22 +- .../BulirschStoerRationalInterpolation.cs | 2 +- src/Numerics/Permutation.cs | 4 +- 8 files changed, 569 insertions(+), 108 deletions(-) diff --git a/build/mathnet-numerics.shfbproj b/build/mathnet-numerics.shfbproj index 8fd332bd..0cf1654c 100644 --- a/build/mathnet-numerics.shfbproj +++ b/build/mathnet-numerics.shfbproj @@ -21,19 +21,26 @@ vs2005 - + + True Numerical Library for .NET Copyright %28c%29 2009-2010 Math.NET - True - Math.NET Numerics + False + Math.NET Numerics Beta 1 API Reference False - + - + + AboveNamespaces + + + + False + HtmlHelp1, MSHelpViewer diff --git a/src/MSUnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/MSUnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index 1652aa29..c6a089dc 100644 --- a/src/MSUnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/MSUnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -73,11 +73,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double /// private readonly IDictionary _matrices = new Dictionary { - { "Singular3x3", new DenseMatrix(new[,] { { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 } }) }, - { "Square3x3", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 }, { -4.4, 5.5, 6.6 } }) }, - { "Square4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { 0.0, 1.1, 2.2, 3.3 }, { 1.0, 2.1, 6.2, 4.3 }, { -4.4, 5.5, 6.6, -7.7 } }) }, - { "Singular4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 } }) }, - { "Tall3x2", new DenseMatrix(new[,] { { -1.1, -2.2 }, { 0.0, 1.1 }, { -4.4, 5.5 } }) }, + { "Singular3x3", new DenseMatrix(new[,] { { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 } }) }, + { "Square3x3", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 }, { -4.4, 5.5, 6.6 } }) }, + { "Square4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { 0.0, 1.1, 2.2, 3.3 }, { 1.0, 2.1, 6.2, 4.3 }, { -4.4, 5.5, 6.6, -7.7 } }) }, + { "Singular4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 } }) }, + { "Tall3x2", new DenseMatrix(new[,] { { -1.1, -2.2 }, { 0.0, 1.1 }, { -4.4, 5.5 } }) }, { "Wide2x3", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 } }) } }; @@ -406,5 +406,34 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double Assert.AreEqual(matrix[14], 0); Assert.AreEqual(matrix[15], 1); } + + /// + /// Can compute the LU factor of a matrix. + /// + [TestMethod] + public void CanComputeLuFactor() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.RowCount]; + Array.Copy(matrix.Data, a, a.Length); + + var ipiv = new int[matrix.RowCount]; + + Provider.LUFactor(a, matrix.RowCount, ipiv); + + AssertHelpers.AlmostEqual(a[0], -4.4, 15); + AssertHelpers.AlmostEqual(a[1], 0.25, 15); + AssertHelpers.AlmostEqual(a[2], 0, 15); + AssertHelpers.AlmostEqual(a[3], 5.5, 15); + AssertHelpers.AlmostEqual(a[4], -3.575, 15); + AssertHelpers.AlmostEqual(a[5], -0.307692307692308, 15); + AssertHelpers.AlmostEqual(a[6], 6.6, 15); + AssertHelpers.AlmostEqual(a[7], -4.95, 15); + AssertHelpers.AlmostEqual(a[8], 0.676923076923077, 15); + + Assert.AreEqual(ipiv[0], 2); + Assert.AreEqual(ipiv[1], 2); + Assert.AreEqual(ipiv[2], 2); + } } } \ No newline at end of file diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 350860c9..2159c0aa 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -4,89 +4,476 @@ extern "C" { - DLLEXPORT int s_cholesky_factor(int n, float a[]){ - char uplo = 'L'; - int info = 0; - SPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; - } - - DLLEXPORT int d_cholesky_factor(int n, double* a){ - char uplo = 'L'; - int info = 0; - DPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; - } - - DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ - char uplo = 'L'; - int info = 0; - Complex8 zero; - zero.real = 0.0; - zero.real = 0.0; - CPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; - } - - DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ - char uplo = 'L'; - int info = 0; - Complex16 zero; - zero.real = 0.0; - zero.real = 0.0; - ZPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; - } - - DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[]) - { - return SLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[]) - { - return DLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[]) - { - return CLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[]) - { - return ZLANGE(&norm, &m, &n, a, &m, work); - } + DLLEXPORT int s_cholesky_factor(int n, float a[]){ + char uplo = 'L'; + int info = 0; + SPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = 0; + } + } + return info; + } + + DLLEXPORT int d_cholesky_factor(int n, double* a){ + char uplo = 'L'; + int info = 0; + DPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = 0; + } + } + return info; + } + + DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ + char uplo = 'L'; + int info = 0; + Complex8 zero; + zero.real = 0.0; + zero.real = 0.0; + CPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + return info; + } + + DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ + char uplo = 'L'; + int info = 0; + Complex16 zero; + zero.real = 0.0; + zero.real = 0.0; + ZPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + return info; + } + + DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[]) + { + return SLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[]) + { + return DLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[]) + { + return CLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[]) + { + return ZLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT void s_lu_factor(int m, float a[], int ipiv[]) + { + int info; + SGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_factor(int m, double a[], int ipiv[]) + { + int info; + DGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_factor(int m, MKL_Complex8 a[], int ipiv[]) + { + int info; + CGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_factor(int m, MKL_Complex16 a[], int ipiv[]) + { + int info; + ZGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void s_lu_inverse(int n, float a[], int ipiv[], float work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + int info; + SGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_inverse(int n, double a[], int ipiv[], double work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + DGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + CGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + ZGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void s_lu_solve(int n, int nrhs, float a[], int ipiv[], float b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans = 'N'; + SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_solve(int n, int nrhs, double a[], int ipiv[], double b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans = 'N'; + DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans = 'N'; + CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans = 'N'; + ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void s_cholesky_solve(int n, int nrhs, float a[], float b[]) + { + char uplo = 'L'; + int info = 0; + SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void d_cholesky_solve(int n, int nrhs, double a[], double b[]) + { + char uplo = 'L'; + int info = 0; + DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) + { + char uplo = 'L'; + int info = 0; + CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) + { + char uplo = 'L'; + int info = 0; + ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len) + { + int info = 0; + SGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + SORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + SORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len) + { + int info = 0; + DGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + DORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + DORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], int len) + { + int info = 0; + CGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + CUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + CUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], int len) + { + int info = 0; + ZGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + ZUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + ZUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + cblas_strsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + cblas_dtrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + MKL_Complex8 one; + one.real = 1.0; + cblas_ctrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + MKL_Complex16 one; + one.real = 1.0; + cblas_ztrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + SGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + } + + DLLEXPORT void d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + DGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + } + + DLLEXPORT void c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], int len, float rwork[]) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + CGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); + } + + DLLEXPORT void z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[], MKL_Complex16 work[], int len, double rwork[]) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + ZGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); + } } \ No newline at end of file diff --git a/src/NativeWrappers/Windows/MKLWrapper32Tests/MKLWrapper32Tests.csproj b/src/NativeWrappers/Windows/MKLWrapper32Tests/MKLWrapper32Tests.csproj index 2c527cc0..fe0d7212 100644 --- a/src/NativeWrappers/Windows/MKLWrapper32Tests/MKLWrapper32Tests.csproj +++ b/src/NativeWrappers/Windows/MKLWrapper32Tests/MKLWrapper32Tests.csproj @@ -55,8 +55,9 @@ AllRules.ruleset - - ..\..\..\..\out\debug\Net40\MathNet.Numerics.dll + + False + ..\..\..\..\out\lib\Net40\MathNet.Numerics.dll diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index 53c88912..a5c9e7a8 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include @@ -154,7 +154,27 @@ /// This is equivalent to the GETRF LAPACK routine. public override void LUFactor(<#=dataType#>[] data, int order, int[] ipiv) { - throw new NotImplementedException(); + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + SafeNativeMethods.<#=prefix#>_lu_factor(order, data, ipiv); } /// @@ -284,6 +304,11 @@ throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + SafeNativeMethods.<#=prefix#>_cholesky_factor(order, a); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index 976a2d66..a5e2f389 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include @@ -101,6 +101,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> #region LAPACK + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern void s_cholesky_factor(int n, [In, Out] float[] a); @@ -112,17 +124,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern void z_cholesky_factor(int n, [In, Out] Complex[] a); - + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); + internal static extern void s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); + internal static extern void d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); + internal static extern void c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); + internal static extern void z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); #endregion LAPACK diff --git a/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs index bcfd4cec..b51e9ce5 100644 --- a/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs @@ -168,7 +168,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms double ho = (_points[i] - t) * d[i] / hp; double den = ho - c[i + 1]; - if (den.AlmostEqual(0.0)) + if (den == 0.0) { return double.NaN; // zero-div, singularity } diff --git a/src/Numerics/Permutation.cs b/src/Numerics/Permutation.cs index 5d821ac4..bc1abf14 100644 --- a/src/Numerics/Permutation.cs +++ b/src/Numerics/Permutation.cs @@ -110,7 +110,7 @@ namespace MathNet.Numerics /// From wikipedia: the permutation 12043 has the inversions (0,2), (1,2) and (3,4). This would be /// encoded using the array [22244]. /// - /// The set of inversions to contruct the permutation from. + /// The set of inversions to construct the permutation from. /// A permutation generated from a sequence of inversions. public static Permutation FromInversions(int[] inv) { @@ -176,7 +176,7 @@ namespace MathNet.Numerics /// /// An array which represents where each integer is permuted too: indices[i] represents that integer i /// is permuted to location indices[i]. - /// True if represents a proper permutation, false otherwise. + /// True if represents a proper permutation, false otherwise. static private bool CheckForProperPermutation(int[] indices) { var idxCheck = new bool[indices.Length];