ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Math.NET Numerics
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1478 Commits
25 Branches
0 Tags
187 MiB
 
 
 
Tree: 1fb3357b80
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '1fb3357b80'
${ noResults }
mathnet-numerics/src/NativeProviders/MKL/lapack.cpp

760 lines
26 KiB
Raw Blame History

#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>
inline MKL_INT lu_factor(MKL_INT m, T a[], MKL_INT ipiv[],
void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*))
{
std::complex<double> x = 5;
MKL_INT info = 0;
getrf(&m, &m, a, &m, ipiv, &info);
shift_ipiv_down(m, ipiv);
return info;
};
template<typename T>
inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork,
void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*),
void (*getri)(const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
MKL_INT* ipiv = new MKL_INT[n];
MKL_INT info = 0;
getrf(&n, &n, a, &n, ipiv, &info);
if (info != 0)
{
delete[] ipiv;
return info;
}
getri(&n, a, &n, ipiv, work, &lwork, &info);
delete[] ipiv;
return info;
};
template<typename T>
inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], MKL_INT lwork,
void (*getri)(const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
shift_ipiv_up(n, ipiv);
MKL_INT info = 0;
getri(&n, a, &n, ipiv, work, &lwork, &info);
shift_ipiv_down(n, ipiv);
return info;
}
template<typename T>
inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[], T b[],
void (*getrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
shift_ipiv_up(n, ipiv);
MKL_INT info = 0;
char trans ='N';
getrs(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
shift_ipiv_down(n, ipiv);
return info;
}
template<typename T>
inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*),
void (*getrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
T* clone = 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>
inline MKL_INT cholesky_factor(MKL_INT n, T* a,
void (*potrf)(const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
char uplo = 'L';
MKL_INT info = 0;
potrf(&uplo, &n, a, &n, &info);
T zero = T();
for (MKL_INT i = 0; i < n; ++i)
{
MKL_INT index = i * n;
for (MKL_INT j = 0; j < n && i > j; ++j)
{
a[index + j] = zero;
}
}
return info;
}
template<typename T>
inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
void (*potrf)(const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*),
void (*potrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
T* clone = 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>
inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[],
void (*potrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
char uplo = 'L';
MKL_INT info = 0;
potrs(&uplo, &n, &nrhs, a, &n, b, &n, &info);
return info;
}
template<typename T>
inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], MKL_INT len,
void (*geqrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*),
void (*orgqr)(const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*))
{
MKL_INT info = 0;
geqrf(&m, &n, r, &m, tau, work, &len, &info);
for (MKL_INT i = 0; i < m; ++i)
{
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>
inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T work[], MKL_INT len,
void (*geqrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*),
void (*orgqr)(const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*))
{
MKL_INT info = 0;
geqrf(&m, &n, q, &m, tau, work, &len, &info);
for (MKL_INT i = 0; i < n; ++i)
{
for (MKL_INT j = 0; j < n; ++j)
{
if (i <= j)
{
r[j * n + i] = q[j * m + i];
}
}
}
orgqr(&m, &n, &n, q, &m, tau, work, &len, &info);
return info;
}
template<typename T>
inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T work[], MKL_INT len,
void (*gels)(const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*, T*,
const MKL_INT*, T* b, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
T* clone_a = 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>
inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len,
void (*ormqr)(const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*,
const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info),
void (*trsm)(const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG,
const MKL_INT, const MKL_INT, const T, const T*, const MKL_INT, T*, const MKL_INT))
{
T* clone_b = 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>
inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len,
void (*unmqr)(const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*,
const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info),
void (*trsm)(const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG,
const MKL_INT, const MKL_INT, const void*, const void*, const MKL_INT, void*, const MKL_INT ldb))
{
T* clone_b = 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>
inline MKL_INT svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len,
void (*gesvd)(const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*,
T*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
{
MKL_INT info = 0;
char job = compute_vectors ? 'A' : 'N';
gesvd(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
return info;
}
template<typename T, typename R>
inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len,
void (*gesvd)(const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*,
R*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, R*, MKL_INT*))
{
MKL_INT info = 0;
MKL_INT dim_s = std::min(m,n);
R* rwork = new R[5 * dim_s];
R* s_local = new R[dim_s];
char job = compute_vectors ? 'A' : 'N';
gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info);
for (MKL_INT index = 0; index < dim_s; ++index)
{
s[index] = s_local[index];
}
delete[] rwork;
delete[] s_local;
return info;
}
template<typename T>
inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
MKL_INT(*gees)(MKL_INT, char, char, int(*)(const T*, const T*), MKL_INT, T* a,
MKL_INT, MKL_INT*, T*, T*, T*, MKL_INT),
MKL_INT(*trevc)(MKL_INT, char, char, lapack_logical*, MKL_INT, const T*,
MKL_INT, T*, MKL_INT, T*, MKL_INT, MKL_INT, MKL_INT*))
{
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] = MKL_Complex16(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>
inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
MKL_INT(*gees)(MKL_INT, char, char, int(*)(const T*), MKL_INT, T* a,
MKL_INT, MKL_INT*, T*, T*, MKL_INT),
MKL_INT(*trevc)(MKL_INT, char, char, const lapack_logical*, MKL_INT, T*,
MKL_INT, T*, MKL_INT, T*, MKL_INT, MKL_INT, MKL_INT*))
{
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 T, typename R>
inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
MKL_INT(*syev)(int, char, char, int, T*, int, R*))
{
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<float>(m, a, ipiv, sgetrf);
}
DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[])
{
return lu_factor<double>(m, a, ipiv, dgetrf);
}
DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[])
{
return lu_factor<MKL_Complex8>(m, a, ipiv, cgetrf);
}
DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[])
{
return lu_factor<MKL_Complex16>(m, a, ipiv, zgetrf);
}
DLLEXPORT MKL_INT s_lu_inverse(MKL_INT n, float a[], float work[], MKL_INT lwork)
{
return lu_inverse<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(n, nrhs, a, b, zgetrf, zgetrs);
}
DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[])
{
return cholesky_factor<float>(n, a, spotrf);
}
DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a)
{
return cholesky_factor<double>(n, a, dpotrf);
}
DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[])
{
return cholesky_factor<MKL_Complex8>(n, a, cpotrf);
}
DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[])
{
return cholesky_factor<MKL_Complex16>(n, a, zpotrf);
}
DLLEXPORT MKL_INT s_cholesky_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[])
{
return cholesky_solve<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<float>(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<double>(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<double>(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<MKL_Complex8>(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<MKL_Complex8>(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<MKL_Complex16>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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<MKL_Complex8>(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<MKL_Complex16>(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<float>(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<double>(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, float>(n, a, vectors, values, d, LAPACKE_ssyev);
}
else
{
return eigen_factor<float>(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, double>(n, a, vectors, values, d, LAPACKE_dsyev);
}
else
{
return eigen_factor<double>(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<MKL_Complex8, float>(n, a, vectors, values, d, LAPACKE_cheev);
}
else
{
return -1;
//return eigen_factor<MKL_Complex16, LAPACK_Z_SELECT1>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
}
}
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<MKL_Complex16, double>(n, a, vectors, values, d, LAPACKE_zheev);
}
else
{
return eigen_complex_factor<MKL_Complex16>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
}
}
DLLEXPORT void SetImprovedConsistency(void)
{
// set improved consistency for mkl and vector functions
mkl_cbwr_set(MKL_CBWR_COMPATIBLE);
vmlSetMode(VML_HA|VML_DOUBLE_CONSISTENT);
}
}