forked from tsai/mathnet-numerics
Christoph Ruegg
11 years ago
50 changed files with 2844 additions and 8745 deletions
@ -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 |
||||
@ -0,0 +1,744 @@ |
|||||
|
#include "wrapper_common.h" |
||||
|
|
||||
|
#include "lapack.h" |
||||
|
#include "lapack_common.h" |
||||
|
#include <algorithm> |
||||
|
|
||||
|
template<typename T, typename GETRF> |
||||
|
inline lapack_int lu_factor(lapack_int m, T a[], lapack_int ipiv[], GETRF getrf) |
||||
|
{ |
||||
|
auto info = getrf(LAPACK_COL_MAJOR, m, m, a, m, ipiv); |
||||
|
shift_ipiv_down(m, ipiv); |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GETRF, typename GETRI> |
||||
|
inline lapack_int lu_inverse(lapack_int n, T a[], GETRF getrf, GETRI getri) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto ipiv = array_new<lapack_int>(n); |
||||
|
auto info = getrf(LAPACK_COL_MAJOR, n, n, a, n, ipiv.get()); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
info = getri(LAPACK_COL_MAJOR, n, a, n, ipiv.get()); |
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GETRI> |
||||
|
inline lapack_int lu_inverse_factored(lapack_int n, T a[], lapack_int ipiv[], GETRI getri) |
||||
|
{ |
||||
|
shift_ipiv_up(n, ipiv); |
||||
|
auto info = getri(LAPACK_COL_MAJOR, n, a, n, ipiv); |
||||
|
shift_ipiv_down(n, ipiv); |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GETRS> |
||||
|
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); |
||||
|
auto info = getrs(LAPACK_COL_MAJOR, 'N', n, nrhs, a, n, ipiv, b, n); |
||||
|
shift_ipiv_down(n, ipiv); |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GETRF, typename GETRS> |
||||
|
inline lapack_int lu_solve(lapack_int n, lapack_int nrhs, T a[], T b[], GETRF getrf, GETRS getrs) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto clone = array_clone(n * n, a); |
||||
|
auto ipiv = array_new<lapack_int>(n); |
||||
|
auto info = getrf(LAPACK_COL_MAJOR, n, n, clone.get(), n, ipiv.get()); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
return getrs(LAPACK_COL_MAJOR, 'N', n, nrhs, clone.get(), n, ipiv.get(), b, n); |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename POTRF> |
||||
|
inline lapack_int cholesky_factor(lapack_int n, T* a, POTRF potrf) |
||||
|
{ |
||||
|
auto info = potrf(LAPACK_COL_MAJOR, 'L', n, a, n); |
||||
|
auto zero = T(); |
||||
|
|
||||
|
for (auto i = 0; i < n; ++i) |
||||
|
{ |
||||
|
auto index = i * n; |
||||
|
|
||||
|
for (auto j = 0; j < n && i > j; ++j) |
||||
|
{ |
||||
|
a[index + j] = zero; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename POTRF, typename POTRS> |
||||
|
inline lapack_int cholesky_solve(lapack_int n, lapack_int nrhs, T a[], T b[], POTRF potrf, POTRS potrs) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto clone = array_clone(n * n, a); |
||||
|
auto info = potrf(LAPACK_COL_MAJOR, 'L', n, clone.get(), n); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
return potrs(LAPACK_COL_MAJOR, 'L', n, nrhs, clone.get(), n, b, n); |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
|
||||
|
template<typename T, typename GEQRF, typename ORGQR> |
||||
|
inline lapack_int qr_factor(lapack_int m, lapack_int n, T r[], T tau[], T q[], GEQRF geqrf, ORGQR orgqr) |
||||
|
{ |
||||
|
auto info = geqrf(LAPACK_COL_MAJOR, m, n, r, m, tau); |
||||
|
|
||||
|
for (auto i = 0; i < m; ++i) |
||||
|
{ |
||||
|
for (auto j = 0; j < m && j < n; ++j) |
||||
|
{ |
||||
|
if (i > j) |
||||
|
{ |
||||
|
q[j * m + i] = r[j * m + i]; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
//compute the q elements explicitly
|
||||
|
if (m <= n) |
||||
|
{ |
||||
|
info = orgqr(LAPACK_COL_MAJOR, m, m, m, q, m, tau); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
info = orgqr(LAPACK_COL_MAJOR, m, m, n, q, m, tau); |
||||
|
} |
||||
|
|
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GEQRF, typename ORGQR> |
||||
|
inline lapack_int qr_thin_factor(lapack_int m, lapack_int n, T q[], T tau[], T r[], GEQRF geqrf, ORGQR orgqr) |
||||
|
{ |
||||
|
auto info = geqrf(LAPACK_COL_MAJOR, m, n, q, m, tau); |
||||
|
|
||||
|
for (auto i = 0; i < n; ++i) |
||||
|
{ |
||||
|
for (auto j = 0; j < n; ++j) |
||||
|
{ |
||||
|
if (i <= j) |
||||
|
{ |
||||
|
r[j * n + i] = q[j * m + i]; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
info = orgqr(LAPACK_COL_MAJOR, m, n, n, q, m, tau); |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GELS> |
||||
|
inline lapack_int qr_solve(lapack_int m, lapack_int n, lapack_int bn, T a[], T b[], T x[], GELS gels) |
||||
|
{ |
||||
|
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); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
copyBtoX(m, n, bn, clone_b.get(), x); |
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename ORMQR, typename TRSM> |
||||
|
inline lapack_int qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], ORMQR ormqr, TRSM trsm) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto clone_b = array_clone(m * bn, b); |
||||
|
auto info = ormqr(LAPACK_COL_MAJOR, 'L', 'T', m, bn, n, r, m, tau, clone_b.get(), m); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b.get(), m); |
||||
|
copyBtoX(m, n, bn, clone_b.get(), x); |
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename R, typename UNMQR, typename TRSM> |
||||
|
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 |
||||
|
{ |
||||
|
auto clone_b = array_clone(m * bn, b); |
||||
|
auto info = unmqr(LAPACK_COL_MAJOR, 'L', 'C', m, bn, n, r, m, tau, clone_b.get(), m); |
||||
|
|
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
T one = 1.0f; |
||||
|
trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, reinterpret_cast<R*>(&one), reinterpret_cast<R*>(r), m, reinterpret_cast<R*>(clone_b.get()), m); |
||||
|
copyBtoX(m, n, bn, clone_b.get(), x); |
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GESVD> |
||||
|
inline lapack_int svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], GESVD gesvd) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto job = compute_vectors ? 'A' : 'N'; |
||||
|
auto dim_s = std::min(m, n); |
||||
|
auto superb = array_new<T>(std::max(2, dim_s) - 1); |
||||
|
return gesvd(LAPACK_COL_MAJOR, job, job, m, n, a, m, s, u, m, v, n, superb.get()); |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename R, typename GESVD> |
||||
|
inline lapack_int complex_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], GESVD gesvd) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto dim_s = std::min(m, n); |
||||
|
auto s_local = array_new<R>(dim_s); |
||||
|
auto superb = array_new<R>(std::max(2, dim_s) - 1); |
||||
|
auto job = compute_vectors ? 'A' : 'N'; |
||||
|
auto info = gesvd(LAPACK_COL_MAJOR, job, job, m, n, a, m, s_local.get(), u, m, v, n, superb.get()); |
||||
|
|
||||
|
for (auto index = 0; index < dim_s; ++index) |
||||
|
{ |
||||
|
s[index] = s_local.get()[index]; |
||||
|
} |
||||
|
|
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename R, typename GEES, typename TREVC> |
||||
|
inline lapack_int eigen_factor(lapack_int n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto clone_a = array_clone(n * n, a); |
||||
|
auto wr = array_new<T>(n); |
||||
|
auto wi = array_new<T>(n); |
||||
|
|
||||
|
lapack_int sdim; |
||||
|
lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a.get(), n, &sdim, wr.get(), wi.get(), vectors, n); |
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
lapack_int m; |
||||
|
info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a.get(), n, nullptr, n, vectors, n, n, &m); |
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
for (auto index = 0; index < n; ++index) |
||||
|
{ |
||||
|
values[index] = R(wr.get()[index], wi.get()[index]); |
||||
|
} |
||||
|
|
||||
|
for (auto i = 0; i < n; ++i) |
||||
|
{ |
||||
|
auto in = i * n; |
||||
|
d[in + i] = wr.get()[i]; |
||||
|
|
||||
|
if (wi.get()[i] > 0) |
||||
|
{ |
||||
|
d[in + n + i] = wi.get()[i]; |
||||
|
} |
||||
|
else if (wi.get()[i] < 0) |
||||
|
{ |
||||
|
d[in - n + i] = wi.get()[i]; |
||||
|
} |
||||
|
} |
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename T, typename GEES, typename 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 |
||||
|
{ |
||||
|
auto clone_a = array_clone(n * n, a); |
||||
|
auto w = array_new<T>(n); |
||||
|
|
||||
|
lapack_int sdim; |
||||
|
lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a.get(), n, &sdim, w.get(), vectors, n); |
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
lapack_int m; |
||||
|
info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a.get(), n, nullptr, n, vectors, n, n, &m); |
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
for (auto i = 0; i < n; ++i) |
||||
|
{ |
||||
|
values[i] = w.get()[i]; |
||||
|
d[i * n + i] = w.get()[i]; |
||||
|
} |
||||
|
|
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template<typename R, typename T, typename SYEV> |
||||
|
inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], SYEV syev) |
||||
|
{ |
||||
|
try |
||||
|
{ |
||||
|
auto clone_a = array_clone(n * n, a); |
||||
|
auto w = array_new<R>(n); |
||||
|
|
||||
|
lapack_int info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a.get(), n, w.get()); |
||||
|
if (info != 0) |
||||
|
{ |
||||
|
return info; |
||||
|
} |
||||
|
|
||||
|
memcpy(vectors, clone_a.get(), n*n*sizeof(T)); |
||||
|
|
||||
|
for (auto index = 0; index < n; ++index) |
||||
|
{ |
||||
|
values[index] = lapack_complex_double(w.get()[index]); |
||||
|
} |
||||
|
|
||||
|
for (auto j = 0; j < n; ++j) |
||||
|
{ |
||||
|
auto jn = j*n; |
||||
|
|
||||
|
for (auto i = 0; i < n; ++i) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
d[jn + i] = w.get()[i]; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return info; |
||||
|
} |
||||
|
catch (std::bad_alloc&) |
||||
|
{ |
||||
|
return INSUFFICIENT_MEMORY; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
extern "C" { |
||||
|
|
||||
|
DLLEXPORT float s_matrix_norm(char norm, lapack_int m, lapack_int n, float a[]) |
||||
|
{ |
||||
|
return LAPACKE_slange(LAPACK_COL_MAJOR, norm, m, n, a, m); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT double d_matrix_norm(char norm, lapack_int m, lapack_int n, double a[]) |
||||
|
{ |
||||
|
return LAPACKE_dlange(LAPACK_COL_MAJOR, norm, m, n, a, m); |
||||
|
} |
||||
|
|
||||
|
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, lapack_complex_double a[]) |
||||
|
{ |
||||
|
return LAPACKE_zlange(LAPACK_COL_MAJOR, norm, m, n, a, m); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_lu_factor(lapack_int m, float a[], lapack_int ipiv[]) |
||||
|
{ |
||||
|
return lu_factor(m, a, ipiv, LAPACKE_sgetrf); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_lu_factor(lapack_int m, double a[], lapack_int ipiv[]) |
||||
|
{ |
||||
|
return lu_factor(m, a, ipiv, LAPACKE_dgetrf); |
||||
|
} |
||||
|
|
||||
|
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, lapack_complex_double a[], lapack_int ipiv[]) |
||||
|
{ |
||||
|
return lu_factor(m, a, ipiv, LAPACKE_zgetrf); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_lu_inverse(lapack_int n, float a[], float work[], lapack_int lwork) |
||||
|
{ |
||||
|
return lu_inverse(n, a, LAPACKE_sgetrf, LAPACKE_sgetri); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_lu_inverse(lapack_int n, double a[], double work[], lapack_int lwork) |
||||
|
{ |
||||
|
return lu_inverse(n, a, LAPACKE_dgetrf, LAPACKE_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, LAPACKE_cgetrf, LAPACKE_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, LAPACKE_zgetrf, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_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, LAPACKE_zgetrs); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_lu_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) |
||||
|
{ |
||||
|
return lu_solve(n, nrhs, a, b, LAPACKE_sgetrf, LAPACKE_sgetrs); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_lu_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) |
||||
|
{ |
||||
|
return lu_solve(n, nrhs, a, b, LAPACKE_dgetrf, LAPACKE_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, LAPACKE_cgetrf, LAPACKE_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, LAPACKE_zgetrf, LAPACKE_zgetrs); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_cholesky_factor(lapack_int n, float a[]) |
||||
|
{ |
||||
|
return cholesky_factor(n, a, LAPACKE_spotrf); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_cholesky_factor(lapack_int n, double* a) |
||||
|
{ |
||||
|
return cholesky_factor(n, a, LAPACKE_dpotrf); |
||||
|
} |
||||
|
|
||||
|
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, lapack_complex_double a[]) |
||||
|
{ |
||||
|
return cholesky_factor(n, a, LAPACKE_zpotrf); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_cholesky_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) |
||||
|
{ |
||||
|
return cholesky_solve(n, nrhs, a, b, LAPACKE_spotrf, LAPACKE_spotrs); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_cholesky_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) |
||||
|
{ |
||||
|
return cholesky_solve(n, nrhs, a, b, LAPACKE_dpotrf, LAPACKE_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, LAPACKE_cpotrf, LAPACKE_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, LAPACKE_zpotrf, LAPACKE_zpotrs); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_cholesky_solve_factored(lapack_int n, lapack_int nrhs, float a[], float b[]) |
||||
|
{ |
||||
|
return LAPACKE_spotrs(LAPACK_COL_MAJOR, 'L', n, nrhs, a, n, b, n); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_cholesky_solve_factored(lapack_int n, lapack_int nrhs, double a[], double b[]) |
||||
|
{ |
||||
|
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, 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, lapack_complex_double a[], lapack_complex_double b[]) |
||||
|
{ |
||||
|
return LAPACKE_zpotrs(LAPACK_COL_MAJOR, 'L', n, nrhs, a, n, b, n); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_qr_factor(lapack_int m, lapack_int n, float r[], float tau[], float q[]) |
||||
|
{ |
||||
|
return qr_factor(m, n, r, tau, q, LAPACKE_sgeqrf, LAPACKE_sorgqr); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_qr_thin_factor(lapack_int m, lapack_int n, float q[], float tau[], float r[]) |
||||
|
{ |
||||
|
return qr_thin_factor(m, n, q, tau, r, LAPACKE_sgeqrf, LAPACKE_sorgqr); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_qr_factor(lapack_int m, lapack_int n, double r[], double tau[], double q[]) |
||||
|
{ |
||||
|
return qr_factor(m, n, r, tau, q, LAPACKE_dgeqrf, LAPACKE_dorgqr); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_qr_thin_factor(lapack_int m, lapack_int n, double q[], double tau[], double r[]) |
||||
|
{ |
||||
|
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, 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, 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, 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, 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); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int s_qr_solve(lapack_int m, lapack_int n, lapack_int bn, float a[], float b[], float x[]) |
||||
|
{ |
||||
|
return qr_solve(m, n, bn, a, b, x, LAPACKE_sgels); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_qr_solve(lapack_int m, lapack_int n, lapack_int bn, double a[], double b[], double x[]) |
||||
|
{ |
||||
|
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, 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, lapack_complex_double a[], lapack_complex_double b[], lapack_complex_double x[]) |
||||
|
{ |
||||
|
return qr_solve(m, n, bn, a, b, x, LAPACKE_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[]) |
||||
|
{ |
||||
|
return qr_solve_factored(m, n, bn, r, b, tau, x, LAPACKE_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[]) |
||||
|
{ |
||||
|
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, lapack_complex_float r[], lapack_complex_float b[], lapack_complex_float tau[], lapack_complex_float x[]) |
||||
|
{ |
||||
|
return complex_qr_solve_factored<lapack_complex_float, float>(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, lapack_complex_double r[], lapack_complex_double b[], lapack_complex_double tau[], lapack_complex_double x[]) |
||||
|
{ |
||||
|
return complex_qr_solve_factored<lapack_complex_double, double>(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[]) |
||||
|
{ |
||||
|
return svd_factor(compute_vectors, m, n, a, s, u, v, LAPACKE_sgesvd); |
||||
|
} |
||||
|
|
||||
|
DLLEXPORT lapack_int d_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, double a[], double s[], double u[], double v[]) |
||||
|
{ |
||||
|
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, lapack_complex_float a[], lapack_complex_float s[], lapack_complex_float u[], lapack_complex_float v[]) |
||||
|
{ |
||||
|
return complex_svd_factor<lapack_complex_float, float>(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, lapack_complex_double a[], lapack_complex_double s[], lapack_complex_double u[], lapack_complex_double v[]) |
||||
|
{ |
||||
|
return complex_svd_factor<lapack_complex_double, double>(compute_vectors, m, n, a, s, u, v, LAPACKE_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<float>(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<double>(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<float>(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<double>(n, a, vectors, values, d, LAPACKE_zheev); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -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 |
|
||||
@ -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 |
||||
@ -1,713 +0,0 @@ |
|||||
#include <algorithm> |
|
||||
#include <complex> |
|
||||
|
|
||||
#define MKL_Complex8 std::complex<float> |
|
||||
#define MKL_Complex16 std::complex<double> |
|
||||
|
|
||||
#include "mkl_lapack.h" |
|
||||
#include "mkl_cblas.h" |
|
||||
#include "lapack_common.h" |
|
||||
#include "wrapper_common.h" |
|
||||
#include "mkl_lapacke.h" |
|
||||
#include "mkl.h" |
|
||||
#include "mkl_trans.h" |
|
||||
|
|
||||
template<typename T, typename GETRF> |
|
||||
inline MKL_INT lu_factor(MKL_INT m, T a[], MKL_INT ipiv[], GETRF getrf) |
|
||||
{ |
|
||||
MKL_INT info = 0; |
|
||||
getrf(&m, &m, a, &m, ipiv, &info); |
|
||||
shift_ipiv_down(m, ipiv); |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GETRF, typename GETRI> |
|
||||
inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork, GETRF getrf, GETRI getri) |
|
||||
{ |
|
||||
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; |
|
||||
} |
|
||||
|
|
||||
getri(&n, a, &n, ipiv, work, &lwork, &info); |
|
||||
delete[] ipiv; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GETRI> |
|
||||
inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], MKL_INT lwork, GETRI getri) |
|
||||
{ |
|
||||
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<typename T, typename GETRS> |
|
||||
inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[], T b[], GETRS getrs) |
|
||||
{ |
|
||||
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<typename T, typename GETRF, typename GETRS> |
|
||||
inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], GETRF getrf, GETRS getrs) |
|
||||
{ |
|
||||
T* clone = Clone(n, n, a); |
|
||||
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<typename T, typename POTRF> |
|
||||
inline MKL_INT cholesky_factor(MKL_INT n, T* a, POTRF potrf) |
|
||||
{ |
|
||||
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 index = i * n; |
|
||||
|
|
||||
for (MKL_INT j = 0; j < n && i > j; ++j) |
|
||||
{ |
|
||||
a[index + j] = zero; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename POTRF, typename POTRS> |
|
||||
inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], POTRF potrf, POTRS potrs) |
|
||||
{ |
|
||||
T* clone = Clone(n, n, a); |
|
||||
char uplo = 'L'; |
|
||||
MKL_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<typename T, typename POTRS> |
|
||||
inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[], POTRS potrs) |
|
||||
{ |
|
||||
char uplo = 'L'; |
|
||||
MKL_INT info = 0; |
|
||||
potrs(&uplo, &n, &nrhs, a, &n, b, &n, &info); |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GEQRF, typename ORGQR> |
|
||||
inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], MKL_INT len, GEQRF geqrf, ORGQR orgqr) |
|
||||
{ |
|
||||
MKL_INT info = 0; |
|
||||
geqrf(&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) |
|
||||
{ |
|
||||
orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); |
|
||||
} |
|
||||
|
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GEQRF, typename ORGQR> |
|
||||
inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T work[], MKL_INT len, GEQRF geqrf, ORGQR orgqr) |
|
||||
{ |
|
||||
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]; |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GELS> |
|
||||
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, GELS gels) |
|
||||
{ |
|
||||
T* clone_a = Clone(m, n, a); |
|
||||
T* clone_b = Clone(m, bn, b); |
|
||||
char N = 'N'; |
|
||||
MKL_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<typename T, typename ORMQR, typename TRSM> |
|
||||
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, ORMQR ormqr, TRSM trsm) |
|
||||
{ |
|
||||
T* clone_b = Clone(m, bn, b); |
|
||||
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); |
|
||||
copyBtoX(m, n, bn, clone_b, x); |
|
||||
delete[] clone_b; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename UNMQR, typename TRSM> |
|
||||
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, UNMQR unmqr, TRSM trsm) |
|
||||
{ |
|
||||
T* clone_b = Clone(m, bn, b); |
|
||||
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; |
|
||||
trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); |
|
||||
copyBtoX(m, n, bn, clone_b, x); |
|
||||
delete[] clone_b; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GESVD> |
|
||||
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, GESVD gesvd) |
|
||||
{ |
|
||||
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<typename T, typename R, typename GESVD> |
|
||||
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, GESVD gesvd) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
s[index] = s_local[index]; |
|
||||
} |
|
||||
|
|
||||
delete[] rwork; |
|
||||
delete[] s_local; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename R, typename GEES, typename TREVC> |
|
||||
inline MKL_INT eigen_factor(MKL_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]; |
|
||||
|
|
||||
MKL_INT sdim; |
|
||||
MKL_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; |
|
||||
} |
|
||||
|
|
||||
MKL_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 (MKL_INT index = 0; index < n; ++index) |
|
||||
{ |
|
||||
values[index] = R(wr[index], wi[index]); |
|
||||
} |
|
||||
|
|
||||
for (MKL_INT i = 0; i < n; ++i) |
|
||||
{ |
|
||||
MKL_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<typename T, typename GEES, typename TREVC> |
|
||||
inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], GEES gees, TREVC trevc) |
|
||||
{ |
|
||||
T* clone_a = Clone(n, n, a); |
|
||||
T* w = new T[n]; |
|
||||
|
|
||||
MKL_INT sdim; |
|
||||
MKL_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; |
|
||||
} |
|
||||
|
|
||||
MKL_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 (MKL_INT i = 0; i < n; ++i) |
|
||||
{ |
|
||||
values[i] = w[i]; |
|
||||
d[i * n + i] = w[i]; |
|
||||
} |
|
||||
|
|
||||
delete[] clone_a; |
|
||||
delete[] w; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename R, typename T, typename SYEV> |
|
||||
inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], SYEV syev) |
|
||||
{ |
|
||||
T* clone_a = Clone(n, n, a); |
|
||||
R* w = new R[n]; |
|
||||
|
|
||||
MKL_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 (MKL_INT index = 0; index < n; ++index) |
|
||||
{ |
|
||||
values[index] = MKL_Complex16(w[index]); |
|
||||
} |
|
||||
|
|
||||
for (MKL_INT j = 0; j < n; ++j) |
|
||||
{ |
|
||||
MKL_INT jn = j*n; |
|
||||
|
|
||||
for (MKL_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, MKL_INT m, MKL_INT n, float a[], float work[]) |
|
||||
{ |
|
||||
return slange(&norm, &m, &n, a, &m, work); |
|
||||
} |
|
||||
|
|
||||
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); |
|
||||
} |
|
||||
|
|
||||
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 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); |
|
||||
} |
|
||||
|
|
||||
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); |
|
||||
} |
|
||||
|
|
||||
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 d_lu_inverse(MKL_INT n, double a[], double work[], MKL_INT lwork) |
|
||||
{ |
|
||||
return lu_inverse(n, a, work, lwork, dgetrf, dgetri); |
|
||||
} |
|
||||
|
|
||||
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); |
|
||||
} |
|
||||
|
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
return lu_solve(n, nrhs, a, b, zgetrf, zgetrs); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[]) |
|
||||
{ |
|
||||
return cholesky_factor(n, a, spotrf); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a) |
|
||||
{ |
|
||||
return cholesky_factor(n, a, dpotrf); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[]) |
|
||||
{ |
|
||||
return cholesky_factor(n, a, cpotrf); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[]) |
|
||||
{ |
|
||||
return cholesky_factor(n, a, zpotrf); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT s_cholesky_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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[]) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
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) |
|
||||
{ |
|
||||
return complex_svd_factor<MKL_Complex8, float>(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) |
|
||||
{ |
|
||||
return complex_svd_factor<MKL_Complex16, double>(compute_vectors, m, n, a, s, u, v, work, len, zgesvd); |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT s_eigen(bool isSymmetric, MKL_INT n, float a[], float vectors[], MKL_Complex16 values[], float d[]) |
|
||||
{ |
|
||||
if (isSymmetric) |
|
||||
{ |
|
||||
return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_ssyev); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT d_eigen(bool isSymmetric, MKL_INT n, double a[], double vectors[], MKL_Complex16 values[], double d[]) |
|
||||
{ |
|
||||
if (isSymmetric) |
|
||||
{ |
|
||||
return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_dsyev); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT c_eigen(bool isSymmetric, MKL_INT n, MKL_Complex8 a[], MKL_Complex8 vectors[], MKL_Complex16 values[], MKL_Complex8 d[]) |
|
||||
{ |
|
||||
if (isSymmetric) |
|
||||
{ |
|
||||
return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_cheev); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
DLLEXPORT MKL_INT z_eigen(bool isSymmetric, MKL_INT n, MKL_Complex16 a[], MKL_Complex16 vectors[], MKL_Complex16 values[], MKL_Complex16 d[]) |
|
||||
{ |
|
||||
if (isSymmetric) |
|
||||
{ |
|
||||
return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_zheev); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
@ -0,0 +1,32 @@ |
|||||
|
#pragma once |
||||
|
|
||||
|
#include <complex> |
||||
|
#include <memory> |
||||
|
|
||||
|
#define MKL_Complex8 std::complex<float> |
||||
|
#define MKL_Complex16 std::complex<double> |
||||
|
#define LAPACK_MEMORY |
||||
|
|
||||
|
#include "mkl.h" |
||||
|
|
||||
|
const int ALIGNMENT = 64; |
||||
|
|
||||
|
struct array_free |
||||
|
{ |
||||
|
void operator()(void* x) { mkl_free(x); } |
||||
|
}; |
||||
|
|
||||
|
template <typename T> using array_ptr = std::unique_ptr<T[], array_free>; |
||||
|
|
||||
|
template<typename T> |
||||
|
inline array_ptr<T> array_new(const int size) |
||||
|
{ |
||||
|
auto ret = static_cast<T*>(mkl_malloc(size * sizeof(T), ALIGNMENT)); |
||||
|
|
||||
|
if (!ret) |
||||
|
{ |
||||
|
throw new std::bad_alloc(); |
||||
|
} |
||||
|
|
||||
|
return array_ptr<T>(ret); |
||||
|
} |
||||
@ -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 |
|
||||
@ -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 |
||||
@ -1,39 +0,0 @@ |
|||||
template <typename _T> |
|
||||
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<typename _Other> inline |
|
||||
complex& operator=(const complex<_Other>& right) |
|
||||
{ |
|
||||
real = (_T)right.real; |
|
||||
imag = (_T)right.imag; |
|
||||
return *this; |
|
||||
} |
|
||||
}; |
|
||||
@ -1,712 +0,0 @@ |
|||||
#include "cblas.h" |
|
||||
|
|
||||
#include "complex.h" |
|
||||
#define LAPACK_COMPLEX_CUSTOM |
|
||||
#define lapack_complex_float complex<float> |
|
||||
#define lapack_complex_double complex<double> |
|
||||
|
|
||||
#include "lapacke.h" |
|
||||
#include "lapack_common.h" |
|
||||
#include "wrapper_common.h" |
|
||||
#include <algorithm> |
|
||||
|
|
||||
template<typename T, typename GETRF> |
|
||||
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<typename T, typename GETRF, typename GETRI> |
|
||||
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<typename T, typename GETRI> |
|
||||
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<typename T, typename GETRS> |
|
||||
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<typename T, typename GETRF, typename GETRS> |
|
||||
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<typename T, typename POTRF> |
|
||||
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<typename T, typename POTRF, typename POTRS> |
|
||||
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<typename T, typename POTRS> |
|
||||
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<typename T, typename GEQRF, typename ORGQR> |
|
||||
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<typename T, typename GEQRF, typename ORGQR> |
|
||||
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<typename T, typename GELS> |
|
||||
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<typename T, typename ORMQR, typename TRSM> |
|
||||
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<typename T, typename R, typename UNMQR, typename TRSM> |
|
||||
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<R*>(&one), reinterpret_cast<R*>(r), m, reinterpret_cast<R*>(clone_b), m); |
|
||||
copyBtoX(m, n, bn, clone_b, x); |
|
||||
delete[] clone_b; |
|
||||
return info; |
|
||||
} |
|
||||
|
|
||||
template<typename T, typename GESVD> |
|
||||
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<typename T, typename R, typename GESVD> |
|
||||
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<typename T, typename R, typename GEES, typename TREVC> |
|
||||
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<typename T, typename GEES, typename TREVC> |
|
||||
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<typename R, typename T, typename SYEV> |
|
||||
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<lapack_complex_float, float>(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<lapack_complex_double, double>(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<lapack_complex_float, float>(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<lapack_complex_double, double>(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<float>(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<double>(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<float>(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<double>(n, a, vectors, values, d, LAPACKE_zheev); |
|
||||
} |
|
||||
else |
|
||||
{ |
|
||||
return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
@ -0,0 +1,11 @@ |
|||||
|
#pragma once |
||||
|
|
||||
|
#define LAPACK_COMPLEX_CUSTOM |
||||
|
#include <complex> |
||||
|
#define lapack_complex_float std::complex<float> |
||||
|
#define lapack_complex_double std::complex<double> |
||||
|
|
||||
|
#include "cblas.h" |
||||
|
#include "lapacke.h" |
||||
|
|
||||
|
|
||||
Loading…
Reference in new issue