From d26d41e848ce8fa3932b214319c42285e85b1cb4 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Thu, 13 Dec 2012 10:23:12 +0200 Subject: [PATCH] removed redundant code in lapack.cpp, using gels for all qr solving (dropping gelsy) --- src/NativeWrappers/MKL/lapack.cpp | 1238 +++++------------ .../MKLWrapperTests/MKLWrapperTests.csproj | 8 +- .../Mkl/MklLinearAlgebraProvider.Complex.cs | 9 +- .../Mkl/MklLinearAlgebraProvider.Complex32.cs | 9 +- .../Mkl/MklLinearAlgebraProvider.double.cs | 9 +- .../Mkl/MklLinearAlgebraProvider.float.cs | 9 +- .../LinearAlgebra/Mkl/SafeNativeMethods.cs | 12 - .../Single/Factorization/QRTests.cs | 2 +- 8 files changed, 391 insertions(+), 905 deletions(-) diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 92f7d827..09bd04ce 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -3,1077 +3,603 @@ #include "wrapper_common.h" #include -extern "C"{ - DLLEXPORT float s_matrix_norm(char norm, MKL_INT m, MKL_INT n, float a[], float work[]) +template +inline void copyBtoX (MKL_INT m, MKL_INT n, MKL_INT bn, T b[], T x[]){ + for (MKL_INT i = 0; i < n; ++i) { - return slange_(&norm, &m, &n, a, &m, work); + for (MKL_INT j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } } +}; - DLLEXPORT double d_matrix_norm(char norm, MKL_INT m, MKL_INT n, double a[], double work[]) - { - return dlange_(&norm, &m, &n, a, &m, work); +inline void shift_ipiv_down(MKL_INT m, MKL_INT ipiv[]){ + for(MKL_INT i = 0; i < m; ++i ){ + ipiv[i] -= 1; } +} - DLLEXPORT float c_matrix_norm(char norm, MKL_INT m, MKL_INT n, MKL_Complex8 a[], float work[]) - { - return clange_(&norm, &m, &n, a, &m, work); +inline void shift_ipiv_up(MKL_INT m, MKL_INT ipiv[]){ + for(MKL_INT i = 0; i < m; ++i ){ + ipiv[i] += 1; } +} - DLLEXPORT double z_matrix_norm(char norm, MKL_INT m, MKL_INT n, MKL_Complex16 a[], double work[]) - { - return zlange_(&norm, &m, &n, a, &m, work); +template +inline MKL_INT lu_factor(MKL_INT m, T a[], MKL_INT ipiv[], + void (*getrf) (const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*)) +{ + MKL_INT info = 0; + getrf(&m, &m, a, &m, ipiv, &info); + shift_ipiv_down(m, ipiv); + return info; +}; + +template +inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork, + void (*getrf) (const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*), + void (*getri) (const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + MKL_INT* ipiv = new MKL_INT[n]; + MKL_INT info = 0; + getrf(&n, &n, a, &n, ipiv, &info); + + if (info != 0){ + delete[] ipiv; + return info; } - DLLEXPORT MKL_INT s_lu_factor(MKL_INT m, float a[], MKL_INT ipiv[]) + getri(&n, a, &n, ipiv, work, &lwork, &info); + delete[] ipiv; + return info; +}; + +template +inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], MKL_INT lwork, + void (*getri) (const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + shift_ipiv_up(n, ipiv); + MKL_INT info = 0; + getri(&n, a, &n, ipiv, work, &lwork, &info); + shift_ipiv_down(n, ipiv); + return info; +} + +template +inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[], T b[], + void (*getrs) (const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + shift_ipiv_up(n, ipiv); + MKL_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 MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], + void (*getrf) (const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*), + void (*getrs) (const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + T* clone = new T[n*n]; + std::memcpy(clone, a, n*n*sizeof(T)); + + MKL_INT* ipiv = new MKL_INT[n]; + MKL_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 MKL_INT cholesky_factor(MKL_INT n, T* a, + void (*potrf) (const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + char uplo = 'L'; + MKL_INT info = 0; + potrf(&uplo, &n, a, &n, &info); + T zero = T(); + for (MKL_INT i = 0; i < n; ++i) { - MKL_INT info = 0; - sgetrf_(&m,&m,a,&m,ipiv,&info); - for(MKL_INT i = 0; i < m; ++i ){ - ipiv[i] -= 1; + MKL_INT index = i * n; + for (MKL_INT j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; } + } + return info; +} + +template +inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], + void (*potrf) (const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*), + void (*potrs) (const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + T* clone = new T[n*n]; + std::memcpy(clone, a, n*n*sizeof(T)); + char uplo = 'L'; + MKL_INT info = 0; + potrf(&uplo, &n, clone, &n, &info); + + if (info != 0){ + delete[] clone; return info; } - DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[]) + potrs(&uplo, &n, &nrhs, clone, &n, b, &n, &info); + delete[] clone; + return info; +} + +template +inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[], + void (*potrs) (const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + char uplo = 'L'; + MKL_INT info = 0; + potrs(&uplo, &n, &nrhs, a, &n, b, &n, &info); + return info; +} + +template +inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], MKL_INT len, + void (*geqrf) (const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*), + void (*orgqr) (const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*)) +{ + MKL_INT info = 0; + geqrf(&m, &n, r, &m, tau, work, &len, &info); + + for (MKL_INT i = 0; i < m; ++i) { - MKL_INT info = 0; - dgetrf_(&m,&m,a,&m,ipiv,&info); - for(MKL_INT i = 0; i < m; ++i ){ - ipiv[i] -= 1; + for (MKL_INT j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } } - return info; } - DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[]) + //compute the q elements explicitly + if (m <= n) + { + orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else { - MKL_INT info = 0; - cgetrf_(&m,&m,a,&m,ipiv,&info); - for(MKL_INT i = 0; i < m; ++i ){ - ipiv[i] -= 1; + orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); + } + + return info; +} + +template +inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T work[], MKL_INT len, + void (*geqrf) (const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*), + void (*orgqr) (const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*)) +{ + MKL_INT info = 0; + geqrf(&m, &n, q, &m, tau, work, &len, &info); + + for (MKL_INT i = 0; i < n; ++i) + { + for (MKL_INT j = 0; j < n; ++j) + { + if( i <= j) { + r[j * n + i] = q[j * m + i]; + } } - return info; } - DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[]) + orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); + + return info; +} + +template +inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T work[], MKL_INT len, + void (*gels) (const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, + const MKL_INT*, T* b, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + T* clone_a = new T[m*n]; + std::memcpy(clone_a, a, m*n*sizeof(T)); + + T* clone_b = new T[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(T)); + + char N = 'N'; + MKL_INT info = 0; + gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); + + for (MKL_INT i = 0; i < n; ++i) { - MKL_INT info = 0; - zgetrf_(&m,&m,a,&m,ipiv,&info); - for(MKL_INT i = 0; i < m; ++i ){ - ipiv[i] -= 1; + for (MKL_INT j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; } - return info; } - DLLEXPORT MKL_INT s_lu_inverse(MKL_INT n, float a[], float work[], MKL_INT lwork) + delete[] clone_a; + delete[] clone_b; + return info; +} + +// combine the next two some how +// the problem is that complex trsm takes void* instead of MKL_COMPLEX +template +inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len, + void (*ormqr) (const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*, + const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info), + void (*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, + const MKL_INT, const MKL_INT, const T, const T*, const MKL_INT, T*, const MKL_INT)) +{ + T* clone_b = new T[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(T)); + + char side ='L'; + char tran = 'T'; + MKL_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); + for (MKL_INT i = 0; i < n; ++i) { - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - sgetrf_(&n,&n,a,&n,ipiv,&info); + for (MKL_INT j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } + } + + delete[] clone_b; + return info; +} + +template +inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len, + void (*unmqr) (const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*, + const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info), + void (*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, + const MKL_INT, const MKL_INT, const void*, const void*, const MKL_INT, void*, const MKL_INT ldb)) +{ + T* clone_b = new T[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(T)); - if (info != 0){ - delete[] ipiv; - return info; + char side ='L'; + char tran = 'C'; + MKL_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, &one, r, m, clone_b, m); + for (MKL_INT i = 0; i < n; ++i) + { + for (MKL_INT j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; } + } - sgetri_(&n,a,&n,ipiv,work,&lwork,&info); - delete[] ipiv; - return info; + delete[] clone_b; + return info; +} + +template +inline MKL_INT svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len, + void (*gesvd) (const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, + T*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*)) +{ + MKL_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 MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len, + void (*gesvd) (const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, + R*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, R*, MKL_INT*)) +{ + MKL_INT info = 0; + MKL_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(MKL_INT index = 0; index < dim_s; ++index){ + T value = {s_local[index], 0.0f}; + s[index] = value; } - DLLEXPORT MKL_INT d_lu_inverse(MKL_INT n, double a[], double work[], MKL_INT lwork) + delete[] rwork; + delete[] s_local; + return info; +} + +extern "C" { + + DLLEXPORT float s_matrix_norm(char norm, MKL_INT m, MKL_INT n, float a[], float work[]) { - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - dgetrf_(&n,&n,a,&n,ipiv,&info); + return slange(&norm, &m, &n, a, &m, work); + } - if (info != 0){ - delete[] ipiv; - return info; - } + DLLEXPORT double d_matrix_norm(char norm, MKL_INT m, MKL_INT n, double a[], double work[]) + { + return dlange(&norm, &m, &n, a, &m, work); + } - dgetri_(&n,a,&n,ipiv,work,&lwork,&info); - delete[] ipiv; - return info; + DLLEXPORT float c_matrix_norm(char norm, MKL_INT m, MKL_INT n, MKL_Complex8 a[], float work[]) + { + return clange(&norm, &m, &n, a, &m, work); } - DLLEXPORT MKL_INT c_lu_inverse(MKL_INT n, MKL_Complex8 a[], MKL_Complex8 work[], MKL_INT lwork) + DLLEXPORT double z_matrix_norm(char norm, MKL_INT m, MKL_INT n, MKL_Complex16 a[], double work[]) { - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - cgetrf_(&n,&n,a,&n,ipiv,&info); + return zlange(&norm, &m, &n, a, &m, work); + } - if (info != 0){ - delete[] ipiv; - return info; - } + DLLEXPORT MKL_INT s_lu_factor(MKL_INT m, float a[], MKL_INT ipiv[]) { + return lu_factor(m, a, ipiv, sgetrf); + } + + DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[]) { + return lu_factor(m, a, ipiv, dgetrf); + } + + DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[]) { + return lu_factor(m, a, ipiv, cgetrf); + } + + DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[]) { + return lu_factor(m, a, ipiv, zgetrf); + } - cgetri_(&n,a,&n,ipiv,work,&lwork,&info); - delete[] ipiv; - return info; + DLLEXPORT MKL_INT s_lu_inverse(MKL_INT n, float a[], float work[], MKL_INT lwork) + { + return lu_inverse(n, a, work, lwork, sgetrf, sgetri); } - DLLEXPORT MKL_INT z_lu_inverse(MKL_INT n, MKL_Complex16 a[], MKL_Complex16 work[], MKL_INT lwork) + DLLEXPORT MKL_INT d_lu_inverse(MKL_INT n, double a[], double work[], MKL_INT lwork) { - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - zgetrf_(&n,&n,a,&n,ipiv,&info); + return lu_inverse(n, a, work, lwork, dgetrf, dgetri); + } - if (info != 0){ - delete[] ipiv; - return info; - } + DLLEXPORT MKL_INT c_lu_inverse(MKL_INT n, MKL_Complex8 a[], MKL_Complex8 work[], MKL_INT lwork) + { + return lu_inverse(n, a, work, lwork, cgetrf, cgetri); + } - zgetri_(&n,a,&n,ipiv,work,&lwork,&info); - delete[] ipiv; - return info; + DLLEXPORT MKL_INT z_lu_inverse(MKL_INT n, MKL_Complex16 a[], MKL_Complex16 work[], MKL_INT lwork) + { + return lu_inverse(n, a, work, lwork, zgetrf, zgetri); } DLLEXPORT MKL_INT s_lu_inverse_factored(MKL_INT n, float a[], MKL_INT ipiv[], float work[], MKL_INT lwork) { - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - MKL_INT info = 0; - sgetri_(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_inverse_factored(n, a, ipiv, work, lwork, sgetri); } DLLEXPORT MKL_INT d_lu_inverse_factored(MKL_INT n, double a[], MKL_INT ipiv[], double work[], MKL_INT lwork) { - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - MKL_INT info = 0; - dgetri_(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_inverse_factored(n, a, ipiv, work, lwork, dgetri); } DLLEXPORT MKL_INT c_lu_inverse_factored(MKL_INT n, MKL_Complex8 a[], MKL_INT ipiv[], MKL_Complex8 work[], MKL_INT lwork) { - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - MKL_INT info = 0; - cgetri_(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_inverse_factored(n, a, ipiv, work, lwork, cgetri); } DLLEXPORT MKL_INT z_lu_inverse_factored(MKL_INT n, MKL_Complex16 a[], MKL_INT ipiv[], MKL_Complex16 work[], MKL_INT lwork) { - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - MKL_INT info = 0; - zgetri_(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_inverse_factored(n, a, ipiv, work, lwork, zgetri); } DLLEXPORT MKL_INT s_lu_solve_factored(MKL_INT n, MKL_INT nrhs, float a[], MKL_INT ipiv[], float b[]) { - MKL_INT info = 0; - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans ='N'; - sgetrs_(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_solve_factored(n, nrhs, a, ipiv, b, sgetrs); } DLLEXPORT MKL_INT d_lu_solve_factored(MKL_INT n, MKL_INT nrhs, double a[], MKL_INT ipiv[], double b[]) { - MKL_INT info = 0; - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans ='N'; - dgetrs_(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_solve_factored(n, nrhs, a, ipiv, b, dgetrs); } DLLEXPORT MKL_INT c_lu_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_INT ipiv[], MKL_Complex8 b[]) { - MKL_INT info = 0; - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans ='N'; - cgetrs_(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_solve_factored(n, nrhs, a, ipiv, b, cgetrs); } DLLEXPORT MKL_INT z_lu_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_INT ipiv[], MKL_Complex16 b[]) { - MKL_INT info = 0; - MKL_INT i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans ='N'; - zgetrs_(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - return info; + return lu_solve_factored(n, nrhs, a, ipiv, b, zgetrs); } DLLEXPORT MKL_INT s_lu_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[]) { - float* clone = new float[n*n]; - std::memcpy(clone, a, n*n*sizeof(float)); - - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - sgetrf_(&n, &n, clone, &n, ipiv, &info); - - if (info != 0){ - delete[] ipiv; - delete[] clone; - return info; - } - - char trans ='N'; - sgetrs_(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); - delete[] ipiv; - delete[] clone; - return info; + return lu_solve(n, nrhs, a, b, sgetrf, sgetrs); } DLLEXPORT MKL_INT d_lu_solve(MKL_INT n, MKL_INT nrhs, double a[], double b[]) { - double* clone = new double[n*n]; - std::memcpy(clone, a, n*n*sizeof(double)); - - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - dgetrf_(&n, &n, clone, &n, ipiv, &info); - - if (info != 0){ - delete[] ipiv; - delete[] clone; - return info; - } - - char trans ='N'; - dgetrs_(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); - delete[] ipiv; - delete[] clone; - return info; + return lu_solve(n, nrhs, a, b, dgetrf, dgetrs); } DLLEXPORT MKL_INT c_lu_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) { - MKL_Complex8* clone = new MKL_Complex8[n*n]; - std::memcpy(clone, a, n*n*sizeof(MKL_Complex8)); - - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - cgetrf_(&n, &n, clone, &n, ipiv, &info); - - if (info != 0){ - delete[] ipiv; - delete[] clone; - return info; - } - - char trans ='N'; - cgetrs_(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); - delete[] ipiv; - delete[] clone; - return info; + return lu_solve(n, nrhs, a, b, cgetrf, cgetrs); } DLLEXPORT MKL_INT z_lu_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) { - MKL_Complex16* clone = new MKL_Complex16[n*n]; - std::memcpy(clone, a, n*n*sizeof(MKL_Complex16)); - - MKL_INT* ipiv = new MKL_INT[n]; - MKL_INT info = 0; - zgetrf_(&n, &n, clone, &n, ipiv, &info); - - if (info != 0){ - delete[] ipiv; - delete[] clone; - return info; - } - - char trans ='N'; - zgetrs_(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); - delete[] ipiv; - delete[] clone; - return info; + return lu_solve(n, nrhs, a, b, zgetrf, zgetrs); } DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[]){ - char uplo = 'L'; - MKL_INT info = 0; - spotrf_(&uplo, &n, a, &n, &info); - for (MKL_INT i = 0; i < n; ++i) - { - MKL_INT index = i * n; - for (MKL_INT j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; + return cholesky_factor(n, a, spotrf); } DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a){ - char uplo = 'L'; - MKL_INT info = 0; - dpotrf_(&uplo, &n, a, &n, &info); - for (MKL_INT i = 0; i < n; ++i) - { - MKL_INT index = i * n; - for (MKL_INT j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; + return cholesky_factor(n, a, dpotrf); } DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[]){ - char uplo = 'L'; - MKL_INT info = 0; - MKL_Complex8 zero = {0.0f, 0.0f}; - cpotrf_(&uplo, &n, a, &n, &info); - for (MKL_INT i = 0; i < n; ++i) - { - MKL_INT index = i * n; - for (MKL_INT j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; + return cholesky_factor(n, a, cpotrf); } DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[]){ - char uplo = 'L'; - MKL_INT info = 0; - MKL_Complex16 zero = {0.0, 0.0}; - zpotrf_(&uplo, &n, a, &n, &info); - for (MKL_INT i = 0; i < n; ++i) - { - MKL_INT index = i * n; - for (MKL_INT j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; + return cholesky_factor(n, a, zpotrf); } DLLEXPORT MKL_INT s_cholesky_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[]) { - float* clone = new float[n*n]; - std::memcpy(clone, a, n*n*sizeof(float)); - char uplo = 'L'; - MKL_INT info = 0; - spotrf_(&uplo, &n, clone, &n, &info); - - if (info != 0){ - delete[] clone; - return info; - } - - spotrs_(&uplo, &n, &nrhs, clone, &n, b, &n, &info); - delete[] clone; - return info; + return cholesky_solve(n, nrhs, a, b, spotrf, spotrs); } DLLEXPORT MKL_INT d_cholesky_solve(MKL_INT n, MKL_INT nrhs, double a[], double b[]) { - double* clone = new double[n*n]; - std::memcpy(clone, a, n*n*sizeof(double)); - char uplo = 'L'; - MKL_INT info = 0; - dpotrf_(&uplo, &n, clone, &n, &info); - - if (info != 0){ - delete[] clone; - return info; - } - - dpotrs_(&uplo, &n, &nrhs, clone, &n, b, &n, &info); - delete[] clone; - return info; + return cholesky_solve(n, nrhs, a, b, dpotrf, dpotrs); } DLLEXPORT MKL_INT c_cholesky_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) { - MKL_Complex8* clone = new MKL_Complex8[n*n]; - std::memcpy(clone, a, n*n*sizeof(MKL_Complex8)); - char uplo = 'L'; - MKL_INT info = 0; - cpotrf_(&uplo, &n, clone, &n, &info); - - if (info != 0){ - delete[] clone; - return info; - } - - cpotrs_(&uplo, &n, &nrhs, clone, &n, b, &n, &info); - delete[] clone; - return info; + return cholesky_solve(n, nrhs, a, b, cpotrf, cpotrs); } DLLEXPORT MKL_INT z_cholesky_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) { - MKL_Complex16* clone = new MKL_Complex16[n*n]; - std::memcpy(clone, a, n*n*sizeof(MKL_Complex16)); - char uplo = 'L'; - MKL_INT info = 0; - zpotrf_(&uplo, &n, clone, &n, &info); - - if (info != 0){ - delete[] clone; - return info; - } - - zpotrs_(&uplo, &n, &nrhs, clone, &n, b, &n, &info); - delete[] clone; - return info; + return cholesky_solve(n, nrhs, a, b, zpotrf, zpotrs); } DLLEXPORT MKL_INT s_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, float a[], float b[]) { - char uplo = 'L'; - MKL_INT info = 0; - spotrs_(&uplo, &n, &nrhs, a, &n, b, &n, &info); - return info; + return cholesky_solve_factored(n, nrhs, a, b, spotrs); } DLLEXPORT MKL_INT d_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, double a[], double b[]) { - char uplo = 'L'; - MKL_INT info = 0; - dpotrs_(&uplo, &n, &nrhs, a, &n, b, &n, &info); - return info; + return cholesky_solve_factored(n, nrhs, a, b, dpotrs); } DLLEXPORT MKL_INT c_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) { - char uplo = 'L'; - MKL_INT info = 0; - cpotrs_(&uplo, &n, &nrhs, a, &n, b, &n, &info); - return info; + return cholesky_solve_factored(n, nrhs, a, b, cpotrs); } DLLEXPORT MKL_INT z_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) { - char uplo = 'L'; - MKL_INT info = 0; - zpotrs_(&uplo, &n, &nrhs, a, &n, b, &n, &info); - return info; + return cholesky_solve_factored(n, nrhs, a, b, zpotrs); } DLLEXPORT MKL_INT s_qr_factor(MKL_INT m, MKL_INT n, float r[], float tau[], float q[], float work[], MKL_INT len) { - MKL_INT info = 0; - sgeqrf_(&m, &n, r, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < m; ++i) - { - for (MKL_INT j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - sorgqr_(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - sorgqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - } - - return info; + return qr_factor(m, n, r, tau, q, work, len, sgeqrf, sorgqr); } DLLEXPORT MKL_INT s_qr_thin_factor(MKL_INT m, MKL_INT n, float q[], float tau[], float r[], float work[], MKL_INT len) { - MKL_INT info = 0; - sgeqrf_(&m, &n, q, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < n; ++j) - { - if( i <= j) { - r[j * n + i] = q[j * m + i]; - } - } - } - - sorgqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - - return info; + return qr_thin_factor(m, n, q, tau, r, work, len, sgeqrf, sorgqr); } - DLLEXPORT MKL_INT d_qr_factor(MKL_INT m, MKL_INT n, double r[], double tau[], double q[], double work[], MKL_INT len) { - MKL_INT info = 0; - dgeqrf_(&m, &n, r, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < m; ++i) - { - for (MKL_INT j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - dorgqr_(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - dorgqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - } - - return info; + return qr_factor(m, n, r, tau, q, work, len, dgeqrf, dorgqr); } DLLEXPORT MKL_INT d_qr_thin_factor(MKL_INT m, MKL_INT n, double q[], double tau[], double r[], double work[], MKL_INT len) { - MKL_INT info = 0; - dgeqrf_(&m, &n, q, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < n; ++j) - { - if( i <= j) { - r[j * n + i] = q[j * m + i]; - } - } - } - - dorgqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - - return info; + return qr_thin_factor(m, n, q, tau, r, work, len, dgeqrf, dorgqr); } DLLEXPORT MKL_INT c_qr_factor(MKL_INT m, MKL_INT n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], MKL_INT len) { - MKL_INT info = 0; - cgeqrf_(&m, &n, r, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < m; ++i) - { - for (MKL_INT j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - cungqr_(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - cungqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - } - - return info; + return qr_factor(m, n, r, tau, q, work, len, cgeqrf, cungqr); } DLLEXPORT MKL_INT c_qr_thin_factor(MKL_INT m, MKL_INT n, MKL_Complex8 q[], MKL_Complex8 tau[], MKL_Complex8 r[], MKL_Complex8 work[], MKL_INT len) { - MKL_INT info = 0; - cgeqrf_(&m, &n, q, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < n; ++j) - { - if( i <= j) { - r[j * n + i] = q[j * m + i]; - } - } - } - - cungqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - - return info; + return qr_thin_factor(m, n, q, tau, r, work, len, cgeqrf, cungqr); } DLLEXPORT MKL_INT z_qr_factor(MKL_INT m, MKL_INT n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], MKL_INT len) { - MKL_INT info = 0; - zgeqrf_(&m, &n, r, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < m; ++i) - { - for (MKL_INT j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - zungqr_(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - zungqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - } - - return info; + return qr_factor(m, n, r, tau, q, work, len, zgeqrf, zungqr); } DLLEXPORT MKL_INT z_qr_thin_factor(MKL_INT m, MKL_INT n, MKL_Complex16 q[], MKL_Complex16 tau[], MKL_Complex16 r[], MKL_Complex16 work[], MKL_INT len) { - MKL_INT info = 0; - zgeqrf_(&m, &n, q, &m, tau, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < n; ++j) - { - if( i <= j) { - r[j * n + i] = q[j * m + i]; - } - } - } - - zungqr_(&m, &n, &n, q, &m, tau, work, &len, &info); - - return info; + return qr_thin_factor(m, n, q, tau, r, work, len, zgeqrf, zungqr); } DLLEXPORT MKL_INT s_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, float a[], float b[], float x[], float work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT* jpvt = new MKL_INT[n]; - MKL_INT rank = 0; - float cond = -1.0; - - float* clone_a = new float[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(float)); - - float* clone_b = new float[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(float)); - - sgelsy_(&m, &n, &bn, clone_a, &m, clone_b, &m, jpvt, &cond, &rank, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] jpvt; - delete[] clone_a; - delete[] clone_b; - return info; + return qr_solve(m, n, bn, a, b, x, work, len, sgels); } DLLEXPORT MKL_INT d_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, double a[], double b[], double x[], double work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT* jpvt = new MKL_INT[n]; - MKL_INT rank = 0; - double cond = -1.0; - - double* clone_a = new double[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(double)); - - double* clone_b = new double[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(double)); - - dgelsy_(&m, &n, &bn, clone_a, &m, clone_b, &m, jpvt, &cond, &rank, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] jpvt; - delete[] clone_a; - delete[] clone_b; - return info; + return qr_solve(m, n, bn, a, b, x, work, len, dgels); } DLLEXPORT MKL_INT c_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex8 a[], MKL_Complex8 b[], MKL_Complex8 x[], MKL_Complex8 work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT* jpvt = new MKL_INT[n]; - float* rwork = new float[2*n]; - MKL_INT rank = 0; - float cond = -1.0; - - MKL_Complex8* clone_a = new MKL_Complex8[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(MKL_Complex8)); - - MKL_Complex8* clone_b = new MKL_Complex8[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex8)); - - cgelsy_(&m, &n, &bn, clone_a, &m, clone_b, &m, jpvt, &cond, &rank, work, &len, rwork, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] jpvt; - delete[] rwork; - delete[] clone_a; - delete[] clone_b; - return info; + return qr_solve(m, n, bn, a, b, x, work, len, cgels); } DLLEXPORT MKL_INT z_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex16 a[], MKL_Complex16 b[], MKL_Complex16 x[], MKL_Complex16 work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT* jpvt = new MKL_INT[n]; - double* rwork = new double[2*n]; - MKL_INT rank = 0; - double cond = -1.0; - - MKL_Complex16* clone_a = new MKL_Complex16[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(MKL_Complex16)); - - MKL_Complex16* clone_b = new MKL_Complex16[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex16)); - - zgelsy_(&m, &n, &bn, clone_a, &m, clone_b, &m, jpvt, &cond, &rank, work, &len, rwork, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] jpvt; - delete[] rwork; - delete[] clone_a; - delete[] clone_b; - return info; - } - - DLLEXPORT MKL_INT s_thin_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, float a[], float b[], float x[], float work[], MKL_INT len) - { - MKL_INT info = 0; - - float* clone_a = new float[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(float)); - - float* clone_b = new float[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(float)); - char N = 'N'; - sgels_(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_a; - delete[] clone_b; - return info; - } - - DLLEXPORT MKL_INT d_thin_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, double a[], double b[], double x[], double work[], MKL_INT len) - { - MKL_INT info = 0; - - double* clone_a = new double[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(double)); - - double* clone_b = new double[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(double)); - - char N = 'N'; - dgels_(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_a; - delete[] clone_b; - return info; - } - - DLLEXPORT MKL_INT c_thin_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex8 a[], MKL_Complex8 b[], MKL_Complex8 x[], MKL_Complex8 work[], MKL_INT len) - { - MKL_INT info = 0; - MKL_Complex8* clone_a = new MKL_Complex8[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(MKL_Complex8)); - - MKL_Complex8* clone_b = new MKL_Complex8[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex8)); - - char N = 'N'; - cgels_(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_a; - delete[] clone_b; - return info; - } - - DLLEXPORT MKL_INT z_thin_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex16 a[], MKL_Complex16 b[], MKL_Complex16 x[], MKL_Complex16 work[], MKL_INT len) - { - MKL_INT info = 0; - - MKL_Complex16* clone_a = new MKL_Complex16[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(MKL_Complex16)); - - MKL_Complex16* clone_b = new MKL_Complex16[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex16)); - - char N = 'N'; - zgels_(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_a; - delete[] clone_b; - return info; + return qr_solve(m, n, bn, a, b, x, work, len, zgels); } DLLEXPORT MKL_INT s_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, float r[], float b[], float tau[], float x[], float work[], MKL_INT len) { - char side ='L'; - char tran = 'T'; - MKL_INT info = 0; - - float* clone_b = new float[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(float)); - - sormqr_(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_b; - return info; + return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, sormqr, cblas_strsm); } DLLEXPORT MKL_INT d_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, double r[], double b[], double tau[], double x[], double work[], MKL_INT len) { - char side ='L'; - char tran = 'T'; - MKL_INT info = 0; - - double* clone_b = new double[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(double)); - - dormqr_(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_b; - return info; + return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, dormqr, cblas_dtrsm); } DLLEXPORT MKL_INT c_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], MKL_INT len) { - char side ='L'; - char tran = 'C'; - MKL_INT info = 0; - - MKL_Complex8* clone_b = new MKL_Complex8[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex8)); - - cunmqr_(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - MKL_Complex8 one = {1.0f, 0.0f}; - cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_b; - return info; + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, cunmqr, cblas_ctrsm); } DLLEXPORT MKL_INT z_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], MKL_INT len) { - char side ='L'; - char tran = 'C'; - MKL_INT info = 0; - - MKL_Complex16* clone_b = new MKL_Complex16[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex16)); - - zunmqr_(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - MKL_Complex16 one = {1.0, 0.0}; - cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); - - for (MKL_INT i = 0; i < n; ++i) - { - for (MKL_INT j = 0; j < bn; ++j) - { - x[j * n + i] = clone_b[j * m + i]; - } - } - - delete[] clone_b; - return info; + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, zunmqr, cblas_ztrsm); } DLLEXPORT MKL_INT s_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, float a[], float s[], float u[], float v[], float work[], MKL_INT len) { - MKL_INT info = 0; - char job = compute_vectors ? 'A' : 'N'; - sgesvd_(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); - return info; + return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, sgesvd); } DLLEXPORT MKL_INT d_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, double a[], double s[], double u[], double v[], double work[], MKL_INT len) { - MKL_INT info = 0; - char job = compute_vectors ? 'A' : 'N'; - dgesvd_(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); - return info; + return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, dgesvd); } DLLEXPORT MKL_INT c_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, MKL_Complex8 a[], MKL_Complex8 s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT dim_s = std::min(m,n); - float* rwork = new float[5 * dim_s]; - float* s_local = new float[dim_s]; - char job = compute_vectors ? 'A' : 'N'; - cgesvd_(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info); - - for(MKL_INT index = 0; index < dim_s; ++index){ - MKL_Complex8 value = {s_local[index], 0.0f}; - s[index] = value; - } - - delete[] rwork; - delete[] s_local; - return info; + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, cgesvd); } - + DLLEXPORT MKL_INT z_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, MKL_Complex16 a[], MKL_Complex16 s[], MKL_Complex16 u[], MKL_Complex16 v[], MKL_Complex16 work[], MKL_INT len) { - MKL_INT info = 0; - MKL_INT dim_s = std::min(m,n); - double* rwork = new double[5 * std::min(m, n)]; - double* s_local = new double[dim_s]; - char job = compute_vectors ? 'A' : 'N'; - zgesvd_(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info); - - for(MKL_INT index = 0; index < dim_s; ++index){ - MKL_Complex16 value = {s_local[index], 0.0f}; - s[index] = value; - } - - delete[] rwork; - delete[] s_local; - return info; + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, zgesvd); } } \ No newline at end of file diff --git a/src/NativeWrappers/Windows/MKLWrapperTests/MKLWrapperTests.csproj b/src/NativeWrappers/Windows/MKLWrapperTests/MKLWrapperTests.csproj index 4ed8ed5b..579b99cc 100644 --- a/src/NativeWrappers/Windows/MKLWrapperTests/MKLWrapperTests.csproj +++ b/src/NativeWrappers/Windows/MKLWrapperTests/MKLWrapperTests.csproj @@ -78,16 +78,16 @@ AssertHelpers.cs - LinearAlgebra\Complex32\LinearAlgebraProviderTests.cs + LinearAlgebraProviderTests\Complex32\LinearAlgebraProviderTests.cs - LinearAlgebra\Complex\LinearAlgebraProviderTests.cs + LinearAlgebraProviderTests\Complex\LinearAlgebraProviderTests.cs - LinearAlgebra\Double\LinearAlgebraProviderTests.cs + LinearAlgebraProviderTests\Double\LinearAlgebraProviderTests.cs - LinearAlgebra\Single\LinearAlgebraProviderTests.cs + LinearAlgebraProviderTests\Single\LinearAlgebraProviderTests.cs LinearAlgebraTests\Complex32\DenseMatrixTests.cs diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs index 5ac470e4..feaf09cb 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs @@ -740,14 +740,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - if (method == QRMethod.Full) - { - SafeNativeMethods.z_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } - else - { - SafeNativeMethods.z_thin_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + SafeNativeMethods.z_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs index d1a0433c..f1455621 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs @@ -739,14 +739,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - if (method == QRMethod.Full) - { - SafeNativeMethods.c_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } - else - { - SafeNativeMethods.c_thin_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + SafeNativeMethods.c_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.cs index 474f2642..cf43e799 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.cs @@ -843,14 +843,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - if (method == QRMethod.Full) - { - SafeNativeMethods.d_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } - else - { - SafeNativeMethods.d_thin_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + SafeNativeMethods.d_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.cs index 6da7054e..b32d698e 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.cs @@ -743,14 +743,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - if (method == QRMethod.Full) - { - SafeNativeMethods.s_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } - else - { - SafeNativeMethods.s_thin_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); - } + SafeNativeMethods.s_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs index eceb2d2a..21851fe5 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs @@ -242,18 +242,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl [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); - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_thin_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_thin_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_thin_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_thin_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); - [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); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs index 29f6d4cd..93d06e9f 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs @@ -446,7 +446,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization // Check the reconstruction. for (var i = 0; i < order; i++) { - AssertHelpers.AlmostEqual(vectorb[i], matrixBReconstruct[i], 9); + AssertHelpers.AlmostEqual(vectorb[i], matrixBReconstruct[i], 3); } // Make sure A didn't change.