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