forked from tsai/mathnet-numerics
10 changed files with 220 additions and 11 deletions
@ -1,4 +1,9 @@ |
|||
#ifndef BLAS_H |
|||
#define BLAS_H |
|||
|
|||
#include "cblas.h" |
|||
|
|||
typedef struct { float real; float imag; } Complex8; |
|||
typedef struct { double real; double imag; } Complex16; |
|||
|
|||
#endif |
|||
|
|||
@ -0,0 +1,64 @@ |
|||
#include "common.h" |
|||
#include "blas.h" |
|||
#include "clapack.h" |
|||
|
|||
extern "C" { |
|||
|
|||
DLLEXPORT int s_cholesky_factor(int n, float a[]){ |
|||
int info = clapack_spotrf(CblasColMajor, CblasLower, n, a, n); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = 0; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int d_cholesky_factor(int n, double* a){ |
|||
int info = clapack_dpotrf(CblasColMajor, CblasLower, n, a, n); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = 0; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ |
|||
int info = clapack_cpotrf(CblasColMajor, CblasLower, n, a, n); |
|||
Complex8 zero; |
|||
zero.real = 0.0; |
|||
zero.real = 0.0; |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = zero; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ |
|||
int info = clapack_zpotrf(CblasColMajor, CblasLower, n, a, n); |
|||
Complex16 zero; |
|||
zero.real = 0.0; |
|||
zero.real = 0.0; |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = zero; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
} |
|||
@ -0,0 +1,10 @@ |
|||
#ifndef COMMON_H |
|||
#define COMMON_H |
|||
|
|||
#ifdef _WINDOWS |
|||
#define DLLEXPORT __declspec( dllexport ) |
|||
#else |
|||
#define DLLEXPORT |
|||
#endif |
|||
|
|||
#endif |
|||
@ -1,4 +1,9 @@ |
|||
#ifndef BLAS_H |
|||
#define BLAS_H |
|||
|
|||
#include "mkl_cblas.h" |
|||
|
|||
typedef MKL_Complex8 Complex8; |
|||
typedef MKL_Complex16 Complex16; |
|||
typedef MKL_Complex16 Complex16; |
|||
|
|||
#endif |
|||
|
|||
@ -0,0 +1,72 @@ |
|||
#include "common.h" |
|||
#include "blas.h" |
|||
#include "mkl_lapack.h" |
|||
|
|||
extern "C" { |
|||
|
|||
DLLEXPORT int s_cholesky_factor(int n, float a[]){ |
|||
char uplo = 'L'; |
|||
int info = 0; |
|||
SPOTRF(&uplo, &n, a, &n, &info); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = 0; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int d_cholesky_factor(int n, double* a){ |
|||
char uplo = 'L'; |
|||
int info = 0; |
|||
DPOTRF(&uplo, &n, a, &n, &info); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = 0; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ |
|||
char uplo = 'L'; |
|||
int info = 0; |
|||
Complex8 zero; |
|||
zero.real = 0.0; |
|||
zero.real = 0.0; |
|||
CPOTRF(&uplo, &n, a, &n, &info); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = zero; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
|
|||
DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ |
|||
char uplo = 'L'; |
|||
int info = 0; |
|||
Complex16 zero; |
|||
zero.real = 0.0; |
|||
zero.real = 0.0; |
|||
ZPOTRF(&uplo, &n, a, &n, &info); |
|||
for (int i = 0; i < n; ++i) |
|||
{ |
|||
int index = i * n; |
|||
for (int j = 0; j < n && i > j; ++j) |
|||
{ |
|||
a[index + j] = zero; |
|||
} |
|||
} |
|||
return info; |
|||
} |
|||
} |
|||
@ -0,0 +1,14 @@ |
|||
#include "mkl_vml.h" |
|||
#include "common.h" |
|||
|
|||
DLLEXPORT void d_vector_add( const int n, const double x[], const double y[], double ret[]){ |
|||
vdAdd( n, x, y, ret ); |
|||
} |
|||
|
|||
DLLEXPORT void d_vector_subtract( const int n, const double x[], const double y[], double ret[]){ |
|||
vdSub( n, x, y, ret ); |
|||
} |
|||
|
|||
DLLEXPORT void d_vector_multiply( const int n, const double x[], const double y[], double ret[]){ |
|||
vdMul( n, x, y, ret ); |
|||
} |
|||
Loading…
Reference in new issue