From c6720b630f267bb274f5621182ae1953a4244d08 Mon Sep 17 00:00:00 2001 From: Kuan Bartel Date: Wed, 15 Jul 2015 01:37:07 +0900 Subject: [PATCH] Updated OpenBLAS LA provider to remove work arrays to match the MKL LA provider changes. Changed MKLWrapper and OpenBLASWrapper projects to use shared blas.c and lapack.cpp files. Each project has individual blas.h and lapack.h header files. --- src/NativeProviders/Common/blas.c | 89 +++ .../{MKL => Common}/lapack.cpp | 116 ++- src/NativeProviders/Common/lapack_common.h | 25 + src/NativeProviders/MKL/blas.c | 89 --- src/NativeProviders/MKL/blas.h | 7 + src/NativeProviders/MKL/lapack.h | 27 +- src/NativeProviders/OpenBLAS/blas.c | 89 --- src/NativeProviders/OpenBLAS/blas.h | 7 + src/NativeProviders/OpenBLAS/complex.h | 39 - src/NativeProviders/OpenBLAS/lapack.cpp | 712 ------------------ src/NativeProviders/OpenBLAS/lapack.h | 9 + .../Windows/MKL/MKLWrapper.vcxproj | 7 +- .../Windows/MKL/MKLWrapper.vcxproj.filters | 21 +- .../Windows/OpenBLAS/OpenBLASWrapper.vcxproj | 8 +- .../OpenBLAS/OpenBLASWrapper.vcxproj.filters | 18 +- src/Numerics/Control.cs | 2 + .../OpenBlasLinearAlgebraProvider.Complex.cs | 493 ++++-------- ...OpenBlasLinearAlgebraProvider.Complex32.cs | 486 ++++-------- .../OpenBlasLinearAlgebraProvider.Double.cs | 505 ++++--------- .../OpenBlasLinearAlgebraProvider.Single.cs | 500 ++++-------- .../OpenBlas/OpenBlasLinearAlgebraProvider.cs | 11 + .../OpenBlas/SafeNativeMethods.cs | 100 +-- 22 files changed, 860 insertions(+), 2500 deletions(-) create mode 100644 src/NativeProviders/Common/blas.c rename src/NativeProviders/{MKL => Common}/lapack.cpp (80%) delete mode 100644 src/NativeProviders/MKL/blas.c create mode 100644 src/NativeProviders/MKL/blas.h delete mode 100644 src/NativeProviders/OpenBLAS/blas.c create mode 100644 src/NativeProviders/OpenBLAS/blas.h delete mode 100644 src/NativeProviders/OpenBLAS/complex.h delete mode 100644 src/NativeProviders/OpenBLAS/lapack.cpp create mode 100644 src/NativeProviders/OpenBLAS/lapack.h diff --git a/src/NativeProviders/Common/blas.c b/src/NativeProviders/Common/blas.c new file mode 100644 index 00000000..4cbcef1f --- /dev/null +++ b/src/NativeProviders/Common/blas.c @@ -0,0 +1,89 @@ +#include "wrapper_common.h" +#include "blas.h" + +#if __cplusplus +extern "C" { +#endif +DLLEXPORT void s_axpy(const blas_int n, const float alpha, const float x[], float y[]){ + cblas_saxpy(n, alpha, x, 1, y, 1); +} + +DLLEXPORT void d_axpy(const blas_int n, const double alpha, const double x[], double y[]){ + cblas_daxpy(n, alpha, x, 1, y, 1); +} + +DLLEXPORT void c_axpy(const blas_int n, const blas_complex_float alpha, const blas_complex_float x[], blas_complex_float y[]){ + cblas_caxpy(n, (float*)&alpha, (float*)x, 1, (float*)y, 1); +} + +DLLEXPORT void z_axpy(const blas_int n, const blas_complex_double alpha, const blas_complex_double x[], blas_complex_double y[]){ + cblas_zaxpy(n, (double*)&alpha, (double*)x, 1, (double*)y, 1); +} + +DLLEXPORT void s_scale(const blas_int n, const float alpha, float x[]){ + cblas_sscal(n, alpha, x, 1); +} + +DLLEXPORT void d_scale(const blas_int n, const double alpha, double x[]){ + cblas_dscal(n, alpha, x, 1); +} + +DLLEXPORT void c_scale(const blas_int n, const blas_complex_float alpha, blas_complex_float x[]){ + cblas_cscal(n, (float*)&alpha, (float*)x, 1); +} + +DLLEXPORT void z_scale(const blas_int n, const blas_complex_double alpha, blas_complex_double x[]){ + cblas_zscal(n, (double*)&alpha, (double*)x, 1); +} + +DLLEXPORT float s_dot_product(const blas_int n, const float x[], const float y[]){ + return cblas_sdot(n, x, 1, y, 1); +} + +DLLEXPORT double d_dot_product(const blas_int n, const double x[], const double y[]){ + return cblas_ddot(n, x, 1, y, 1); +} + +DLLEXPORT blas_complex_float c_dot_product(const blas_int n, const blas_complex_float x[], const blas_complex_float y[]){ + blas_complex_float ret; + cblas_cdotu_sub(n, (float*)x, 1, (float*)y, 1, &ret); + return ret; +} + +DLLEXPORT blas_complex_double z_dot_product(const blas_int n, const blas_complex_double x[], const blas_complex_double y[]){ + blas_complex_double ret; + cblas_zdotu_sub(n, (double*)x, 1, (double*)y, 1, &ret); + return ret; +} + +DLLEXPORT void s_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blas_int m, const blas_int n, const blas_int k, const float alpha, const float x[], const float y[], const float beta, float c[]){ + blas_int lda = transA == CblasNoTrans ? m : k; + blas_int ldb = transB == CblasNoTrans ? k : n; + + cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} + +DLLEXPORT void d_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blas_int m, const blas_int n, const blas_int k, const double alpha, const double x[], const double y[], const double beta, double c[]){ + blas_int lda = transA == CblasNoTrans ? m : k; + blas_int ldb = transB == CblasNoTrans ? k : n; + + cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} + +DLLEXPORT void c_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blas_int m, const blas_int n, const blas_int k, const blas_complex_float alpha, const blas_complex_float x[], const blas_complex_float y[], const blas_complex_float beta, blas_complex_float c[]){ + blas_int lda = transA == CblasNoTrans ? m : k; + blas_int ldb = transB == CblasNoTrans ? k : n; + + cblas_cgemm(CblasColMajor, transA, transB, m, n, k, (float*)&alpha, (float*)x, lda, (float*)y, ldb, (float*)&beta, (float*)c, m); +} + +DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blas_int m, const blas_int n, const blas_int k, const blas_complex_double alpha, const blas_complex_double x[], const blas_complex_double y[], const blas_complex_double beta, blas_complex_double c[]){ + blas_int lda = transA == CblasNoTrans ? m : k; + blas_int ldb = transB == CblasNoTrans ? k : n; + + cblas_zgemm(CblasColMajor, transA, transB, m, n, k, (double*)&alpha, (double*)x, lda, (double*)y, ldb, (double*)&beta, (double*)c, m); +} + +#if __cplusplus +} +#endif diff --git a/src/NativeProviders/MKL/lapack.cpp b/src/NativeProviders/Common/lapack.cpp similarity index 80% rename from src/NativeProviders/MKL/lapack.cpp rename to src/NativeProviders/Common/lapack.cpp index aa500dae..b0955923 100644 --- a/src/NativeProviders/MKL/lapack.cpp +++ b/src/NativeProviders/Common/lapack.cpp @@ -1,35 +1,8 @@ -#include -#include #include "wrapper_common.h" #include "lapack.h" #include "lapack_common.h" - -struct ptrfree -{ - void operator()(void* x) { PTRFREE(x); } -}; - -template using ptr = std::unique_ptr < T, ptrfree >; - -template -inline ptr array_new(const int size, int alignment = ALIGNMENT) { - auto ret = static_cast(PTRALLOC(size * sizeof(T), alignment)); - - if (!ret) - { - throw new std::bad_alloc(); - } - - return ptr(ret); -} - -template -inline ptr array_clone(const int size, const T* array){ - auto clone = array_new(size); - memcpy(clone.get(), array, size * sizeof(T)); - return clone; -} +#include template inline lapack_int lu_factor(lapack_int m, T a[], lapack_int ipiv[], GETRF getrf) @@ -204,7 +177,8 @@ inline lapack_int qr_thin_factor(lapack_int m, lapack_int n, T q[], T tau[], T r template inline lapack_int qr_solve(lapack_int m, lapack_int n, lapack_int bn, T a[], T b[], T x[], GELS gels) { - try{ + try + { auto clone_a = array_clone(m * n, a); auto clone_b = array_clone(m * bn, b); auto info = gels(LAPACK_COL_MAJOR, 'N', m, n, bn, clone_a.get(), m, clone_b.get(), m); @@ -221,7 +195,6 @@ inline lapack_int qr_solve(lapack_int m, lapack_int n, lapack_int bn, T a[], T b { return INSUFFICIENT_MEMORY; } - } template @@ -245,10 +218,9 @@ inline lapack_int qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T { return INSUFFICIENT_MEMORY; } - } -template +template inline lapack_int complex_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], UNMQR unmqr, TRSM trsm) { try @@ -262,7 +234,7 @@ inline lapack_int complex_qr_solve_factored(lapack_int m, lapack_int n, lapack_i } T one = 1.0f; - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b.get(), m); + trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, reinterpret_cast(&one), reinterpret_cast(r), m, reinterpret_cast(clone_b.get()), m); copyBtoX(m, n, bn, clone_b.get(), x); return info; } @@ -363,7 +335,7 @@ inline lapack_int eigen_factor(lapack_int n, T a[], T vectors[], R values[], T d } template -inline lapack_int eigen_complex_factor(lapack_int n, T a[], T vectors[], std::complex values[], T d[], GEES gees, TREVC trevc) +inline lapack_int eigen_complex_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], GEES gees, TREVC trevc) { try { @@ -399,7 +371,7 @@ inline lapack_int eigen_complex_factor(lapack_int n, T a[], T vectors[], std::co } template -inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], std::complex values[], T d[], SYEV syev) +inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], SYEV syev) { try { @@ -416,7 +388,7 @@ inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], std::comple for (auto index = 0; index < n; ++index) { - values[index] = std::complex(w.get()[index]); + values[index] = lapack_complex_double(w.get()[index]); } for (auto j = 0; j < n; ++j) @@ -452,12 +424,12 @@ extern "C" { return LAPACKE_dlange(LAPACK_COL_MAJOR, norm, m, n, a, m); } - DLLEXPORT float c_matrix_norm(char norm, lapack_int m, lapack_int n, std::complex a[]) + DLLEXPORT float c_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_float a[]) { return LAPACKE_clange(LAPACK_COL_MAJOR, norm, m, n, a, m); } - DLLEXPORT double z_matrix_norm(char norm, lapack_int m, lapack_int n, std::complex a[]) + DLLEXPORT double z_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_double a[]) { return LAPACKE_zlange(LAPACK_COL_MAJOR, norm, m, n, a, m); } @@ -472,12 +444,12 @@ extern "C" { return lu_factor(m, a, ipiv, LAPACKE_dgetrf); } - DLLEXPORT lapack_int c_lu_factor(lapack_int m, std::complex a[], lapack_int ipiv[]) + DLLEXPORT lapack_int c_lu_factor(lapack_int m, lapack_complex_float a[], lapack_int ipiv[]) { return lu_factor(m, a, ipiv, LAPACKE_cgetrf); } - DLLEXPORT lapack_int z_lu_factor(lapack_int m, std::complex a[], lapack_int ipiv[]) + DLLEXPORT lapack_int z_lu_factor(lapack_int m, lapack_complex_double a[], lapack_int ipiv[]) { return lu_factor(m, a, ipiv, LAPACKE_zgetrf); } @@ -492,12 +464,12 @@ extern "C" { return lu_inverse(n, a, LAPACKE_dgetrf, LAPACKE_dgetri); } - DLLEXPORT lapack_int c_lu_inverse(lapack_int n, std::complex a[], std::complex work[], lapack_int lwork) + DLLEXPORT lapack_int c_lu_inverse(lapack_int n, lapack_complex_float a[], lapack_complex_float work[], lapack_int lwork) { return lu_inverse(n, a, LAPACKE_cgetrf, LAPACKE_cgetri); } - DLLEXPORT lapack_int z_lu_inverse(lapack_int n, std::complex a[], std::complex work[], lapack_int lwork) + DLLEXPORT lapack_int z_lu_inverse(lapack_int n, lapack_complex_double a[], lapack_complex_double work[], lapack_int lwork) { return lu_inverse(n, a, LAPACKE_zgetrf, LAPACKE_zgetri); } @@ -512,12 +484,12 @@ extern "C" { return lu_inverse_factored(n, a, ipiv, LAPACKE_dgetri); } - DLLEXPORT lapack_int c_lu_inverse_factored(lapack_int n, std::complex a[], lapack_int ipiv[], std::complex work[], lapack_int lwork) + DLLEXPORT lapack_int c_lu_inverse_factored(lapack_int n, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float work[], lapack_int lwork) { return lu_inverse_factored(n, a, ipiv, LAPACKE_cgetri); } - DLLEXPORT lapack_int z_lu_inverse_factored(lapack_int n, std::complex a[], lapack_int ipiv[], std::complex work[], lapack_int lwork) + DLLEXPORT lapack_int z_lu_inverse_factored(lapack_int n, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double work[], lapack_int lwork) { return lu_inverse_factored(n, a, ipiv, LAPACKE_zgetri); } @@ -532,12 +504,12 @@ extern "C" { return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACKE_dgetrs); } - DLLEXPORT lapack_int c_lu_solve_factored(lapack_int n, lapack_int nrhs, std::complex a[], lapack_int ipiv[], std::complex b[]) + DLLEXPORT lapack_int c_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float b[]) { return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACKE_cgetrs); } - DLLEXPORT lapack_int z_lu_solve_factored(lapack_int n, lapack_int nrhs, std::complex a[], lapack_int ipiv[], std::complex b[]) + DLLEXPORT lapack_int z_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double b[]) { return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACKE_zgetrs); } @@ -552,12 +524,12 @@ extern "C" { return lu_solve(n, nrhs, a, b, LAPACKE_dgetrf, LAPACKE_dgetrs); } - DLLEXPORT lapack_int c_lu_solve(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int c_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) { return lu_solve(n, nrhs, a, b, LAPACKE_cgetrf, LAPACKE_cgetrs); } - DLLEXPORT lapack_int z_lu_solve(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int z_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) { return lu_solve(n, nrhs, a, b, LAPACKE_zgetrf, LAPACKE_zgetrs); } @@ -572,12 +544,12 @@ extern "C" { return cholesky_factor(n, a, LAPACKE_dpotrf); } - DLLEXPORT lapack_int c_cholesky_factor(lapack_int n, std::complex a[]) + DLLEXPORT lapack_int c_cholesky_factor(lapack_int n, lapack_complex_float a[]) { return cholesky_factor(n, a, LAPACKE_cpotrf); } - DLLEXPORT lapack_int z_cholesky_factor(lapack_int n, std::complex a[]) + DLLEXPORT lapack_int z_cholesky_factor(lapack_int n, lapack_complex_double a[]) { return cholesky_factor(n, a, LAPACKE_zpotrf); } @@ -592,12 +564,12 @@ extern "C" { return cholesky_solve(n, nrhs, a, b, LAPACKE_dpotrf, LAPACKE_dpotrs); } - DLLEXPORT lapack_int c_cholesky_solve(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int c_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) { return cholesky_solve(n, nrhs, a, b, LAPACKE_cpotrf, LAPACKE_cpotrs); } - DLLEXPORT lapack_int z_cholesky_solve(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int z_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) { return cholesky_solve(n, nrhs, a, b, LAPACKE_zpotrf, LAPACKE_zpotrs); } @@ -612,12 +584,12 @@ extern "C" { return LAPACKE_dpotrs(LAPACK_COL_MAJOR, 'L', n, nrhs, a, n, b, n); } - DLLEXPORT lapack_int c_cholesky_solve_factored(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int c_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) { return LAPACKE_cpotrs(LAPACK_COL_MAJOR, 'L', n, nrhs, a, n, b, n); } - DLLEXPORT lapack_int z_cholesky_solve_factored(lapack_int n, lapack_int nrhs, std::complex a[], std::complex b[]) + DLLEXPORT lapack_int z_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) { return LAPACKE_zpotrs(LAPACK_COL_MAJOR, 'L', n, nrhs, a, n, b, n); } @@ -642,22 +614,22 @@ extern "C" { return qr_thin_factor(m, n, q, tau, r, LAPACKE_dgeqrf, LAPACKE_dorgqr); } - DLLEXPORT lapack_int c_qr_factor(lapack_int m, lapack_int n, std::complex r[], std::complex tau[], std::complex q[]) + DLLEXPORT lapack_int c_qr_factor(lapack_int m, lapack_int n, lapack_complex_float r[], lapack_complex_float tau[], lapack_complex_float q[]) { return qr_factor(m, n, r, tau, q, LAPACKE_cgeqrf, LAPACKE_cungqr); } - DLLEXPORT lapack_int c_qr_thin_factor(lapack_int m, lapack_int n, std::complex q[], std::complex tau[], std::complex r[]) + DLLEXPORT lapack_int c_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_float q[], lapack_complex_float tau[], lapack_complex_float r[]) { return qr_thin_factor(m, n, q, tau, r, LAPACKE_cgeqrf, LAPACKE_cungqr); } - DLLEXPORT lapack_int z_qr_factor(lapack_int m, lapack_int n, std::complex r[], std::complex tau[], std::complex q[]) + DLLEXPORT lapack_int z_qr_factor(lapack_int m, lapack_int n, lapack_complex_double r[], lapack_complex_double tau[], lapack_complex_double q[]) { return qr_factor(m, n, r, tau, q, LAPACKE_zgeqrf, LAPACKE_zungqr); } - DLLEXPORT lapack_int z_qr_thin_factor(lapack_int m, lapack_int n, std::complex q[], std::complex tau[], std::complex r[]) + DLLEXPORT lapack_int z_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_double q[], lapack_complex_double tau[], lapack_complex_double r[]) { return qr_thin_factor(m, n, q, tau, r, LAPACKE_zgeqrf, LAPACKE_zungqr); } @@ -672,12 +644,12 @@ extern "C" { return qr_solve(m, n, bn, a, b, x, LAPACKE_dgels); } - DLLEXPORT lapack_int c_qr_solve(lapack_int m, lapack_int n, lapack_int bn, std::complex a[], std::complex b[], std::complex x[]) + DLLEXPORT lapack_int c_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float a[], lapack_complex_float b[], lapack_complex_float x[]) { return qr_solve(m, n, bn, a, b, x, LAPACKE_cgels); } - DLLEXPORT lapack_int z_qr_solve(lapack_int m, lapack_int n, lapack_int bn, std::complex a[], std::complex b[], std::complex x[]) + DLLEXPORT lapack_int z_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double a[], lapack_complex_double b[], lapack_complex_double x[]) { return qr_solve(m, n, bn, a, b, x, LAPACKE_zgels); } @@ -692,14 +664,14 @@ extern "C" { return qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_dormqr, cblas_dtrsm); } - DLLEXPORT lapack_int c_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, std::complex r[], std::complex b[], std::complex tau[], std::complex x[]) + DLLEXPORT lapack_int c_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float r[], lapack_complex_float b[], lapack_complex_float tau[], lapack_complex_float x[]) { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_cunmqr, cblas_ctrsm); + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_cunmqr, cblas_ctrsm); } - DLLEXPORT lapack_int z_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, std::complex r[], std::complex b[], std::complex tau[], std::complex x[]) + DLLEXPORT lapack_int z_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double r[], lapack_complex_double b[], lapack_complex_double tau[], lapack_complex_double x[]) { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_zunmqr, cblas_ztrsm); + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_zunmqr, cblas_ztrsm); } DLLEXPORT lapack_int s_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, float a[], float s[], float u[], float v[]) @@ -712,17 +684,17 @@ extern "C" { return svd_factor(compute_vectors, m, n, a, s, u, v, LAPACKE_dgesvd); } - DLLEXPORT lapack_int c_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, std::complex a[], std::complex s[], std::complex u[], std::complex v[]) + DLLEXPORT lapack_int c_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_float a[], lapack_complex_float s[], lapack_complex_float u[], lapack_complex_float v[]) { - return complex_svd_factor, float>(compute_vectors, m, n, a, s, u, v, LAPACKE_cgesvd); + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, LAPACKE_cgesvd); } - DLLEXPORT lapack_int z_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, std::complex a[], std::complex s[], std::complex u[], std::complex v[]) + DLLEXPORT lapack_int z_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_double a[], lapack_complex_double s[], lapack_complex_double u[], lapack_complex_double v[]) { - return complex_svd_factor, double>(compute_vectors, m, n, a, s, u, v, LAPACKE_zgesvd); + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, LAPACKE_zgesvd); } - DLLEXPORT lapack_int s_eigen(bool isSymmetric, lapack_int n, float a[], float vectors[], std::complex values[], float d[]) + DLLEXPORT lapack_int s_eigen(bool isSymmetric, lapack_int n, float a[], float vectors[], lapack_complex_double values[], float d[]) { if (isSymmetric) { @@ -734,7 +706,7 @@ extern "C" { } } - DLLEXPORT lapack_int d_eigen(bool isSymmetric, lapack_int n, double a[], double vectors[], std::complex values[], double d[]) + DLLEXPORT lapack_int d_eigen(bool isSymmetric, lapack_int n, double a[], double vectors[], lapack_complex_double values[], double d[]) { if (isSymmetric) { @@ -746,7 +718,7 @@ extern "C" { } } - DLLEXPORT lapack_int c_eigen(bool isSymmetric, lapack_int n, std::complex a[], std::complex vectors[], std::complex values[], std::complex d[]) + DLLEXPORT lapack_int c_eigen(bool isSymmetric, lapack_int n, lapack_complex_float a[], lapack_complex_float vectors[], lapack_complex_double values[], lapack_complex_float d[]) { if (isSymmetric) { @@ -758,7 +730,7 @@ extern "C" { } } - DLLEXPORT lapack_int z_eigen(bool isSymmetric, lapack_int n, std::complex a[], std::complex vectors[], std::complex values[], std::complex d[]) + DLLEXPORT lapack_int z_eigen(bool isSymmetric, lapack_int n, lapack_complex_double a[], lapack_complex_double vectors[], lapack_complex_double values[], lapack_complex_double d[]) { if (isSymmetric) { diff --git a/src/NativeProviders/Common/lapack_common.h b/src/NativeProviders/Common/lapack_common.h index 46677a4b..5441d599 100644 --- a/src/NativeProviders/Common/lapack_common.h +++ b/src/NativeProviders/Common/lapack_common.h @@ -35,4 +35,29 @@ inline void copyBtoX (int m, int n, int bn, T b[], T x[]){ } } +#ifndef LAPACK_MEMORY +#define LAPACK_MEMORY + +#include + +const int INSUFFICIENT_MEMORY = -999999; +const int ALIGNMENT = 64; + +template using array_ptr = std::unique_ptr; + +template +inline array_ptr array_new(const int size, int alignment = ALIGNMENT) +{ + return array_ptr(new T[size]); +} + #endif + +template +inline array_ptr array_clone(const int size, const T* array) { + auto clone = array_new(size); + memcpy(clone.get(), array, size * sizeof(T)); + return clone; +} + +#endif \ No newline at end of file diff --git a/src/NativeProviders/MKL/blas.c b/src/NativeProviders/MKL/blas.c deleted file mode 100644 index 45fa4c55..00000000 --- a/src/NativeProviders/MKL/blas.c +++ /dev/null @@ -1,89 +0,0 @@ -#include "mkl_cblas.h" -#include "wrapper_common.h" - -#if __cplusplus -extern "C" { -#endif -DLLEXPORT void s_axpy(const MKL_INT n, const float alpha, const float x[], float y[]){ - cblas_saxpy(n, alpha, x, 1, y, 1); -} - -DLLEXPORT void d_axpy(const MKL_INT n, const double alpha, const double x[], double y[]){ - cblas_daxpy(n, alpha, x, 1, y, 1); -} - -DLLEXPORT void c_axpy(const MKL_INT n, const MKL_Complex8 alpha, const MKL_Complex8 x[], MKL_Complex8 y[]){ - cblas_caxpy(n, &alpha, x, 1, y, 1); -} - -DLLEXPORT void z_axpy(const MKL_INT n, const MKL_Complex16 alpha, const MKL_Complex16 x[], MKL_Complex16 y[]){ - cblas_zaxpy(n, &alpha, x, 1, y, 1); -} - -DLLEXPORT void s_scale(const MKL_INT n, const float alpha, float x[]){ - cblas_sscal(n, alpha, x, 1); -} - -DLLEXPORT void d_scale(const MKL_INT n, const double alpha, double x[]){ - cblas_dscal(n, alpha, x, 1); -} - -DLLEXPORT void c_scale(const MKL_INT n, const MKL_Complex8 alpha, MKL_Complex8 x[]){ - cblas_cscal(n, &alpha, x, 1); -} - -DLLEXPORT void z_scale(const MKL_INT n, const MKL_Complex16 alpha, MKL_Complex16 x[]){ - cblas_zscal(n, &alpha, x, 1); -} - -DLLEXPORT float s_dot_product(const MKL_INT n, const float x[], const float y[]){ - return cblas_sdot(n, x, 1, y, 1); -} - -DLLEXPORT double d_dot_product(const MKL_INT n, const double x[], const double y[]){ - return cblas_ddot(n, x, 1, y, 1); -} - -DLLEXPORT MKL_Complex8 c_dot_product(const MKL_INT n, const MKL_Complex8 x[], const MKL_Complex8 y[]){ - MKL_Complex8 ret; - cblas_cdotu_sub(n, x, 1, y, 1, &ret); - return ret; -} - -DLLEXPORT MKL_Complex16 z_dot_product(const MKL_INT n, const MKL_Complex16 x[], const MKL_Complex16 y[]){ - MKL_Complex16 ret; - cblas_zdotu_sub(n, x, 1, y, 1, &ret); - return ret; -} - -DLLEXPORT void s_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const MKL_INT m, const MKL_INT n, const MKL_INT k, const float alpha, const float x[], const float y[], const float beta, float c[]){ - MKL_INT lda = transA == CblasNoTrans ? m : k; - MKL_INT ldb = transB == CblasNoTrans ? k : n; - - cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); -} - -DLLEXPORT void d_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const MKL_INT m, const MKL_INT n, const MKL_INT k, const double alpha, const double x[], const double y[], const double beta, double c[]){ - MKL_INT lda = transA == CblasNoTrans ? m : k; - MKL_INT ldb = transB == CblasNoTrans ? k : n; - - cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); -} - -DLLEXPORT void c_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const MKL_INT m, const MKL_INT n, const MKL_INT k, const MKL_Complex8 alpha, const MKL_Complex8 x[], const MKL_Complex8 y[], const MKL_Complex8 beta, MKL_Complex8 c[]){ - MKL_INT lda = transA == CblasNoTrans ? m : k; - MKL_INT ldb = transB == CblasNoTrans ? k : n; - - cblas_cgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); -} - -DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const MKL_INT m, const MKL_INT n, const MKL_INT k, const MKL_Complex16 alpha, const MKL_Complex16 x[], const MKL_Complex16 y[], const MKL_Complex16 beta, MKL_Complex16 c[]){ - MKL_INT lda = transA == CblasNoTrans ? m : k; - MKL_INT ldb = transB == CblasNoTrans ? k : n; - - cblas_zgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); -} - -#if __cplusplus -} -#endif diff --git a/src/NativeProviders/MKL/blas.h b/src/NativeProviders/MKL/blas.h new file mode 100644 index 00000000..555d466c --- /dev/null +++ b/src/NativeProviders/MKL/blas.h @@ -0,0 +1,7 @@ +#pragma once + +#include "mkl_cblas.h" + +#define blas_int MKL_INT +#define blas_complex_float MKL_Complex8 +#define blas_complex_double MKL_Complex16 diff --git a/src/NativeProviders/MKL/lapack.h b/src/NativeProviders/MKL/lapack.h index dce0c23a..b57ca9ee 100644 --- a/src/NativeProviders/MKL/lapack.h +++ b/src/NativeProviders/MKL/lapack.h @@ -1,17 +1,34 @@ #pragma once -#include +#include #define MKL_Complex8 std::complex #define MKL_Complex16 std::complex #include "mkl.h" +#define LAPACK_MEMORY +#include + const int INSUFFICIENT_MEMORY = -999999; const int ALIGNMENT = 64; -//#define PTRALLOC( size, alignment ) LAPACKE_malloc( size ) -//#define PTRFREE( p ) LAPACKE_free( p ) +struct array_free +{ + void operator()(void* x) { mkl_free(x); } +}; + +template using array_ptr = std::unique_ptr; + +template +inline array_ptr array_new(const int size, int alignment = ALIGNMENT) +{ + auto ret = static_cast(mkl_malloc(size * sizeof(T), alignment)); + + if (!ret) + { + throw new std::bad_alloc(); + } -#define PTRALLOC( size, alignment ) mkl_malloc( size, alignment ) -#define PTRFREE( p ) mkl_free( p ) + return array_ptr(ret); +} diff --git a/src/NativeProviders/OpenBLAS/blas.c b/src/NativeProviders/OpenBLAS/blas.c deleted file mode 100644 index 3bf9290e..00000000 --- a/src/NativeProviders/OpenBLAS/blas.c +++ /dev/null @@ -1,89 +0,0 @@ -#include "cblas.h" -#include "wrapper_common.h" - -#if __cplusplus -extern "C" { -#endif -DLLEXPORT void s_axpy(const blasint n, const float alpha, const float x[], float y[]){ - cblas_saxpy(n, alpha, x, 1, y, 1); -} - -DLLEXPORT void d_axpy(const blasint n, const double alpha, const double x[], double y[]){ - cblas_daxpy(n, alpha, x, 1, y, 1); -} - -DLLEXPORT void c_axpy(const blasint n, const openblas_complex_float alpha, const openblas_complex_float x[], openblas_complex_float y[]){ - cblas_caxpy(n, (float*)&alpha, (float*)x, 1, (float*)y, 1); -} - -DLLEXPORT void z_axpy(const blasint n, const openblas_complex_double alpha, const openblas_complex_double x[], openblas_complex_double y[]){ - cblas_zaxpy(n, (double*)&alpha, (double*)x, 1, (double*)y, 1); -} - -DLLEXPORT void s_scale(const blasint n, const float alpha, float x[]){ - cblas_sscal(n, alpha, x, 1); -} - -DLLEXPORT void d_scale(const blasint n, const double alpha, double x[]){ - cblas_dscal(n, alpha, x, 1); -} - -DLLEXPORT void c_scale(const blasint n, const openblas_complex_float alpha, openblas_complex_float x[]){ - cblas_cscal(n, (float*)&alpha, (float*)x, 1); -} - -DLLEXPORT void z_scale(const blasint n, const openblas_complex_double alpha, openblas_complex_double x[]){ - cblas_zscal(n, (double*)&alpha, (double*)x, 1); -} - -DLLEXPORT float s_dot_product(const blasint n, const float x[], const float y[]){ - return cblas_sdot(n, x, 1, y, 1); -} - -DLLEXPORT double d_dot_product(const blasint n, const double x[], const double y[]){ - return cblas_ddot(n, x, 1, y, 1); -} - -DLLEXPORT openblas_complex_float c_dot_product(const blasint n, const openblas_complex_float x[], const openblas_complex_float y[]){ - openblas_complex_float ret; - cblas_cdotu_sub(n, (float*)x, 1, (float*)y, 1, &ret); - return ret; -} - -DLLEXPORT openblas_complex_double z_dot_product(const blasint n, const openblas_complex_double x[], const openblas_complex_double y[]){ - openblas_complex_double ret; - cblas_zdotu_sub(n, (double*)x, 1, (double*)y, 1, &ret); - return ret; -} - -DLLEXPORT void s_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const float alpha, const float x[], const float y[], const float beta, float c[]){ - blasint lda = transA == CblasNoTrans ? m : k; - blasint ldb = transB == CblasNoTrans ? k : n; - - cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); -} - -DLLEXPORT void d_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const double alpha, const double x[], const double y[], const double beta, double c[]){ - blasint lda = transA == CblasNoTrans ? m : k; - blasint ldb = transB == CblasNoTrans ? k : n; - - cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); -} - -DLLEXPORT void c_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const openblas_complex_float alpha, const openblas_complex_float x[], const openblas_complex_float y[], const openblas_complex_float beta, openblas_complex_float c[]){ - blasint lda = transA == CblasNoTrans ? m : k; - blasint ldb = transB == CblasNoTrans ? k : n; - - cblas_cgemm(CblasColMajor, transA, transB, m, n, k, (float*)&alpha, (float*)x, lda, (float*)y, ldb, (float*)&beta, (float*)c, m); -} - -DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const openblas_complex_double alpha, const openblas_complex_double x[], const openblas_complex_double y[], const openblas_complex_double beta, openblas_complex_double c[]){ - blasint lda = transA == CblasNoTrans ? m : k; - blasint ldb = transB == CblasNoTrans ? k : n; - - cblas_zgemm(CblasColMajor, transA, transB, m, n, k, (double*)&alpha, (double*)x, lda, (double*)y, ldb, (double*)&beta, (double*)c, m); -} - -#if __cplusplus -} -#endif diff --git a/src/NativeProviders/OpenBLAS/blas.h b/src/NativeProviders/OpenBLAS/blas.h new file mode 100644 index 00000000..c5c0a7e0 --- /dev/null +++ b/src/NativeProviders/OpenBLAS/blas.h @@ -0,0 +1,7 @@ +#pragma once + +#include "cblas.h" + +#define blas_int blasint +#define blas_complex_float openblas_complex_float +#define blas_complex_double openblas_complex_double diff --git a/src/NativeProviders/OpenBLAS/complex.h b/src/NativeProviders/OpenBLAS/complex.h deleted file mode 100644 index d107e400..00000000 --- a/src/NativeProviders/OpenBLAS/complex.h +++ /dev/null @@ -1,39 +0,0 @@ -template -struct complex -{ - _T real, imag; - - complex(_T _real = 0, _T _imag = 0) - { - real = _real; - imag = _imag; - } - - complex(const complex<_T>& right) - { - real = right.real; - imag = right.imag; - } - - complex& operator=(const complex& right) - { - real = right.real; - imag = right.imag; - return *this; - } - - complex& operator=(const _T& right) - { - real = right; - imag = 0; - return *this; - } - - template inline - complex& operator=(const complex<_Other>& right) - { - real = (_T)right.real; - imag = (_T)right.imag; - return *this; - } -}; \ No newline at end of file diff --git a/src/NativeProviders/OpenBLAS/lapack.cpp b/src/NativeProviders/OpenBLAS/lapack.cpp deleted file mode 100644 index fee9abec..00000000 --- a/src/NativeProviders/OpenBLAS/lapack.cpp +++ /dev/null @@ -1,712 +0,0 @@ -#include "cblas.h" - -#include "complex.h" -#define LAPACK_COMPLEX_CUSTOM -#define lapack_complex_float complex -#define lapack_complex_double complex - -#include "lapacke.h" -#include "lapack_common.h" -#include "wrapper_common.h" -#include - -template -inline lapack_int lu_factor(lapack_int m, T a[], lapack_int ipiv[], GETRF getrf) -{ - lapack_int info = 0; - getrf(&m, &m, a, &m, ipiv, &info); - shift_ipiv_down(m, ipiv); - return info; -}; - -template -inline lapack_int lu_inverse(lapack_int n, T a[], T work[], lapack_int lwork, GETRF getrf, GETRI getri) -{ - lapack_int* ipiv = new lapack_int[n]; - lapack_int info = 0; - getrf(&n, &n, a, &n, ipiv, &info); - - if (info != 0) - { - delete[] ipiv; - return info; - } - - getri(&n, a, &n, ipiv, work, &lwork, &info); - delete[] ipiv; - return info; -}; - -template -inline lapack_int lu_inverse_factored(lapack_int n, T a[], lapack_int ipiv[], T work[], lapack_int lwork, GETRI getri) -{ - shift_ipiv_up(n, ipiv); - lapack_int info = 0; - getri(&n, a, &n, ipiv, work, &lwork, &info); - shift_ipiv_down(n, ipiv); - return info; -} - -template -inline lapack_int lu_solve_factored(lapack_int n, lapack_int nrhs, T a[], lapack_int ipiv[], T b[], GETRS getrs) -{ - shift_ipiv_up(n, ipiv); - lapack_int info = 0; - char trans ='N'; - getrs(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - shift_ipiv_down(n, ipiv); - return info; -} - -template -inline lapack_int lu_solve(lapack_int n, lapack_int nrhs, T a[], T b[], GETRF getrf, GETRS getrs) -{ - T* clone = Clone(n, n, a); - lapack_int* ipiv = new lapack_int[n]; - lapack_int info = 0; - getrf(&n, &n, clone, &n, ipiv, &info); - - if (info != 0) - { - delete[] ipiv; - delete[] clone; - return info; - } - - char trans ='N'; - getrs(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); - delete[] ipiv; - delete[] clone; - return info; -} - - -template -inline lapack_int cholesky_factor(lapack_int n, T* a, POTRF potrf) -{ - char uplo = 'L'; - lapack_int info = 0; - potrf(&uplo, &n, a, &n, &info); - T zero = T(); - - for (lapack_int i = 0; i < n; ++i) - { - lapack_int index = i * n; - - for (lapack_int j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - - return info; -} - -template -inline lapack_int cholesky_solve(lapack_int n, lapack_int nrhs, T a[], T b[], POTRF potrf, POTRS potrs) -{ - T* clone = Clone(n, n, a); - char uplo = 'L'; - lapack_int info = 0; - potrf(&uplo, &n, clone, &n, &info); - - if (info != 0) - { - delete[] clone; - return info; - } - - potrs(&uplo, &n, &nrhs, clone, &n, b, &n, &info); - delete[] clone; - return info; -} - -template -inline lapack_int cholesky_solve_factored(lapack_int n, lapack_int nrhs, T a[], T b[], POTRS potrs) -{ - char uplo = 'L'; - lapack_int info = 0; - potrs(&uplo, &n, &nrhs, a, &n, b, &n, &info); - return info; -} - -template -inline lapack_int qr_factor(lapack_int m, lapack_int n, T r[], T tau[], T q[], T work[], lapack_int len, GEQRF geqrf, ORGQR orgqr) -{ - lapack_int info = 0; - geqrf(&m, &n, r, &m, tau, work, &len, &info); - - for (lapack_int i = 0; i < m; ++i) - { - for (lapack_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) - { - orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); - } - - return info; -} - -template -inline lapack_int qr_thin_factor(lapack_int m, lapack_int n, T q[], T tau[], T r[], T work[], lapack_int len, GEQRF geqrf, ORGQR orgqr) -{ - lapack_int info = 0; - geqrf(&m, &n, q, &m, tau, work, &len, &info); - - for (lapack_int i = 0; i < n; ++i) - { - for (lapack_int j = 0; j < n; ++j) - { - if (i <= j) - { - r[j * n + i] = q[j * m + i]; - } - } - } - - orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); - return info; -} - -template -inline lapack_int qr_solve(lapack_int m, lapack_int n, lapack_int bn, T a[], T b[], T x[], T work[], lapack_int len, GELS gels) -{ - T* clone_a = Clone(m, n, a); - T* clone_b = Clone(m, bn, b); - char N = 'N'; - lapack_int info = 0; - gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - copyBtoX(m, n, bn, clone_b, x); - delete[] clone_a; - delete[] clone_b; - return info; -} - -template -inline lapack_int qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], T work[], lapack_int len, ORMQR ormqr, TRSM trsm) -{ - T* clone_b = Clone(m, bn, b); - char side ='L'; - char tran = 'T'; - lapack_int info = 0; - ormqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); - copyBtoX(m, n, bn, clone_b, x); - delete[] clone_b; - return info; -} - -template -inline lapack_int complex_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], T work[], lapack_int len, UNMQR unmqr, TRSM trsm) -{ - T* clone_b = Clone(m, bn, b); - char side ='L'; - char tran = 'C'; - lapack_int info = 0; - unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - T one = { 1.0f, 0.0f }; - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, reinterpret_cast(&one), reinterpret_cast(r), m, reinterpret_cast(clone_b), m); - copyBtoX(m, n, bn, clone_b, x); - delete[] clone_b; - return info; -} - -template -inline lapack_int svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], T work[], lapack_int len, GESVD gesvd) -{ - lapack_int info = 0; - char job = compute_vectors ? 'A' : 'N'; - gesvd(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); - return info; -} - -template -inline lapack_int complex_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], T work[], lapack_int len, GESVD gesvd) -{ - lapack_int info = 0; - lapack_int dim_s = std::min(m,n); - R* rwork = new R[5 * dim_s]; - R* s_local = new R[dim_s]; - char job = compute_vectors ? 'A' : 'N'; - gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info); - - for (lapack_int index = 0; index < dim_s; ++index) - { - s[index] = s_local[index]; - } - - delete[] rwork; - delete[] s_local; - return info; -} - -template -inline lapack_int eigen_factor(lapack_int n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) -{ - T* clone_a = Clone(n, n, a); - T* wr = new T[n]; - T* wi = new T[n]; - - lapack_int sdim; - lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, wr, wi, vectors, n); - if (info != 0) - { - delete[] clone_a; - delete[] wr; - delete[] wi; - return info; - } - - lapack_int m; - info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); - if (info != 0) - { - delete[] clone_a; - delete[] wr; - delete[] wi; - return info; - } - - for (lapack_int index = 0; index < n; ++index) - { - values[index] = R(wr[index], wi[index]); - } - - for (lapack_int i = 0; i < n; ++i) - { - lapack_int in = i * n; - d[in + i] = wr[i]; - - if (wi[i] > 0) - { - d[in + n + i] = wi[i]; - } - else if (wi[i] < 0) - { - d[in - n + i] = wi[i]; - } - } - - delete[] clone_a; - delete[] wr; - delete[] wi; - return info; -} - -template -inline lapack_int eigen_complex_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], GEES gees, TREVC trevc) -{ - T* clone_a = Clone(n, n, a); - T* w = new T[n]; - - lapack_int sdim; - lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, w, vectors, n); - if (info != 0) - { - delete[] clone_a; - delete[] w; - return info; - } - - lapack_int m; - info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); - if (info != 0) - { - delete[] clone_a; - delete[] w; - return info; - } - - for (lapack_int i = 0; i < n; ++i) - { - values[i] = w[i]; - d[i * n + i] = w[i]; - } - - delete[] clone_a; - delete[] w; - return info; -} - -template -inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], SYEV syev) -{ - T* clone_a = Clone(n, n, a); - R* w = new R[n]; - - lapack_int info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w); - if (info != 0) - { - delete[] clone_a; - delete[] w; - return info; - } - - memcpy(vectors, clone_a, n*n*sizeof(T)); - - for (lapack_int index = 0; index < n; ++index) - { - values[index] = lapack_complex_double(w[index]); - } - - for (lapack_int j = 0; j < n; ++j) - { - lapack_int jn = j*n; - - for (lapack_int i = 0; i < n; ++i) - { - if (i == j) - { - d[jn + i] = w[i]; - } - } - } - - delete[] clone_a; - delete[] w; - return info; -} - -extern "C" { - - DLLEXPORT float s_matrix_norm(char norm, lapack_int m, lapack_int n, float a[], float work[]) - { - return LAPACKE_slange_work(CblasColMajor, norm, m, n, a, m, work); - } - - DLLEXPORT double d_matrix_norm(char norm, lapack_int m, lapack_int n, double a[], double work[]) - { - return LAPACKE_dlange_work(CblasColMajor, norm, m, n, a, m, work); - } - - DLLEXPORT float c_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_float a[], float work[]) - { - return LAPACKE_clange_work(CblasColMajor, norm, m, n, a, m, work); - } - - DLLEXPORT double z_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_double a[], double work[]) - { - return LAPACKE_zlange_work(CblasColMajor, norm, m, n, a, m, work); - } - - DLLEXPORT lapack_int s_lu_factor(lapack_int m, float a[], lapack_int ipiv[]) - { - return lu_factor(m, a, ipiv, LAPACK_sgetrf); - } - - DLLEXPORT lapack_int d_lu_factor(lapack_int m, double a[], lapack_int ipiv[]) - { - return lu_factor(m, a, ipiv, LAPACK_dgetrf); - } - - DLLEXPORT lapack_int c_lu_factor(lapack_int m, lapack_complex_float a[], lapack_int ipiv[]) - { - return lu_factor(m, a, ipiv, LAPACK_cgetrf); - } - - DLLEXPORT lapack_int z_lu_factor(lapack_int m, lapack_complex_double a[], lapack_int ipiv[]) - { - return lu_factor(m, a, ipiv, LAPACK_zgetrf); - } - - DLLEXPORT lapack_int s_lu_inverse(lapack_int n, float a[], float work[], lapack_int lwork) - { - return lu_inverse(n, a, work, lwork, LAPACK_sgetrf, LAPACK_sgetri); - } - - DLLEXPORT lapack_int d_lu_inverse(lapack_int n, double a[], double work[], lapack_int lwork) - { - return lu_inverse(n, a, work, lwork, LAPACK_dgetrf, LAPACK_dgetri); - } - - DLLEXPORT lapack_int c_lu_inverse(lapack_int n, lapack_complex_float a[], lapack_complex_float work[], lapack_int lwork) - { - return lu_inverse(n, a, work, lwork, LAPACK_cgetrf, LAPACK_cgetri); - } - - DLLEXPORT lapack_int z_lu_inverse(lapack_int n, lapack_complex_double a[], lapack_complex_double work[], lapack_int lwork) - { - return lu_inverse(n, a, work, lwork, LAPACK_zgetrf, LAPACK_zgetri); - } - - DLLEXPORT lapack_int s_lu_inverse_factored(lapack_int n, float a[], lapack_int ipiv[], float work[], lapack_int lwork) - { - return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_sgetri); - } - - DLLEXPORT lapack_int d_lu_inverse_factored(lapack_int n, double a[], lapack_int ipiv[], double work[], lapack_int lwork) - { - return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_dgetri); - } - - DLLEXPORT lapack_int c_lu_inverse_factored(lapack_int n, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float work[], lapack_int lwork) - { - return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_cgetri); - } - - DLLEXPORT lapack_int z_lu_inverse_factored(lapack_int n, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double work[], lapack_int lwork) - { - return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_zgetri); - } - - DLLEXPORT lapack_int s_lu_solve_factored(lapack_int n, lapack_int nrhs, float a[], lapack_int ipiv[], float b[]) - { - return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_sgetrs); - } - - DLLEXPORT lapack_int d_lu_solve_factored(lapack_int n, lapack_int nrhs, double a[], lapack_int ipiv[], double b[]) - { - return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_dgetrs); - } - - DLLEXPORT lapack_int c_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float b[]) - { - return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_cgetrs); - } - - DLLEXPORT lapack_int z_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double b[]) - { - return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_zgetrs); - } - - DLLEXPORT lapack_int s_lu_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) - { - return lu_solve(n, nrhs, a, b, LAPACK_sgetrf, LAPACK_sgetrs); - } - - DLLEXPORT lapack_int d_lu_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) - { - return lu_solve(n, nrhs, a, b, LAPACK_dgetrf, LAPACK_dgetrs); - } - - DLLEXPORT lapack_int c_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) - { - return lu_solve(n, nrhs, a, b, LAPACK_cgetrf, LAPACK_cgetrs); - } - - DLLEXPORT lapack_int z_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) - { - return lu_solve(n, nrhs, a, b, LAPACK_zgetrf, LAPACK_zgetrs); - } - - DLLEXPORT lapack_int s_cholesky_factor(lapack_int n, float a[]) - { - return cholesky_factor(n, a, LAPACK_spotrf); - } - - DLLEXPORT lapack_int d_cholesky_factor(lapack_int n, double* a) - { - return cholesky_factor(n, a, LAPACK_dpotrf); - } - - DLLEXPORT lapack_int c_cholesky_factor(lapack_int n, lapack_complex_float a[]) - { - return cholesky_factor(n, a, LAPACK_cpotrf); - } - - DLLEXPORT lapack_int z_cholesky_factor(lapack_int n, lapack_complex_double a[]) - { - return cholesky_factor(n, a, LAPACK_zpotrf); - } - - DLLEXPORT lapack_int s_cholesky_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) - { - return cholesky_solve(n, nrhs, a, b, LAPACK_spotrf, LAPACK_spotrs); - } - - DLLEXPORT lapack_int d_cholesky_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) - { - return cholesky_solve(n, nrhs, a, b, LAPACK_dpotrf, LAPACK_dpotrs); - } - - DLLEXPORT lapack_int c_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) - { - return cholesky_solve(n, nrhs, a, b, LAPACK_cpotrf, LAPACK_cpotrs); - } - - DLLEXPORT lapack_int z_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) - { - return cholesky_solve(n, nrhs, a, b, LAPACK_zpotrf, LAPACK_zpotrs); - } - - DLLEXPORT lapack_int s_cholesky_solve_factored(lapack_int n, lapack_int nrhs, float a[], float b[]) - { - return cholesky_solve_factored(n, nrhs, a, b, LAPACK_spotrs); - } - - DLLEXPORT lapack_int d_cholesky_solve_factored(lapack_int n, lapack_int nrhs, double a[], double b[]) - { - return cholesky_solve_factored(n, nrhs, a, b, LAPACK_dpotrs); - } - - DLLEXPORT lapack_int c_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) - { - return cholesky_solve_factored(n, nrhs, a, b, LAPACK_cpotrs); - } - - DLLEXPORT lapack_int z_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) - { - return cholesky_solve_factored(n, nrhs, a, b, LAPACK_zpotrs); - } - - DLLEXPORT lapack_int s_qr_factor(lapack_int m, lapack_int n, float r[], float tau[], float q[], float work[], lapack_int len) - { - return qr_factor(m, n, r, tau, q, work, len, LAPACK_sgeqrf, LAPACK_sorgqr); - } - - DLLEXPORT lapack_int s_qr_thin_factor(lapack_int m, lapack_int n, float q[], float tau[], float r[], float work[], lapack_int len) - { - return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_sgeqrf, LAPACK_sorgqr); - } - - DLLEXPORT lapack_int d_qr_factor(lapack_int m, lapack_int n, double r[], double tau[], double q[], double work[], lapack_int len) - { - return qr_factor(m, n, r, tau, q, work, len, LAPACK_dgeqrf, LAPACK_dorgqr); - } - - DLLEXPORT lapack_int d_qr_thin_factor(lapack_int m, lapack_int n, double q[], double tau[], double r[], double work[], lapack_int len) - { - return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_dgeqrf, LAPACK_dorgqr); - } - - DLLEXPORT lapack_int c_qr_factor(lapack_int m, lapack_int n, lapack_complex_float r[], lapack_complex_float tau[], lapack_complex_float q[], lapack_complex_float work[], lapack_int len) - { - return qr_factor(m, n, r, tau, q, work, len, LAPACK_cgeqrf, LAPACK_cungqr); - } - - DLLEXPORT lapack_int c_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_float q[], lapack_complex_float tau[], lapack_complex_float r[], lapack_complex_float work[], lapack_int len) - { - return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_cgeqrf, LAPACK_cungqr); - } - - DLLEXPORT lapack_int z_qr_factor(lapack_int m, lapack_int n, lapack_complex_double r[], lapack_complex_double tau[], lapack_complex_double q[], lapack_complex_double work[], lapack_int len) - { - return qr_factor(m, n, r, tau, q, work, len, LAPACK_zgeqrf, LAPACK_zungqr); - } - - DLLEXPORT lapack_int z_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_double q[], lapack_complex_double tau[], lapack_complex_double r[], lapack_complex_double work[], lapack_int len) - { - return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_zgeqrf, LAPACK_zungqr); - } - - DLLEXPORT lapack_int s_qr_solve(lapack_int m, lapack_int n, lapack_int bn, float a[], float b[], float x[], float work[], lapack_int len) - { - return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_sgels); - } - - DLLEXPORT lapack_int d_qr_solve(lapack_int m, lapack_int n, lapack_int bn, double a[], double b[], double x[], double work[], lapack_int len) - { - return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_dgels); - } - - DLLEXPORT lapack_int c_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float a[], lapack_complex_float b[], lapack_complex_float x[], lapack_complex_float work[], lapack_int len) - { - return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_cgels); - } - - DLLEXPORT lapack_int z_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double a[], lapack_complex_double b[], lapack_complex_double x[], lapack_complex_double work[], lapack_int len) - { - return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_zgels); - } - - DLLEXPORT lapack_int s_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, float r[], float b[], float tau[], float x[], float work[], lapack_int len) - { - return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_sormqr, cblas_strsm); - } - - DLLEXPORT lapack_int d_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, double r[], double b[], double tau[], double x[], double work[], lapack_int len) - { - return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_dormqr, cblas_dtrsm); - } - - DLLEXPORT lapack_int c_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float r[], lapack_complex_float b[], lapack_complex_float tau[], lapack_complex_float x[], lapack_complex_float work[], lapack_int len) - { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_cunmqr, cblas_ctrsm); - } - - DLLEXPORT lapack_int z_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double r[], lapack_complex_double b[], lapack_complex_double tau[], lapack_complex_double x[], lapack_complex_double work[], lapack_int len) - { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_zunmqr, cblas_ztrsm); - } - - DLLEXPORT lapack_int s_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, float a[], float s[], float u[], float v[], float work[], lapack_int len) - { - return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_sgesvd); - } - - DLLEXPORT lapack_int d_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, double a[], double s[], double u[], double v[], double work[], lapack_int len) - { - return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_dgesvd); - } - - DLLEXPORT lapack_int c_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_float a[], lapack_complex_float s[], lapack_complex_float u[], lapack_complex_float v[], lapack_complex_float work[], lapack_int len) - { - return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_cgesvd); - } - - DLLEXPORT lapack_int z_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_double a[], lapack_complex_double s[], lapack_complex_double u[], lapack_complex_double v[], lapack_complex_double work[], lapack_int len) - { - return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_zgesvd); - } - - DLLEXPORT lapack_int s_eigen(bool isSymmetric, lapack_int n, float a[], float vectors[], lapack_complex_double values[], float d[]) - { - if (isSymmetric) - { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); - } - else - { - return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc); - } - } - - DLLEXPORT lapack_int d_eigen(bool isSymmetric, lapack_int n, double a[], double vectors[], lapack_complex_double values[], double d[]) - { - if (isSymmetric) - { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); - } - else - { - return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); - } - } - - DLLEXPORT lapack_int c_eigen(bool isSymmetric, lapack_int n, lapack_complex_float a[], lapack_complex_float vectors[], lapack_complex_double values[], lapack_complex_float d[]) - { - if (isSymmetric) - { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); - } - else - { - return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); - } - } - - DLLEXPORT lapack_int z_eigen(bool isSymmetric, lapack_int n, lapack_complex_double a[], lapack_complex_double vectors[], lapack_complex_double values[], lapack_complex_double d[]) - { - if (isSymmetric) - { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); - } - else - { - return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); - } - } -} diff --git a/src/NativeProviders/OpenBLAS/lapack.h b/src/NativeProviders/OpenBLAS/lapack.h new file mode 100644 index 00000000..30d40bbe --- /dev/null +++ b/src/NativeProviders/OpenBLAS/lapack.h @@ -0,0 +1,9 @@ +#pragma once + +#define LAPACK_COMPLEX_CUSTOM +#include +#define lapack_complex_float std::complex +#define lapack_complex_double std::complex + +#include "cblas.h" +#include "lapacke.h" diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj index 504a3317..002e17f4 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj @@ -289,10 +289,10 @@ + + - - @@ -300,6 +300,9 @@ + + + diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters index c543369a..ce59f97e 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters @@ -15,12 +15,6 @@ - - Source Files - - - Source Files - Source Files @@ -33,6 +27,12 @@ Source Files + + Source Files + + + Source Files + @@ -43,5 +43,14 @@ Header Files + + Header Files + + + Header Files + + + Header Files + \ No newline at end of file diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj index eb3db333..ad28493e 100644 --- a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj @@ -203,16 +203,18 @@ + + - - - + + + diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters index 3a4caaa1..4f847d3b 100644 --- a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters @@ -15,16 +15,16 @@ - + Source Files - + Source Files - + Source Files - + Source Files @@ -34,10 +34,16 @@ - + Header Files - + + Header Files + + + Header Files + + Header Files diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index 4fd3fee5..825f1463 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -88,7 +88,9 @@ namespace MathNet.Numerics break; default: +#if NATIVE if (!TryUseNative()) +#endif { UseManaged(); } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs index a683b7bf..9b387546 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs @@ -76,8 +76,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); } - var work = new double[rows]; - return SafeNativeMethods.z_matrix_norm((byte)norm, rows, columns, matrix, work); + return SafeNativeMethods.z_matrix_norm((byte)norm, rows, columns, matrix); } /// @@ -229,7 +228,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m * n) + if (c.Length != m*n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -264,7 +263,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order * order) + if (data.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -274,7 +273,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - SafeNativeMethods.z_lu_factor(order, data, ipiv); + var info = SafeNativeMethods.z_lu_factor(order, data, ipiv); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -291,82 +295,27 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - var work = new Complex[order]; - SafeNativeMethods.z_lu_inverse(order, a, work, work.Length); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// This is equivalent to the GETRI LAPACK routine. - [SecuritySafeCritical] - public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - var work = new Complex[order]; - SafeNativeMethods.z_lu_inverse_factored(order, a, ipiv, work, order); - } - - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [SecuritySafeCritical] - public override void LUInverse(Complex[] a, int order, Complex[] work) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } + var info = SafeNativeMethods.z_lu_inverse(order, a); - if (a.Length != order * order) + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new MemoryAllocationException(); } - if (work == null) + if (info < 0) { - throw new ArgumentNullException("work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (work.Length < order) + if (info > 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new SingularUMatrixException(info); } - - SafeNativeMethods.z_lu_inverse(order, a, work, work.Length); } /// @@ -375,12 +324,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// The LU factored N by N matrix. Contains the inverse On exit. /// The order of the square matrix . /// The pivot indices of . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) + public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) { if (a == null) { @@ -392,7 +338,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -402,17 +348,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } + var info = SafeNativeMethods.z_lu_inverse_factored(order, a, ipiv); - if (work.Length < order) + if (info < 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - SafeNativeMethods.z_lu_inverse_factored(order, a, ipiv, work, order); + if (info > 0) + { + throw new SingularUMatrixException(info); + } } /// @@ -431,12 +377,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -446,7 +392,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.z_lu_solve(order, columnsOfB, a, b); + var info = SafeNativeMethods.z_lu_solve(order, columnsOfB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -471,7 +427,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -481,7 +437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -491,7 +447,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.z_lu_solve_factored(order, columnsOfB, a, ipiv, b); + var info = SafeNativeMethods.z_lu_solve_factored(order, columnsOfB, a, ipiv, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -514,13 +480,18 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } var info = SafeNativeMethods.z_cholesky_factor(order, a); + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + if (info > 0) { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); @@ -549,7 +520,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -559,7 +530,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.z_cholesky_solve(orderA, columnsB, a, b); + var info = SafeNativeMethods.z_cholesky_solve(orderA, columnsB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -583,7 +564,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -593,51 +574,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.z_cholesky_solve_factored(orderA, columnsB, a, b); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } + var info = SafeNativeMethods.z_cholesky_solve_factored(orderA, columnsB, a, b); - if (r.Length != rowsR * columnsR) + if (info < 0) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + throw new InvalidParameterException(Math.Abs(info)); } - - if (tau.Length < Math.Min(rowsR, columnsR)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } - - if (q.Length != rowsR * rowsR) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); - } - - var work = new Complex[columnsR * Control.BlockSize]; - SafeNativeMethods.z_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -647,16 +589,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// it is overwritten with the R matrix of the QR factorization. /// The number of rows in the A matrix. /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau, Complex[] work) + public override void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) { if (r == null) { @@ -668,12 +607,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -683,18 +617,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR * rowsR) + if (q.Length != rowsR*rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR * Control.BlockSize) + var info = SafeNativeMethods.z_qr_factor(rowsR, columnsR, r, tau, q); + + if (info < 0) { - work[0] = columnsR * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.z_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -737,66 +670,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); } - var work = new Complex[columnsA * Control.BlockSize]; - SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); - } - - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau, Complex[] work) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (work == null) - { - throw new ArgumentNullException("q"); - } - - if (q.Length != rowsA * columnsA) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); - } - - if (tau.Length < Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } - - if (r.Length != columnsA * columnsA) - { - throw new ArgumentException( - string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); - } + var info = SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r); - if (work.Length < columnsA * Control.BlockSize) + if (info < 0) { - work[0] = columnsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); } /// @@ -812,27 +691,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex[columns * Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B 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, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) { if (a == null) { @@ -849,11 +707,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); @@ -874,56 +727,40 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.RowsLessThanColumns); } - if (work.Length < 1) + var info = SafeNativeMethods.z_qr_solve(rows, columns, columnsB, a, b, x); + + if (info == (int)NativeError.MemoryAllocation) { - work[0] = rows * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new MemoryAllocationException(); } - SafeNativeMethods.z_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex[columnsR * Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); + if (info > 0) + { + throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient, "a"); + } } /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The Q matrix obtained by calling . /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. + /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - /// The work array - only used in the native provider. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) + public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -945,11 +782,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - int rowsQ, columnsQ, rowsR, columnsR; if (method == QRMethod.Full) { @@ -982,15 +814,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } - if (work.Length < 1) - { - work[0] = rowsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - if (method == QRMethod.Full) { - SafeNativeMethods.z_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + var info = SafeNativeMethods.z_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } else { @@ -1000,61 +836,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas } } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (s == null) - { - throw new ArgumentNullException("s"); - } - - if (u == null) - { - throw new ArgumentNullException("u"); - } - - if (vt == null) - { - throw new ArgumentNullException("vt"); - } - - if (u.Length != rowsA * rowsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); - } - - if (vt.Length != columnsA * columnsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); - } - - if (s.Length != Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); - } - - var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -1081,24 +862,23 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex[Math.Min(rowsA, columnsA)]; - var u = new Complex[rowsA * rowsA]; - var vt = new Complex[columnsA * columnsA]; + var u = new Complex[rowsA*rowsA]; + var vt = new Complex[columnsA*columnsA]; var clone = new Complex[a.Length]; a.Copy(clone); - SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } @@ -1114,12 +894,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// singular vectors. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// The work array. For real matrices, the work array should be at least - /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). - /// On exit, work[0] contains the optimal work size value. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) { if (a == null) { @@ -1141,17 +918,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("vt"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (u.Length != rowsA * rowsA) + if (u.Length != rowsA*rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA * columnsA) + if (vt.Length != columnsA*columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1161,18 +933,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - if (work.Length == 0) + var info = SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt); + + if (info == (int) NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + throw new MemoryAllocationException(); } - if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (info < 0) { - work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + if (info > 0) { throw new NonConvergenceException(); } @@ -1229,7 +1002,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); } - if (SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + var info = SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + + if (info > 0) { throw new NonConvergenceException(); } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs index 9f118068..44692f36 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs @@ -76,8 +76,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); } - var work = new float[rows]; - return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix, work); + return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix); } /// @@ -229,7 +228,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m * n) + if (c.Length != m*n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -264,7 +263,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order * order) + if (data.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -274,7 +273,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - SafeNativeMethods.c_lu_factor(order, data, ipiv); + var info = SafeNativeMethods.c_lu_factor(order, data, ipiv); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -291,82 +295,27 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - var work = new Complex32[order]; - SafeNativeMethods.c_lu_inverse(order, a, work, work.Length); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// This is equivalent to the GETRI LAPACK routine. - [SecuritySafeCritical] - public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - var work = new Complex32[order]; - SafeNativeMethods.c_lu_inverse_factored(order, a, ipiv, work, order); - } - - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [SecuritySafeCritical] - public override void LUInverse(Complex32[] a, int order, Complex32[] work) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } + var info = SafeNativeMethods.c_lu_inverse(order, a); - if (a.Length != order * order) + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new MemoryAllocationException(); } - if (work == null) + if (info < 0) { - throw new ArgumentNullException("work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (work.Length < order) + if (info > 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new SingularUMatrixException(info); } - - SafeNativeMethods.c_lu_inverse(order, a, work, work.Length); } /// @@ -375,12 +324,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// The LU factored N by N matrix. Contains the inverse On exit. /// The order of the square matrix . /// The pivot indices of . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) + public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { if (a == null) { @@ -392,7 +338,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -402,17 +348,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } + var info = SafeNativeMethods.c_lu_inverse_factored(order, a, ipiv); - if (work.Length < order) + if (info < 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - SafeNativeMethods.c_lu_inverse_factored(order, a, ipiv, work, order); + if (info > 0) + { + throw new SingularUMatrixException(info); + } } /// @@ -431,12 +377,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -446,7 +392,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.c_lu_solve(order, columnsOfB, a, b); + var info = SafeNativeMethods.c_lu_solve(order, columnsOfB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -471,7 +427,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -481,7 +437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -491,7 +447,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.c_lu_solve_factored(order, columnsOfB, a, ipiv, b); + var info = SafeNativeMethods.c_lu_solve_factored(order, columnsOfB, a, ipiv, b); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -514,13 +475,18 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } var info = SafeNativeMethods.c_cholesky_factor(order, a); + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + if (info > 0) { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); @@ -549,7 +515,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -559,7 +525,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.c_cholesky_solve(orderA, columnsB, a, b); + var info = SafeNativeMethods.c_cholesky_solve(orderA, columnsB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -583,7 +559,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -593,51 +569,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.c_cholesky_solve_factored(orderA, columnsB, a, b); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR * columnsR) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); - } + var info = SafeNativeMethods.c_cholesky_solve_factored(orderA, columnsB, a, b); - if (tau.Length < Math.Min(rowsR, columnsR)) + if (info < 0) { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + throw new InvalidParameterException(Math.Abs(info)); } - - if (q.Length != rowsR * rowsR) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); - } - - var work = new Complex32[columnsR * Control.BlockSize]; - SafeNativeMethods.c_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -647,16 +584,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// it is overwritten with the R matrix of the QR factorization. /// The number of rows in the A matrix. /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau, Complex32[] work) + public override void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) { if (r == null) { @@ -668,12 +602,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -683,18 +612,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR * rowsR) + if (q.Length != rowsR*rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR * Control.BlockSize) + var info = SafeNativeMethods.c_qr_factor(rowsR, columnsR, r, tau, q); + + if (info < 0) { - work[0] = columnsR * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.c_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -737,66 +665,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); } - var work = new Complex32[columnsA * Control.BlockSize]; - SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); - } + var info = SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r); - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau, Complex32[] work) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (info < 0) { - throw new ArgumentNullException("q"); + throw new InvalidParameterException(Math.Abs(info)); } - - if (work == null) - { - throw new ArgumentNullException("q"); - } - - if (q.Length != rowsA * columnsA) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); - } - - if (tau.Length < Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } - - if (r.Length != columnsA * columnsA) - { - throw new ArgumentException( - string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); - } - - if (work.Length < columnsA * Control.BlockSize) - { - work[0] = columnsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - - SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); } /// @@ -812,27 +686,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex32[columns * Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B 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, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) { if (a == null) { @@ -849,11 +702,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); @@ -874,56 +722,40 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.RowsLessThanColumns); } - if (work.Length < 1) + var info = SafeNativeMethods.c_qr_solve(rows, columns, columnsB, a, b, x); + + if (info == (int)NativeError.MemoryAllocation) { - work[0] = rows * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new MemoryAllocationException(); } - SafeNativeMethods.c_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex32[columnsR * Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); + if (info > 0) + { + throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient, "a"); + } } /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The Q matrix obtained by calling . /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. + /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - /// The work array - only used in the native provider. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) + public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -945,11 +777,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - int rowsQ, columnsQ, rowsR, columnsR; if (method == QRMethod.Full) { @@ -982,15 +809,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } - if (work.Length < 1) - { - work[0] = rowsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - if (method == QRMethod.Full) { - SafeNativeMethods.c_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + var info = SafeNativeMethods.c_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } else { @@ -1000,61 +831,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas } } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (s == null) - { - throw new ArgumentNullException("s"); - } - - if (u == null) - { - throw new ArgumentNullException("u"); - } - - if (vt == null) - { - throw new ArgumentNullException("vt"); - } - - if (u.Length != rowsA * rowsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); - } - - if (vt.Length != columnsA * columnsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); - } - - if (s.Length != Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); - } - - var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -1081,24 +857,23 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex32[Math.Min(rowsA, columnsA)]; - var u = new Complex32[rowsA * rowsA]; - var vt = new Complex32[columnsA * columnsA]; + var u = new Complex32[rowsA*rowsA]; + var vt = new Complex32[columnsA*columnsA]; var clone = new Complex32[a.Length]; a.Copy(clone); - SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } @@ -1114,12 +889,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// singular vectors. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// The work array. For real matrices, the work array should be at least - /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). - /// On exit, work[0] contains the optimal work size value. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) { if (a == null) { @@ -1141,11 +913,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("vt"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); @@ -1161,18 +928,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - if (work.Length == 0) + var info = SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt); + + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + throw new MemoryAllocationException(); } - if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (info < 0) { - work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + if (info > 0) { throw new NonConvergenceException(); } @@ -1226,10 +994,22 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (matrixD.Length != order * order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); } - if (SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + var info = SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + + if (info > 0) { throw new NonConvergenceException(); } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs index 0143e0be..637acb01 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs @@ -76,8 +76,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); } - var work = new double[rows]; - return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix, work); + return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix); } /// @@ -229,7 +228,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m * n) + if (c.Length != m*n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -264,7 +263,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order * order) + if (data.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -274,7 +273,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - SafeNativeMethods.d_lu_factor(order, data, ipiv); + var info = SafeNativeMethods.d_lu_factor(order, data, ipiv); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -291,82 +295,27 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - var work = new double[order]; - SafeNativeMethods.d_lu_inverse(order, a, work, work.Length); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// This is equivalent to the GETRI LAPACK routine. - [SecuritySafeCritical] - public override void LUInverseFactored(double[] a, int order, int[] ipiv) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - var work = new double[order]; - SafeNativeMethods.d_lu_inverse_factored(order, a, ipiv, work, order); - } - - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [SecuritySafeCritical] - public override void LUInverse(double[] a, int order, double[] work) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } + var info = SafeNativeMethods.d_lu_inverse(order, a); - if (a.Length != order * order) + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new MemoryAllocationException(); } - if (work == null) + if (info < 0) { - throw new ArgumentNullException("work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (work.Length < order) + if (info > 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new SingularUMatrixException(info); } - - SafeNativeMethods.d_lu_inverse(order, a, work, work.Length); } /// @@ -375,12 +324,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// The LU factored N by N matrix. Contains the inverse On exit. /// The order of the square matrix . /// The pivot indices of . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) + public override void LUInverseFactored(double[] a, int order, int[] ipiv) { if (a == null) { @@ -392,7 +338,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -402,17 +348,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } + var info = SafeNativeMethods.d_lu_inverse_factored(order, a, ipiv); - if (work.Length < order) + if (info < 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - SafeNativeMethods.d_lu_inverse_factored(order, a, ipiv, work, order); + if (info > 0) + { + throw new SingularUMatrixException(info); + } } /// @@ -431,12 +377,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -446,7 +392,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.d_lu_solve(order, columnsOfB, a, b); + var info = SafeNativeMethods.d_lu_solve(order, columnsOfB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -471,7 +427,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -481,7 +437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -491,7 +447,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.d_lu_solve_factored(order, columnsOfB, a, ipiv, b); + var info = SafeNativeMethods.d_lu_solve_factored(order, columnsOfB, a, ipiv, b); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -514,13 +475,23 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } var info = SafeNativeMethods.d_cholesky_factor(order, a); + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + if (info > 0) { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); @@ -549,7 +520,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -559,7 +530,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.d_cholesky_solve(orderA, columnsB, a, b); + var info = SafeNativeMethods.d_cholesky_solve(orderA, columnsB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -583,7 +564,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -593,51 +574,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.d_cholesky_solve_factored(orderA, columnsB, a, b); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } + var info = SafeNativeMethods.d_cholesky_solve_factored(orderA, columnsB, a, b); - if (r.Length != rowsR * columnsR) + if (info < 0) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + throw new InvalidParameterException(Math.Abs(info)); } - - if (tau.Length < Math.Min(rowsR, columnsR)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } - - if (q.Length != rowsR * rowsR) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); - } - - var work = new double[columnsR * Control.BlockSize]; - SafeNativeMethods.d_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -647,16 +589,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// it is overwritten with the R matrix of the QR factorization. /// The number of rows in the A matrix. /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau, double[] work) + public override void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) { if (r == null) { @@ -668,12 +607,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -683,18 +617,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR * rowsR) + if (q.Length != rowsR*rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR * Control.BlockSize) + var info = SafeNativeMethods.d_qr_factor(rowsR, columnsR, r, tau, q); + + if (info < 0) { - work[0] = columnsR * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.d_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -722,7 +655,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (q.Length != rowsA * columnsA) + if (q.Length != rowsA*columnsA) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } @@ -732,71 +665,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (r.Length != columnsA * columnsA) + if (r.Length != columnsA*columnsA) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); } - var work = new double[columnsA * Control.BlockSize]; - SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); - } - - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau, double[] work) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (work == null) - { - throw new ArgumentNullException("q"); - } - - if (q.Length != rowsA * columnsA) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); - } - - if (tau.Length < Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } - - if (r.Length != columnsA * columnsA) - { - throw new ArgumentException( - string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); - } + var info = SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r); - if (work.Length < columnsA * Control.BlockSize) + if (info < 0) { - work[0] = columnsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); } /// @@ -812,27 +691,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) - { - var work = new double[columns * Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B 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, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) { if (a == null) { @@ -849,11 +707,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); @@ -874,56 +727,40 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.RowsLessThanColumns); } - if (work.Length < 1) + var info = SafeNativeMethods.d_qr_solve(rows, columns, columnsB, a, b, x); + + if (info == (int)NativeError.MemoryAllocation) { - work[0] = rows * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new MemoryAllocationException(); } - SafeNativeMethods.d_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) - { - var work = new double[columnsR * Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); + if (info > 0) + { + throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient, "a"); + } } /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The Q matrix obtained by calling . /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. + /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - /// The work array - only used in the native provider. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) + public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -945,11 +782,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - int rowsQ, columnsQ, rowsR, columnsR; if (method == QRMethod.Full) { @@ -982,15 +814,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } - if (work.Length < 1) - { - work[0] = rowsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - if (method == QRMethod.Full) { - SafeNativeMethods.d_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + var info = SafeNativeMethods.d_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } else { @@ -1000,61 +836,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas } } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (s == null) - { - throw new ArgumentNullException("s"); - } - - if (u == null) - { - throw new ArgumentNullException("u"); - } - - if (vt == null) - { - throw new ArgumentNullException("vt"); - } - - if (u.Length != rowsA * rowsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); - } - - if (vt.Length != columnsA * columnsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); - } - - if (s.Length != Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); - } - - var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -1081,24 +862,23 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; var s = new double[Math.Min(rowsA, columnsA)]; - var u = new double[rowsA * rowsA]; - var vt = new double[columnsA * columnsA]; + var u = new double[rowsA*rowsA]; + var vt = new double[columnsA*columnsA]; var clone = new double[a.Length]; a.Copy(clone); - SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } @@ -1114,12 +894,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// singular vectors. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// The work array. For real matrices, the work array should be at least - /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). - /// On exit, work[0] contains the optimal work size value. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work) + public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) { if (a == null) { @@ -1141,11 +918,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("vt"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); @@ -1161,18 +933,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - if (work.Length == 0) + var info = SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt); + + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + throw new MemoryAllocationException(); } - if (work.Length < Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))) + if (info < 0) { - work[0] = Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA)); - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + if (info > 0) { throw new NonConvergenceException(); } @@ -1194,9 +967,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrix"); } - if (matrix.Length != order * order) + if (matrix.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrix"); } if (matrixEv == null) @@ -1204,9 +977,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrixEv"); } - if (matrixEv.Length != order * order) + if (matrixEv.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixEv"); } if (vectorEv == null) @@ -1224,12 +997,24 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrixD"); } - if (matrixD.Length != order * order) + if (matrixD.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixD"); } - if (SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + var info = SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + + if (info > 0) { throw new NonConvergenceException(); } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs index 61c934c1..c228e039 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs @@ -76,8 +76,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); } - var work = new float[rows]; - return SafeNativeMethods.s_matrix_norm((byte)norm, rows, columns, matrix, work); + return SafeNativeMethods.s_matrix_norm((byte)norm, rows, columns, matrix); } /// @@ -229,7 +228,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m * n) + if (c.Length != m*n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -264,7 +263,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order * order) + if (data.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -274,7 +273,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - SafeNativeMethods.s_lu_factor(order, data, ipiv); + var info = SafeNativeMethods.s_lu_factor(order, data, ipiv); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -291,82 +295,27 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - var work = new float[order]; - SafeNativeMethods.s_lu_inverse(order, a, work, work.Length); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// This is equivalent to the GETRI LAPACK routine. - [SecuritySafeCritical] - public override void LUInverseFactored(float[] a, int order, int[] ipiv) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } + var info = SafeNativeMethods.s_lu_inverse(order, a); - if (a.Length != order * order) + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new MemoryAllocationException(); } - if (ipiv.Length != order) + if (info < 0) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + throw new InvalidParameterException(Math.Abs(info)); } - var work = new float[order]; - SafeNativeMethods.s_lu_inverse_factored(order, a, ipiv, work, order); - } - - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [SecuritySafeCritical] - public override void LUInverse(float[] a, int order, float[] work) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (work.Length < order) + if (info > 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new SingularUMatrixException(info); } - - SafeNativeMethods.s_lu_inverse(order, a, work, work.Length); } /// @@ -375,12 +324,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// The LU factored N by N matrix. Contains the inverse On exit. /// The order of the square matrix . /// The pivot indices of . - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) + public override void LUInverseFactored(float[] a, int order, int[] ipiv) { if (a == null) { @@ -392,7 +338,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -402,17 +348,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } + var info = SafeNativeMethods.s_lu_inverse_factored(order, a, ipiv); - if (work.Length < order) + if (info < 0) { - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - SafeNativeMethods.s_lu_inverse_factored(order, a, ipiv, work, order); + if (info > 0) + { + throw new SingularUMatrixException(info); + } } /// @@ -431,12 +377,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("a"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -446,7 +392,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.s_lu_solve(order, columnsOfB, a, b); + var info = SafeNativeMethods.s_lu_solve(order, columnsOfB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -471,7 +427,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -481,7 +437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB * order) + if (b.Length != columnsOfB*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -491,7 +447,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.s_lu_solve_factored(order, columnsOfB, a, ipiv, b); + var info = SafeNativeMethods.s_lu_solve_factored(order, columnsOfB, a, ipiv, b); + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -514,13 +475,18 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order * order) + if (a.Length != order*order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } var info = SafeNativeMethods.s_cholesky_factor(order, a); + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + if (info > 0) { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); @@ -549,7 +515,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -559,7 +525,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.s_cholesky_solve(orderA, columnsB, a, b); + var info = SafeNativeMethods.s_cholesky_solve(orderA, columnsB, a, b); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } /// @@ -583,7 +559,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA * columnsB) + if (b.Length != orderA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -593,51 +569,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - SafeNativeMethods.s_cholesky_solve_factored(orderA, columnsB, a, b); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR * columnsR) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); - } - - if (tau.Length < Math.Min(rowsR, columnsR)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } + var info = SafeNativeMethods.s_cholesky_solve_factored(orderA, columnsB, a, b); - if (q.Length != rowsR * rowsR) + if (info < 0) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); + throw new InvalidParameterException(Math.Abs(info)); } - - var work = new float[columnsR * Control.BlockSize]; - SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -647,16 +584,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// it is overwritten with the R matrix of the QR factorization. /// The number of rows in the A matrix. /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau, float[] work) + public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) { if (r == null) { @@ -668,12 +602,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -683,18 +612,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR * rowsR) + if (q.Length != rowsR*rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR * Control.BlockSize) + var info = SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q); + + if (info < 0) { - work[0] = columnsR * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - - SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -737,66 +665,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); } - var work = new float[columnsA * Control.BlockSize]; - SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); - } - - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// The work array. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - [SecuritySafeCritical] - public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau, float[] work) - { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (work == null) - { - throw new ArgumentNullException("q"); - } - - if (q.Length != rowsA * columnsA) - { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); - } - - if (tau.Length < Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); - } + var info = SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r); - if (r.Length != columnsA * columnsA) + if (info < 0) { - throw new ArgumentException( - string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + throw new InvalidParameterException(Math.Abs(info)); } - - if (work.Length < columnsA * Control.BlockSize) - { - work[0] = columnsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - - SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); } /// @@ -812,27 +686,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) - { - var work = new float[columns * Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B 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, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) { if (a == null) { @@ -849,11 +702,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); @@ -874,56 +722,40 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.RowsLessThanColumns); } - if (work.Length < 1) + var info = SafeNativeMethods.s_qr_solve(rows, columns, columnsB, a, b, x); + + if (info == (int)NativeError.MemoryAllocation) { - work[0] = rows * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new MemoryAllocationException(); } - SafeNativeMethods.s_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - [SecuritySafeCritical] - public override void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) - { - var work = new float[columnsR * Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); + if (info > 0) + { + throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient, "a"); + } } /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The Q matrix obtained by calling . /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. + /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - /// The work array - only used in the native provider. The array must have a length of at least N, - /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal - /// work size value. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) + public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -945,11 +777,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("q"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - int rowsQ, columnsQ, rowsR, columnsR; if (method == QRMethod.Full) { @@ -982,15 +809,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } - if (work.Length < 1) - { - work[0] = rowsA * Control.BlockSize; - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); - } - if (method == QRMethod.Full) { - SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + var info = SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } } else { @@ -1000,61 +831,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas } } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (s == null) - { - throw new ArgumentNullException("s"); - } - - if (u == null) - { - throw new ArgumentNullException("u"); - } - - if (vt == null) - { - throw new ArgumentNullException("vt"); - } - - if (u.Length != rowsA * rowsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); - } - - if (vt.Length != columnsA * columnsA) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); - } - - if (s.Length != Math.Min(rowsA, columnsA)) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); - } - - var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -1081,24 +857,23 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; var s = new float[Math.Min(rowsA, columnsA)]; - var u = new float[rowsA * rowsA]; - var vt = new float[columnsA * columnsA]; + var u = new float[rowsA*rowsA]; + var vt = new float[columnsA*columnsA]; var clone = new float[a.Length]; a.Copy(clone); - SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } @@ -1114,12 +889,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// singular vectors. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// The work array. For real matrices, the work array should be at least - /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). - /// On exit, work[0] contains the optimal work size value. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work) + public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) { if (a == null) { @@ -1141,17 +913,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("vt"); } - if (work == null) - { - throw new ArgumentNullException("work"); - } - - if (u.Length != rowsA * rowsA) + if (u.Length != rowsA*rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA * columnsA) + if (vt.Length != columnsA*columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1161,18 +928,19 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - if (work.Length == 0) + var info = SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt); + + if (info == (int)NativeError.MemoryAllocation) { - throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + throw new MemoryAllocationException(); } - if (work.Length < Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))) + if (info < 0) { - work[0] = Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA)); - throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + throw new InvalidParameterException(Math.Abs(info)); } - if (SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + if (info > 0) { throw new NonConvergenceException(); } @@ -1194,9 +962,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrix"); } - if (matrix.Length != order * order) + if (matrix.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrix"); } if (matrixEv == null) @@ -1204,9 +972,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrixEv"); } - if (matrixEv.Length != order * order) + if (matrixEv.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixEv"); } if (vectorEv == null) @@ -1224,12 +992,24 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas throw new ArgumentNullException("matrixD"); } - if (matrixD.Length != order * order) + if (matrixD.Length != order*order) { - throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixD"); } - if (SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + var info = SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + + if (info == (int)NativeError.MemoryAllocation) + { + throw new MemoryAllocationException(); + } + + if (info < 0) + { + throw new InvalidParameterException(Math.Abs(info)); + } + + if (info > 0) { throw new NonConvergenceException(); } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs index 5e4f5f61..ba7d8876 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -34,6 +34,17 @@ using System; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { + /// + /// Error codes return from the native OpenBLAS provider. + /// + public enum NativeError : int + { + /// + /// Unable to allocate memory. + /// + MemoryAllocation = -999999 + } + internal enum ParallelType : int { Sequential = 0, diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs index fa59a5d7..fc024d90 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -117,16 +117,28 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas #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); + internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a); [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 double d_matrix_norm(byte norm, int rows, int columns, [In] double[] a); [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 float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a); [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 double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_factor(int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern int s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); @@ -141,28 +153,28 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas internal static extern int z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork); + internal static extern int s_lu_inverse(int n, [In, Out] float[] a); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork); + internal static extern int d_lu_inverse(int n, [In, Out] double[] a); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork); + internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork); + internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork); + internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork); + internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork); + internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); + internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern int s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); @@ -188,18 +200,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas [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_factor(int n, [In, Out] float[] a); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); - [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); @@ -225,77 +225,77 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); + internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_thin_factor(int m, int n, [In, Out] float[] q, [In, Out] float[] tau, [In, Out] float[] r, [In, Out] float[] work, int len); + internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); + internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_thin_factor(int m, int n, [In, Out] double[] q, [In, Out] double[] tau, [In, Out] double[] r, [In, Out] double[] work, int len); + internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); + internal static extern int s_qr_thin_factor(int m, int n, [In, Out] float[] q, [In, Out] float[] tau, [In, Out] float[] r); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_thin_factor(int m, int n, [In, Out] Complex32[] q, [In, Out] Complex32[] tau, [In, Out] Complex32[] r, [In, Out] Complex32[] work, int len); + internal static extern int d_qr_thin_factor(int m, int n, [In, Out] double[] q, [In, Out] double[] tau, [In, Out] double[] r); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); + internal static extern int c_qr_thin_factor(int m, int n, [In, Out] Complex32[] q, [In, Out] Complex32[] tau, [In, Out] Complex32[] r); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_thin_factor(int m, int n, [In, Out] Complex[] q, [In, Out] Complex[] tau, [In, Out] Complex[] r, [In, Out] Complex[] work, int len); + internal static extern int z_qr_thin_factor(int m, int n, [In, Out] Complex[] q, [In, Out] Complex[] tau, [In, Out] Complex[] r); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); + internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); + internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); + internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); + internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_svd_factor(bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len); + internal static extern int s_svd_factor([MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_svd_factor(bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len); + internal static extern int d_svd_factor([MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len); + internal static extern int c_svd_factor([MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len); + internal static extern int z_svd_factor([MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_eigen(bool isSymmetric, int n, [In, Out] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d); + internal static extern int s_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_eigen(bool isSymmetric, int n, [In, Out] double[] a, [In, Out] double[] vectors, [In, Out] Complex[] values, [In, Out] double[] d); + internal static extern int d_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] double[] a, [In, Out] double[] vectors, [In, Out] Complex[] values, [In, Out] double[] d); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_eigen(bool isSymmetric, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out] Complex[] values, [In, Out] Complex32[] d); + internal static extern int c_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out] Complex[] values, [In, Out] Complex32[] d); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_eigen(bool isSymmetric, int n, [In, Out] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); - + internal static extern int z_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); + #endregion LAPACK } }