/* -*- c -*- */
/*
*****************************************************************************
** INCLUDES **
*****************************************************************************
*/
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
#include "Python.h"
#include "numpy/arrayobject.h"
#include "numpy/ufuncobject.h"
#include "npy_pycompat.h"
#include "npy_config.h"
#include <stddef.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
static const char* umath_linalg_version_string = "0.1.5";
/*
****************************************************************************
* Debugging support *
****************************************************************************
*/
#define TRACE_TXT(...) do { fprintf (stderr, __VA_ARGS__); } while (0)
#define STACK_TRACE do {} while (0)
#define TRACE\
do { \
fprintf (stderr, \
"%s:%d:%s\n", \
__FILE__, \
__LINE__, \
__FUNCTION__); \
STACK_TRACE; \
} while (0)
#if 0
#include <execinfo.h>
void
dbg_stack_trace()
{
void *trace[32];
size_t size;
size = backtrace(trace, sizeof(trace)/sizeof(trace[0]));
backtrace_symbols_fd(trace, size, 1);
}
#undef STACK_TRACE
#define STACK_TRACE do { dbg_stack_trace(); } while (0)
#endif
/*
*****************************************************************************
* BLAS/LAPACK calling macros *
*****************************************************************************
*/
#ifdef NO_APPEND_FORTRAN
# define FNAME(x) x
#else
# define FNAME(x) x##_
#endif
typedef struct { float r, i; } f2c_complex;
typedef struct { double r, i; } f2c_doublecomplex;
/* typedef long int (*L_fp)(); */
extern int
FNAME(sgeev)(char *jobvl, char *jobvr, int *n,
float a[], int *lda, float wr[], float wi[],
float vl[], int *ldvl, float vr[], int *ldvr,
float work[], int lwork[],
int *info);
extern int
FNAME(dgeev)(char *jobvl, char *jobvr, int *n,
double a[], int *lda, double wr[], double wi[],
double vl[], int *ldvl, double vr[], int *ldvr,
double work[], int lwork[],
int *info);
extern int
FNAME(cgeev)(char *jobvl, char *jobvr, int *n,
f2c_doublecomplex a[], int *lda,
f2c_doublecomplex w[],
f2c_doublecomplex vl[], int *ldvl,
f2c_doublecomplex vr[], int *ldvr,
f2c_doublecomplex work[], int *lwork,
double rwork[],
int *info);
extern int
FNAME(zgeev)(char *jobvl, char *jobvr, int *n,
f2c_doublecomplex a[], int *lda,
f2c_doublecomplex w[],
f2c_doublecomplex vl[], int *ldvl,
f2c_doublecomplex vr[], int *ldvr,
f2c_doublecomplex work[], int *lwork,
double rwork[],
int *info);
extern int
FNAME(ssyevd)(char *jobz, char *uplo, int *n,
float a[], int *lda, float w[], float work[],
int *lwork, int iwork[], int *liwork,
int *info);
extern int
FNAME(dsyevd)(char *jobz, char *uplo, int *n,
double a[], int *lda, double w[], double work[],
int *lwork, int iwork[], int *liwork,
int *info);
extern int
FNAME(cheevd)(char *jobz, char *uplo, int *n,
f2c_complex a[], int *lda,
float w[], f2c_complex work[],
int *lwork, float rwork[], int *lrwork, int iwork[],
int *liwork,
int *info);
extern int
FNAME(zheevd)(char *jobz, char *uplo, int *n,
f2c_doublecomplex a[], int *lda,
double w[], f2c_doublecomplex work[],
int *lwork, double rwork[], int *lrwork, int iwork[],
int *liwork,
int *info);
extern int
FNAME(dgelsd)(int *m, int *n, int *nrhs,
double a[], int *lda, double b[], int *ldb,
double s[], double *rcond, int *rank,
double work[], int *lwork, int iwork[],
int *info);
extern int
FNAME(zgelsd)(int *m, int *n, int *nrhs,
f2c_doublecomplex a[], int *lda,
f2c_doublecomplex b[], int *ldb,
double s[], double *rcond, int *rank,
f2c_doublecomplex work[], int *lwork,
double rwork[], int iwork[],
int *info);
extern int
FNAME(sgesv)(int *n, int *nrhs,
float a[], int *lda,
int ipiv[],
float b[], int *ldb,
int *info);
extern int
FNAME(dgesv)(int *n, int *nrhs,
double a[], int *lda,
int ipiv[],
double b[], int *ldb,
int *info);
extern int
FNAME(cgesv)(int *n, int *nrhs,
f2c_complex a[], int *lda,
int ipiv[],
f2c_complex b[], int *ldb,
int *info);
extern int
FNAME(zgesv)(int *n, int *nrhs,
f2c_doublecomplex a[], int *lda,
int ipiv[],
f2c_doublecomplex b[], int *ldb,
int *info);
extern int
FNAME(sgetrf)(int *m, int *n,
float a[], int *lda,
int ipiv[],
int *info);
extern int
FNAME(dgetrf)(int *m, int *n,
double a[], int *lda,
int ipiv[],
int *info);
extern int
FNAME(cgetrf)(int *m, int *n,
f2c_complex a[], int *lda,
int ipiv[],
int *info);
extern int
FNAME(zgetrf)(int *m, int *n,
f2c_doublecomplex a[], int *lda,
int ipiv[],
int *info);
extern int
FNAME(spotrf)(char *uplo, int *n,
float a[], int *lda,
int *info);
extern int
FNAME(dpotrf)(char *uplo, int *n,
double a[], int *lda,
int *info);
extern int
FNAME(cpotrf)(char *uplo, int *n,
f2c_complex a[], int *lda,
int *info);
extern int
FNAME(zpotrf)(char *uplo, int *n,
f2c_doublecomplex a[], int *lda,
int *info);
extern int
FNAME(sgesdd)(char *jobz, int *m, int *n,
float a[], int *lda, float s[], float u[],
int *ldu, float vt[], int *ldvt, float work[],
int *lwork, int iwork[], int *info);
extern int
FNAME(dgesdd)(char *jobz, int *m, int *n,
double a[], int *lda, double s[], double u[],
int *ldu, double vt[], int *ldvt, double work[],
int *lwork, int iwork[], int *info);
extern int
FNAME(cgesdd)(char *jobz, int *m, int *n,
f2c_complex a[], int *lda,
float s[], f2c_complex u[], int *ldu,
f2c_complex vt[], int *ldvt,
f2c_complex work[], int *lwork,
float rwork[], int iwork[], int *info);
extern int
FNAME(zgesdd)(char *jobz, int *m, int *n,
f2c_doublecomplex a[], int *lda,
double s[], f2c_doublecomplex u[], int *ldu,
f2c_doublecomplex vt[], int *ldvt,
f2c_doublecomplex work[], int *lwork,
double rwork[], int iwork[], int *info);
extern int
FNAME(spotrs)(char *uplo, int *n, int *nrhs,
float a[], int *lda,
float b[], int *ldb,
int *info);
extern int
FNAME(dpotrs)(char *uplo, int *n, int *nrhs,
double a[], int *lda,
double b[], int *ldb,
int *info);
extern int
FNAME(cpotrs)(char *uplo, int *n, int *nrhs,
f2c_complex a[], int *lda,
f2c_complex b[], int *ldb,
int *info);
extern int
FNAME(zpotrs)(char *uplo, int *n, int *nrhs,
f2c_doublecomplex a[], int *lda,
f2c_doublecomplex b[], int *ldb,
int *info);
extern int
FNAME(spotri)(char *uplo, int *n,
float a[], int *lda,
int *info);
extern int
FNAME(dpotri)(char *uplo, int *n,
double a[], int *lda,
int *info);
extern int
FNAME(cpotri)(char *uplo, int *n,
f2c_complex a[], int *lda,
int *info);
extern int
FNAME(zpotri)(char *uplo, int *n,
f2c_doublecomplex a[], int *lda,
int *info);
extern int
FNAME(scopy)(int *n,
float *sx, int *incx,
float *sy, int *incy);
extern int
FNAME(dcopy)(int *n,
double *sx, int *incx,
double *sy, int *incy);
extern int
FNAME(ccopy)(int *n,
f2c_complex *sx, int *incx,
f2c_complex *sy, int *incy);
extern int
FNAME(zcopy)(int *n,
f2c_doublecomplex *sx, int *incx,
f2c_doublecomplex *sy, int *incy);
extern float
FNAME(sdot)(int *n,
float *sx, int *incx,
float *sy, int *incy);
extern double
FNAME(ddot)(int *n,
double *sx, int *incx,
double *sy, int *incy);
extern void
FNAME(cdotu)(f2c_complex *ret, int *n,
f2c_complex *sx, int *incx,
f2c_complex *sy, int *incy);
extern void
FNAME(zdotu)(f2c_doublecomplex *ret, int *n,
f2c_doublecomplex *sx, int *incx,
f2c_doublecomplex *sy, int *incy);
extern void
FNAME(cdotc)(f2c_complex *ret, int *n,
f2c_complex *sx, int *incx,
f2c_complex *sy, int *incy);
extern void
FNAME(zdotc)(f2c_doublecomplex *ret, int *n,
f2c_doublecomplex *sx, int *incx,
f2c_doublecomplex *sy, int *incy);
extern int
FNAME(sgemm)(char *transa, char *transb,
int *m, int *n, int *k,
float *alpha,
float *a, int *lda,
float *b, int *ldb,
float *beta,
float *c, int *ldc);
extern int
FNAME(dgemm)(char *transa, char *transb,
int *m, int *n, int *k,
double *alpha,
double *a, int *lda,
double *b, int *ldb,
double *beta,
double *c, int *ldc);
extern int
FNAME(cgemm)(char *transa, char *transb,
int *m, int *n, int *k,
f2c_complex *alpha,
f2c_complex *a, int *lda,
f2c_complex *b, int *ldb,
f2c_complex *beta,
f2c_complex *c, int *ldc);
extern int
FNAME(zgemm)(char *transa, char *transb,
int *m, int *n, int *k,
f2c_doublecomplex *alpha,
f2c_doublecomplex *a, int *lda,
f2c_doublecomplex *b, int *ldb,
f2c_doublecomplex *beta,
f2c_doublecomplex *c, int *ldc);
#define LAPACK_T(FUNC) \
TRACE_TXT("Calling LAPACK ( " # FUNC " )\n"); \
FNAME(FUNC)
#define BLAS(FUNC) \
FNAME(FUNC)
#define LAPACK(FUNC) \
FNAME(FUNC)
typedef int fortran_int;
typedef float fortran_real;
typedef double fortran_doublereal;
typedef f2c_complex fortran_complex;
typedef f2c_doublecomplex fortran_doublecomplex;
/*
*****************************************************************************
** Some handy functions **
*****************************************************************************
*/
static NPY_INLINE void *
offset_ptr(void* ptr, ptrdiff_t offset)
{
return (void*)((npy_uint8*)ptr + offset);
}
static NPY_INLINE int
get_fp_invalid_and_clear(void)
{
int status;
status = npy_clear_floatstatus();
return !!(status & NPY_FPE_INVALID);
}
static NPY_INLINE void
set_fp_invalid_or_clear(int error_occurred)
{
if (error_occurred) {
npy_set_floatstatus_invalid();
}
else {
npy_clear_floatstatus();
}
}
/*
*****************************************************************************
** Some handy constants **
*****************************************************************************
*/
#define UMATH_LINALG_MODULE_NAME "_umath_linalg"
typedef union {
fortran_complex f;
npy_cfloat npy;
float array[2];
} COMPLEX_t;
typedef union {
fortran_doublecomplex f;
npy_cdouble npy;
double array[2];
} DOUBLECOMPLEX_t;
static float s_one;
static float s_zero;
static float s_minus_one;
static float s_ninf;
static float s_nan;
static double d_one;
static double d_zero;
static double d_minus_one;
static double d_ninf;
static double d_nan;
static COMPLEX_t c_one;
static COMPLEX_t c_zero;
static COMPLEX_t c_minus_one;
static COMPLEX_t c_ninf;
static COMPLEX_t c_nan;
static DOUBLECOMPLEX_t z_one;
static DOUBLECOMPLEX_t z_zero;
static DOUBLECOMPLEX_t z_minus_one;
static DOUBLECOMPLEX_t z_ninf;
static DOUBLECOMPLEX_t z_nan;
static void init_constants(void)
{
/*
this is needed as NPY_INFINITY and NPY_NAN macros
can't be used as initializers. I prefer to just set
all the constants the same way.
*/
s_one = 1.0f;
s_zero = 0.0f;
s_minus_one = -1.0f;
s_ninf = -NPY_INFINITYF;
s_nan = NPY_NANF;
d_one = 1.0;
d_zero = 0.0;
d_minus_one = -1.0;
d_ninf = -NPY_INFINITY;
d_nan = NPY_NAN;
c_one.array[0] = 1.0f;
c_one.array[1] = 0.0f;
c_zero.array[0] = 0.0f;
c_zero.array[1] = 0.0f;
c_minus_one.array[0] = -1.0f;
c_minus_one.array[1] = 0.0f;
c_ninf.array[0] = -NPY_INFINITYF;
c_ninf.array[1] = 0.0f;
c_nan.array[0] = NPY_NANF;
c_nan.array[1] = NPY_NANF;
z_one.array[0] = 1.0;
z_one.array[1] = 0.0;
z_zero.array[0] = 0.0;
z_zero.array[1] = 0.0;
z_minus_one.array[0] = -1.0;
z_minus_one.array[1] = 0.0;
z_ninf.array[0] = -NPY_INFINITY;
z_ninf.array[1] = 0.0;
z_nan.array[0] = NPY_NAN;
z_nan.array[1] = NPY_NAN;
}
/*
*****************************************************************************
** Structs used for data rearrangement **
*****************************************************************************
*/
/* this struct contains information about how to linearize in a local buffer
a matrix so that it can be used by blas functions.
All strides are specified in number of elements (similar to what blas
expects)
dst_row_strides: number of elements between different row. Matrix is
considered row-major
dst_column_strides: number of elements between different columns in the
destination buffer
rows: number of rows of the matrix
columns: number of columns of the matrix
src_row_strides: strides needed to access the next row in the source matrix
src_column_strides: strides needed to access the next column in the source
matrix
*/
typedef struct linearize_data_struct
{
size_t rows;
size_t columns;
ptrdiff_t row_strides;
ptrdiff_t column_strides;
} LINEARIZE_DATA_t;
static NPY_INLINE void
init_linearize_data(LINEARIZE_DATA_t *lin_data,
int rows,
int columns,
ptrdiff_t row_strides,
ptrdiff_t column_strides)
{
lin_data->rows = rows;
lin_data->columns = columns;
lin_data->row_strides = row_strides;
lin_data->column_strides = column_strides;
}
static NPY_INLINE void
dump_ufunc_object(PyUFuncObject* ufunc)
{
TRACE_TXT("\n\n%s '%s' (%d input(s), %d output(s), %d specialization(s).\n",
ufunc->core_enabled? "generalized ufunc" : "scalar ufunc",
ufunc->name, ufunc->nin, ufunc->nout, ufunc->ntypes);
if (ufunc->core_enabled) {
int arg;
int dim;
TRACE_TXT("\t%s (%d dimension(s) detected).\n",
ufunc->core_signature, ufunc->core_num_dim_ix);
for (arg = 0; arg < ufunc->nargs; arg++){
int * arg_dim_ix = ufunc->core_dim_ixs + ufunc->core_offsets[arg];
TRACE_TXT("\t\targ %d (%s) has %d dimension(s): (",
arg, arg < ufunc->nin? "INPUT" : "OUTPUT",
ufunc->core_num_dims[arg]);
for (dim = 0; dim < ufunc->core_num_dims[arg]; dim ++) {
TRACE_TXT(" %d", arg_dim_ix[dim]);
}
TRACE_TXT(" )\n");
}
}
}
static NPY_INLINE void
dump_linearize_data(const char* name, const LINEARIZE_DATA_t* params)
{
TRACE_TXT("\n\t%s rows: %zd columns: %zd"\
"\n\t\trow_strides: %td column_strides: %td"\
"\n", name, params->rows, params->columns,
params->row_strides, params->column_strides);
}
static NPY_INLINE float
FLOAT_add(float op1, float op2)
{
return op1 + op2;
}
static NPY_INLINE double
DOUBLE_add(double op1, double op2)
{
return op1 + op2;
}
static NPY_INLINE COMPLEX_t
CFLOAT_add(COMPLEX_t op1, COMPLEX_t op2)
{
COMPLEX_t result;
result.array[0] = op1.array[0] + op2.array[0];
result.array[1] = op1.array[1] + op2.array[1];
return result;
}
static NPY_INLINE DOUBLECOMPLEX_t
CDOUBLE_add(DOUBLECOMPLEX_t op1, DOUBLECOMPLEX_t op2)
{
DOUBLECOMPLEX_t result;
result.array[0] = op1.array[0] + op2.array[0];
result.array[1] = op1.array[1] + op2.array[1];
return result;
}
static NPY_INLINE float
FLOAT_mul(float op1, float op2)
{
return op1*op2;
}
static NPY_INLINE double
DOUBLE_mul(double op1, double op2)
{
return op1*op2;
}
static NPY_INLINE COMPLEX_t
CFLOAT_mul(COMPLEX_t op1, COMPLEX_t op2)
{
COMPLEX_t result;
result.array[0] = op1.array[0]*op2.array[0] - op1.array[1]*op2.array[1];
result.array[1] = op1.array[1]*op2.array[0] + op1.array[0]*op2.array[1];
return result;
}
static NPY_INLINE DOUBLECOMPLEX_t
CDOUBLE_mul(DOUBLECOMPLEX_t op1, DOUBLECOMPLEX_t op2)
{
DOUBLECOMPLEX_t result;
result.array[0] = op1.array[0]*op2.array[0] - op1.array[1]*op2.array[1];
result.array[1] = op1.array[1]*op2.array[0] + op1.array[0]*op2.array[1];
return result;
}
static NPY_INLINE float
FLOAT_mulc(float op1, float op2)
{
return op1*op2;
}
static NPY_INLINE double
DOUBLE_mulc(float op1, float op2)
{
return op1*op2;
}
static NPY_INLINE COMPLEX_t
CFLOAT_mulc(COMPLEX_t op1, COMPLEX_t op2)
{
COMPLEX_t result;
result.array[0] = op1.array[0]*op2.array[0] + op1.array[1]*op2.array[1];
result.array[1] = op1.array[0]*op2.array[1] - op1.array[1]*op2.array[0];
return result;
}
static NPY_INLINE DOUBLECOMPLEX_t
CDOUBLE_mulc(DOUBLECOMPLEX_t op1, DOUBLECOMPLEX_t op2)
{
DOUBLECOMPLEX_t result;
result.array[0] = op1.array[0]*op2.array[0] + op1.array[1]*op2.array[1];
result.array[1] = op1.array[0]*op2.array[1] - op1.array[1]*op2.array[0];
return result;
}
static NPY_INLINE void
print_FLOAT(npy_float s)
{
TRACE_TXT(" %8.4f", s);
}
static NPY_INLINE void
print_DOUBLE(npy_double d)
{
TRACE_TXT(" %10.6f", d);
}
static NPY_INLINE void
print_CFLOAT(npy_cfloat c)
{
float* c_parts = (float*)&c;
TRACE_TXT("(%8.4f, %8.4fj)", c_parts[0], c_parts[1]);
}
static NPY_INLINE void
print_CDOUBLE(npy_cdouble z)
{
double* z_parts = (double*)&z;
TRACE_TXT("(%8.4f, %8.4fj)", z_parts[0], z_parts[1]);
}
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = npy_float, npy_double, npy_cfloat, npy_cdouble#
*/
static NPY_INLINE void
dump_@TYPE@_matrix(const char* name,
size_t rows, size_t columns,
const @typ@* ptr)
{
size_t i, j;
TRACE_TXT("\n%s %p (%zd, %zd)\n", name, ptr, rows, columns);
for (i = 0; i < rows; i++)
{
TRACE_TXT("| ");
for (j = 0; j < columns; j++)
{
print_@TYPE@(ptr[j*rows + i]);
TRACE_TXT(", ");
}
TRACE_TXT(" |\n");
}
}
/**end repeat**/
/*
*****************************************************************************
** Basics **
*****************************************************************************
*/
static NPY_INLINE fortran_int
fortran_int_min(fortran_int x, fortran_int y) {
return x < y ? x : y;
}
static NPY_INLINE fortran_int
fortran_int_max(fortran_int x, fortran_int y) {
return x > y ? x : y;
}
#define INIT_OUTER_LOOP_1 \
npy_intp dN = *dimensions++;\
npy_intp N_;\
npy_intp s0 = *steps++;
#define INIT_OUTER_LOOP_2 \
INIT_OUTER_LOOP_1\
npy_intp s1 = *steps++;
#define INIT_OUTER_LOOP_3 \
INIT_OUTER_LOOP_2\
npy_intp s2 = *steps++;
#define INIT_OUTER_LOOP_4 \
INIT_OUTER_LOOP_3\
npy_intp s3 = *steps++;
#define INIT_OUTER_LOOP_5 \
INIT_OUTER_LOOP_4\
npy_intp s4 = *steps++;
#define INIT_OUTER_LOOP_6 \
INIT_OUTER_LOOP_5\
npy_intp s5 = *steps++;
#define BEGIN_OUTER_LOOP_2 \
for (N_ = 0;\
N_ < dN;\
N_++, args[0] += s0,\
args[1] += s1) {
#define BEGIN_OUTER_LOOP_3 \
for (N_ = 0;\
N_ < dN;\
N_++, args[0] += s0,\
args[1] += s1,\
args[2] += s2) {
#define BEGIN_OUTER_LOOP_4 \
for (N_ = 0;\
N_ < dN;\
N_++, args[0] += s0,\
args[1] += s1,\
args[2] += s2,\
args[3] += s3) {
#define BEGIN_OUTER_LOOP_5 \
for (N_ = 0;\
N_ < dN;\
N_++, args[0] += s0,\
args[1] += s1,\
args[2] += s2,\
args[3] += s3,\
args[4] += s4) {
#define BEGIN_OUTER_LOOP_6 \
for (N_ = 0;\
N_ < dN;\
N_++, args[0] += s0,\
args[1] += s1,\
args[2] += s2,\
args[3] += s3,\
args[4] += s4,\
args[5] += s5) {
#define END_OUTER_LOOP }
static NPY_INLINE void
update_pointers(npy_uint8** bases, ptrdiff_t* offsets, size_t count)
{
size_t i;
for (i = 0; i < count; ++i) {
bases[i] += offsets[i];
}
}
/* disable -Wmaybe-uninitialized as there is some code that generate false
positives with this warning
*/
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
/*
*****************************************************************************
** HELPER FUNCS **
*****************************************************************************
*/
/* rearranging of 2D matrices using blas */
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = float, double, COMPLEX_t, DOUBLECOMPLEX_t#
#copy = scopy, dcopy, ccopy, zcopy#
#nan = s_nan, d_nan, c_nan, z_nan#
*/
static NPY_INLINE void *
linearize_@TYPE@_matrix(void *dst_in,
void *src_in,
const LINEARIZE_DATA_t* data)
{
@typ@ *src = (@typ@ *) src_in;
@typ@ *dst = (@typ@ *) dst_in;
if (dst) {
int i, j;
@typ@* rv = dst;
fortran_int columns = (fortran_int)data->columns;
fortran_int column_strides =
(fortran_int)(data->column_strides/sizeof(@typ@));
fortran_int one = 1;
for (i = 0; i < data->rows; i++) {
if (column_strides > 0) {
FNAME(@copy@)(&columns,
(void*)src, &column_strides,
(void*)dst, &one);
}
else if (column_strides < 0) {
FNAME(@copy@)(&columns,
(void*)((@typ@*)src + (columns-1)*column_strides),
&column_strides,
(void*)dst, &one);
}
else {
/*
* Zero stride has undefined behavior in some BLAS
* implementations (e.g. OSX Accelerate), so do it
* manually
*/
for (j = 0; j < columns; ++j) {
memcpy((@typ@*)dst + j, (@typ@*)src, sizeof(@typ@));
}
}
src += data->row_strides/sizeof(@typ@);
dst += data->columns;
}
return rv;
} else {
return src;
}
}
static NPY_INLINE void *
delinearize_@TYPE@_matrix(void *dst_in,
void *src_in,
const LINEARIZE_DATA_t* data)
{
@typ@ *src = (@typ@ *) src_in;
@typ@ *dst = (@typ@ *) dst_in;
if (src) {
int i;
@typ@ *rv = src;
fortran_int columns = (fortran_int)data->columns;
fortran_int column_strides =
(fortran_int)(data->column_strides/sizeof(@typ@));
fortran_int one = 1;
for (i = 0; i < data->rows; i++) {
if (column_strides > 0) {
FNAME(@copy@)(&columns,
(void*)src, &one,
(void*)dst, &column_strides);
}
else if (column_strides < 0) {
FNAME(@copy@)(&columns,
(void*)src, &one,
(void*)((@typ@*)dst + (columns-1)*column_strides),
&column_strides);
}
else {
/*
* Zero stride has undefined behavior in some BLAS
* implementations (e.g. OSX Accelerate), so do it
* manually
*/
if (columns > 0) {
memcpy((@typ@*)dst,
(@typ@*)src + (columns-1),
sizeof(@typ@));
}
}
src += data->columns;
dst += data->row_strides/sizeof(@typ@);
}
return rv;
} else {
return src;
}
}
static NPY_INLINE void
nan_@TYPE@_matrix(void *dst_in, const LINEARIZE_DATA_t* data)
{
@typ@ *dst = (@typ@ *) dst_in;
int i, j;
for (i = 0; i < data->rows; i++) {
@typ@ *cp = dst;
ptrdiff_t cs = data->column_strides/sizeof(@typ@);
for (j = 0; j < data->columns; ++j) {
*cp = @nan@;
cp += cs;
}
dst += data->row_strides/sizeof(@typ@);
}
}
/**end repeat**/
/* identity square matrix generation */
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = float, double, COMPLEX_t, DOUBLECOMPLEX_t#
#cblas_type = s, d, c, z#
*/
static NPY_INLINE void
identity_@TYPE@_matrix(void *ptr, size_t n)
{
size_t i;
@typ@ *matrix = (@typ@*) ptr;
/* in IEEE floating point, zeroes are represented as bitwise 0 */
memset(matrix, 0, n*n*sizeof(@typ@));
for (i = 0; i < n; ++i)
{
*matrix = @cblas_type@_one;
matrix += n+1;
}
}
/**end repeat**/
/* lower/upper triangular matrix using blas (in place) */
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = float, double, COMPLEX_t, DOUBLECOMPLEX_t#
#cblas_type = s, d, c, z#
*/
static NPY_INLINE void
triu_@TYPE@_matrix(void *ptr, size_t n)
{
size_t i, j;
@typ@ *matrix = (@typ@*)ptr;
matrix += n;
for (i = 1; i < n; ++i) {
for (j = 0; j < i; ++j) {
matrix[j] = @cblas_type@_zero;
}
matrix += n;
}
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* Determinants */
/**begin repeat
#TYPE = FLOAT, DOUBLE#
#typ = npy_float, npy_double#
#log_func = npy_logf, npy_log#
#exp_func = npy_expf, npy_exp#
#zero = 0.0f, 0.0#
*/
static NPY_INLINE void
@TYPE@_slogdet_from_factored_diagonal(@typ@* src,
fortran_int m,
@typ@ *sign,
@typ@ *logdet)
{
@typ@ acc_sign = *sign;
@typ@ acc_logdet = @zero@;
int i;
for (i = 0; i < m; i++) {
@typ@ abs_element = *src;
if (abs_element < @zero@) {
acc_sign = -acc_sign;
abs_element = -abs_element;
}
acc_logdet += @log_func@(abs_element);
src += m+1;
}
*sign = acc_sign;
*logdet = acc_logdet;
}
static NPY_INLINE @typ@
@TYPE@_det_from_slogdet(@typ@ sign, @typ@ logdet)
{
@typ@ result = sign * @exp_func@(logdet);
return result;
}
/**end repeat**/
/**begin repeat
#TYPE = CFLOAT, CDOUBLE#
#typ = npy_cfloat, npy_cdouble#
#basetyp = npy_float, npy_double#
#abs_func = npy_cabsf, npy_cabs#
#log_func = npy_logf, npy_log#
#exp_func = npy_expf, npy_exp#
#zero = 0.0f, 0.0#
*/
#define RE(COMPLEX) (((@basetyp@*)(&COMPLEX))[0])
#define IM(COMPLEX) (((@basetyp@*)(&COMPLEX))[1])
static NPY_INLINE @typ@
@TYPE@_mult(@typ@ op1, @typ@ op2)
{
@typ@ rv;
RE(rv) = RE(op1)*RE(op2) - IM(op1)*IM(op2);
IM(rv) = RE(op1)*IM(op2) + IM(op1)*RE(op2);
return rv;
}
static NPY_INLINE void
@TYPE@_slogdet_from_factored_diagonal(@typ@* src,
fortran_int m,
@typ@ *sign,
@basetyp@ *logdet)
{
int i;
@typ@ sign_acc = *sign;
@basetyp@ logdet_acc = @zero@;
for (i = 0; i < m; i++)
{
@basetyp@ abs_element = @abs_func@(*src);
@typ@ sign_element;
RE(sign_element) = RE(*src) / abs_element;
IM(sign_element) = IM(*src) / abs_element;
sign_acc = @TYPE@_mult(sign_acc, sign_element);
logdet_acc += @log_func@(abs_element);
src += m + 1;
}
*sign = sign_acc;
*logdet = logdet_acc;
}
static NPY_INLINE @typ@
@TYPE@_det_from_slogdet(@typ@ sign, @basetyp@ logdet)
{
@typ@ tmp;
RE(tmp) = @exp_func@(logdet);
IM(tmp) = @zero@;
return @TYPE@_mult(sign, tmp);
}
#undef RE
#undef IM
/**end repeat**/
/* As in the linalg package, the determinant is computed via LU factorization
* using LAPACK.
* slogdet computes sign + log(determinant).
* det computes sign * exp(slogdet).
*/
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = npy_float, npy_double, npy_cfloat, npy_cdouble#
#basetyp = npy_float, npy_double, npy_float, npy_double#
#cblas_type = s, d, c, z#
*/
static NPY_INLINE void
@TYPE@_slogdet_single_element(fortran_int m,
void* src,
fortran_int* pivots,
@typ@ *sign,
@basetyp@ *logdet)
{
fortran_int info = 0;
fortran_int lda = fortran_int_max(m, 1);
int i;
/* note: done in place */
LAPACK(@cblas_type@getrf)(&m, &m, (void *)src, &lda, pivots, &info);
if (info == 0) {
int change_sign = 0;
/* note: fortran uses 1 based indexing */
for (i = 0; i < m; i++)
{
change_sign += (pivots[i] != (i+1));
}
memcpy(sign,
(change_sign % 2)?
&@cblas_type@_minus_one :
&@cblas_type@_one
, sizeof(*sign));
@TYPE@_slogdet_from_factored_diagonal(src, m, sign, logdet);
} else {
/*
if getrf fails, use 0 as sign and -inf as logdet
*/
memcpy(sign, &@cblas_type@_zero, sizeof(*sign));
memcpy(logdet, &@cblas_type@_ninf, sizeof(*logdet));
}
}
static void
@TYPE@_slogdet(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
fortran_int m;
npy_uint8 *tmp_buff = NULL;
size_t matrix_size;
size_t pivot_size;
size_t safe_m;
/* notes:
* matrix will need to be copied always, as factorization in lapack is
* made inplace
* matrix will need to be in column-major order, as expected by lapack
* code (fortran)
* always a square matrix
* need to allocate memory for both, matrix_buffer and pivot buffer
*/
INIT_OUTER_LOOP_3
m = (fortran_int) dimensions[0];
safe_m = m;
matrix_size = safe_m * safe_m * sizeof(@typ@);
pivot_size = safe_m * sizeof(fortran_int);
tmp_buff = (npy_uint8 *)malloc(matrix_size + pivot_size);
if (tmp_buff) {
LINEARIZE_DATA_t lin_data;
/* swapped steps to get matrix in FORTRAN order */
init_linearize_data(&lin_data, m, m,
(ptrdiff_t)steps[1],
(ptrdiff_t)steps[0]);
BEGIN_OUTER_LOOP_3
linearize_@TYPE@_matrix(tmp_buff, args[0], &lin_data);
@TYPE@_slogdet_single_element(m,
(void*)tmp_buff,
(fortran_int*)(tmp_buff+matrix_size),
(@typ@*)args[1],
(@basetyp@*)args[2]);
END_OUTER_LOOP
free(tmp_buff);
}
}
static void
@TYPE@_det(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
fortran_int m;
npy_uint8 *tmp_buff;
size_t matrix_size;
size_t pivot_size;
size_t safe_m;
/* notes:
* matrix will need to be copied always, as factorization in lapack is
* made inplace
* matrix will need to be in column-major order, as expected by lapack
* code (fortran)
* always a square matrix
* need to allocate memory for both, matrix_buffer and pivot buffer
*/
INIT_OUTER_LOOP_2
m = (fortran_int) dimensions[0];
safe_m = m;
matrix_size = safe_m * safe_m * sizeof(@typ@);
pivot_size = safe_m * sizeof(fortran_int);
tmp_buff = (npy_uint8 *)malloc(matrix_size + pivot_size);
if (tmp_buff) {
LINEARIZE_DATA_t lin_data;
@typ@ sign;
@basetyp@ logdet;
/* swapped steps to get matrix in FORTRAN order */
init_linearize_data(&lin_data, m, m,
(ptrdiff_t)steps[1],
(ptrdiff_t)steps[0]);
BEGIN_OUTER_LOOP_2
linearize_@TYPE@_matrix(tmp_buff, args[0], &lin_data);
@TYPE@_slogdet_single_element(m,
(void*)tmp_buff,
(fortran_int*)(tmp_buff+matrix_size),
&sign,
&logdet);
*(@typ@ *)args[1] = @TYPE@_det_from_slogdet(sign, logdet);
END_OUTER_LOOP
free(tmp_buff);
}
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* Eigh family */
typedef struct eigh_params_struct {
void *A; /* matrix */
void *W; /* eigenvalue vector */
void *WORK; /* main work buffer */
void *RWORK; /* secondary work buffer (for complex versions) */
void *IWORK;
fortran_int N;
fortran_int LWORK;
fortran_int LRWORK;
fortran_int LIWORK;
char JOBZ;
char UPLO;
fortran_int LDA;
} EIGH_PARAMS_t;
/**begin repeat
#TYPE = FLOAT, DOUBLE#
#typ = npy_float, npy_double#
#ftyp = fortran_real, fortran_doublereal#
#lapack_func = ssyevd, dsyevd#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(EIGH_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBZ, ¶ms->UPLO, ¶ms->N,
params->A, ¶ms->LDA, params->W,
params->WORK, ¶ms->LWORK,
params->IWORK, ¶ms->LIWORK,
&rv);
return rv;
}
/*
* Initialize the parameters to use in for the lapack function _syevd
* Handles buffer allocation
*/
static NPY_INLINE int
init_@lapack_func@(EIGH_PARAMS_t* params, char JOBZ, char UPLO,
fortran_int N)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *mem_buff2 = NULL;
fortran_int lwork;
fortran_int liwork;
npy_uint8 *a, *w, *work, *iwork;
size_t safe_N = N;
size_t alloc_size = safe_N * (safe_N + 1) * sizeof(@typ@);
fortran_int lda = fortran_int_max(N, 1);
mem_buff = malloc(alloc_size);
if (!mem_buff) {
goto error;
}
a = mem_buff;
w = mem_buff + safe_N * safe_N * sizeof(@typ@);
params->A = a;
params->W = w;
params->RWORK = NULL; /* unused */
params->N = N;
params->LRWORK = 0; /* unused */
params->JOBZ = JOBZ;
params->UPLO = UPLO;
params->LDA = lda;
/* Work size query */
{
@typ@ query_work_size;
fortran_int query_iwork_size;
params->LWORK = -1;
params->LIWORK = -1;
params->WORK = &query_work_size;
params->IWORK = &query_iwork_size;
if (call_@lapack_func@(params) != 0) {
goto error;
}
lwork = (fortran_int)query_work_size;
liwork = query_iwork_size;
}
mem_buff2 = malloc(lwork*sizeof(@typ@) + liwork*sizeof(fortran_int));
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
iwork = mem_buff2 + lwork*sizeof(@typ@);
params->LWORK = lwork;
params->WORK = work;
params->LIWORK = liwork;
params->IWORK = iwork;
return 1;
error:
/* something failed */
memset(params, 0, sizeof(*params));
free(mem_buff2);
free(mem_buff);
return 0;
}
/**end repeat**/
/**begin repeat
#TYPE = CFLOAT, CDOUBLE#
#typ = npy_cfloat, npy_cdouble#
#basetyp = npy_float, npy_double#
#ftyp = fortran_complex, fortran_doublecomplex#
#fbasetyp = fortran_real, fortran_doublereal#
#lapack_func = cheevd, zheevd#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(EIGH_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBZ, ¶ms->UPLO, ¶ms->N,
params->A, ¶ms->LDA, params->W,
params->WORK, ¶ms->LWORK,
params->RWORK, ¶ms->LRWORK,
params->IWORK, ¶ms->LIWORK,
&rv);
return rv;
}
/*
* Initialize the parameters to use in for the lapack function _heev
* Handles buffer allocation
*/
static NPY_INLINE int
init_@lapack_func@(EIGH_PARAMS_t *params,
char JOBZ,
char UPLO,
fortran_int N)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *mem_buff2 = NULL;
fortran_int lwork;
fortran_int lrwork;
fortran_int liwork;
npy_uint8 *a, *w, *work, *rwork, *iwork;
size_t safe_N = N;
fortran_int lda = fortran_int_max(N, 1);
mem_buff = malloc(safe_N * safe_N * sizeof(@typ@) +
safe_N * sizeof(@basetyp@));
if (!mem_buff) {
goto error;
}
a = mem_buff;
w = mem_buff + safe_N * safe_N * sizeof(@typ@);
params->A = a;
params->W = w;
params->N = N;
params->JOBZ = JOBZ;
params->UPLO = UPLO;
params->LDA = lda;
/* Work size query */
{
@ftyp@ query_work_size;
@fbasetyp@ query_rwork_size;
fortran_int query_iwork_size;
params->LWORK = -1;
params->LRWORK = -1;
params->LIWORK = -1;
params->WORK = &query_work_size;
params->RWORK = &query_rwork_size;
params->IWORK = &query_iwork_size;
if (call_@lapack_func@(params) != 0) {
goto error;
}
lwork = (fortran_int)*(@fbasetyp@*)&query_work_size;
lrwork = (fortran_int)query_rwork_size;
liwork = query_iwork_size;
}
mem_buff2 = malloc(lwork*sizeof(@typ@) +
lrwork*sizeof(@basetyp@) +
liwork*sizeof(fortran_int));
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
rwork = work + lwork*sizeof(@typ@);
iwork = rwork + lrwork*sizeof(@basetyp@);
params->WORK = work;
params->RWORK = rwork;
params->IWORK = iwork;
params->LWORK = lwork;
params->LRWORK = lrwork;
params->LIWORK = liwork;
return 1;
/* something failed */
error:
memset(params, 0, sizeof(*params));
free(mem_buff2);
free(mem_buff);
return 0;
}
/**end repeat**/
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#BASETYPE = FLOAT, DOUBLE, FLOAT, DOUBLE#
#typ = npy_float, npy_double, npy_cfloat, npy_cdouble#
#basetyp = npy_float, npy_double, npy_float, npy_double#
#lapack_func = ssyevd, dsyevd, cheevd, zheevd#
**/
/*
* (M, M)->(M,)(M, M)
* dimensions[1] -> M
* args[0] -> A[in]
* args[1] -> W
* args[2] -> A[out]
*/
static NPY_INLINE void
release_@lapack_func@(EIGH_PARAMS_t *params)
{
/* allocated memory in A and WORK */
free(params->A);
free(params->WORK);
memset(params, 0, sizeof(*params));
}
static NPY_INLINE void
@TYPE@_eigh_wrapper(char JOBZ,
char UPLO,
char**args,
npy_intp* dimensions,
npy_intp* steps)
{
ptrdiff_t outer_steps[3];
size_t iter;
size_t outer_dim = *dimensions++;
size_t op_count = (JOBZ=='N')?2:3;
EIGH_PARAMS_t eigh_params;
int error_occurred = get_fp_invalid_and_clear();
for (iter = 0; iter < op_count; ++iter) {
outer_steps[iter] = (ptrdiff_t) steps[iter];
}
steps += op_count;
if (init_@lapack_func@(&eigh_params,
JOBZ,
UPLO,
(fortran_int)dimensions[0])) {
LINEARIZE_DATA_t matrix_in_ld;
LINEARIZE_DATA_t eigenvectors_out_ld;
LINEARIZE_DATA_t eigenvalues_out_ld;
init_linearize_data(&matrix_in_ld,
eigh_params.N, eigh_params.N,
steps[1], steps[0]);
init_linearize_data(&eigenvalues_out_ld,
1, eigh_params.N,
0, steps[2]);
if ('V' == eigh_params.JOBZ) {
init_linearize_data(&eigenvectors_out_ld,
eigh_params.N, eigh_params.N,
steps[4], steps[3]);
}
for (iter = 0; iter < outer_dim; ++iter) {
int not_ok;
/* copy the matrix in */
linearize_@TYPE@_matrix(eigh_params.A, args[0], &matrix_in_ld);
not_ok = call_@lapack_func@(&eigh_params);
if (!not_ok) {
/* lapack ok, copy result out */
delinearize_@BASETYPE@_matrix(args[1],
eigh_params.W,
&eigenvalues_out_ld);
if ('V' == eigh_params.JOBZ) {
delinearize_@TYPE@_matrix(args[2],
eigh_params.A,
&eigenvectors_out_ld);
}
} else {
/* lapack fail, set result to nan */
error_occurred = 1;
nan_@BASETYPE@_matrix(args[1], &eigenvalues_out_ld);
if ('V' == eigh_params.JOBZ) {
nan_@TYPE@_matrix(args[2], &eigenvectors_out_ld);
}
}
update_pointers((npy_uint8**)args, outer_steps, op_count);
}
release_@lapack_func@(&eigh_params);
}
set_fp_invalid_or_clear(error_occurred);
}
/**end repeat**/
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
*/
static void
@TYPE@_eighlo(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_eigh_wrapper('V', 'L', args, dimensions, steps);
}
static void
@TYPE@_eighup(char **args,
npy_intp *dimensions,
npy_intp *steps,
void* NPY_UNUSED(func))
{
@TYPE@_eigh_wrapper('V', 'U', args, dimensions, steps);
}
static void
@TYPE@_eigvalshlo(char **args,
npy_intp *dimensions,
npy_intp *steps,
void* NPY_UNUSED(func))
{
@TYPE@_eigh_wrapper('N', 'L', args, dimensions, steps);
}
static void
@TYPE@_eigvalshup(char **args,
npy_intp *dimensions,
npy_intp *steps,
void* NPY_UNUSED(func))
{
@TYPE@_eigh_wrapper('N', 'U', args, dimensions, steps);
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* Solve family (includes inv) */
typedef struct gesv_params_struct
{
void *A; /* A is (N, N) of base type */
void *B; /* B is (N, NRHS) of base type */
fortran_int * IPIV; /* IPIV is (N) */
fortran_int N;
fortran_int NRHS;
fortran_int LDA;
fortran_int LDB;
} GESV_PARAMS_t;
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#typ = npy_float, npy_double, npy_cfloat, npy_cdouble#
#ftyp = fortran_real, fortran_doublereal,
fortran_complex, fortran_doublecomplex#
#lapack_func = sgesv, dgesv, cgesv, zgesv#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(GESV_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->N, ¶ms->NRHS,
params->A, ¶ms->LDA,
params->IPIV,
params->B, ¶ms->LDB,
&rv);
return rv;
}
/*
* Initialize the parameters to use in for the lapack function _heev
* Handles buffer allocation
*/
static NPY_INLINE int
init_@lapack_func@(GESV_PARAMS_t *params, fortran_int N, fortran_int NRHS)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *a, *b, *ipiv;
size_t safe_N = N;
size_t safe_NRHS = NRHS;
fortran_int ld = fortran_int_max(N, 1);
mem_buff = malloc(safe_N * safe_N * sizeof(@ftyp@) +
safe_N * safe_NRHS*sizeof(@ftyp@) +
safe_N * sizeof(fortran_int));
if (!mem_buff) {
goto error;
}
a = mem_buff;
b = a + safe_N * safe_N * sizeof(@ftyp@);
ipiv = b + safe_N * safe_NRHS * sizeof(@ftyp@);
params->A = a;
params->B = b;
params->IPIV = (fortran_int*)ipiv;
params->N = N;
params->NRHS = NRHS;
params->LDA = ld;
params->LDB = ld;
return 1;
error:
free(mem_buff);
memset(params, 0, sizeof(*params));
return 0;
}
static NPY_INLINE void
release_@lapack_func@(GESV_PARAMS_t *params)
{
/* memory block base is in A */
free(params->A);
memset(params, 0, sizeof(*params));
}
static void
@TYPE@_solve(char **args, npy_intp *dimensions, npy_intp *steps,
void *NPY_UNUSED(func))
{
GESV_PARAMS_t params;
fortran_int n, nrhs;
int error_occurred = get_fp_invalid_and_clear();
INIT_OUTER_LOOP_3
n = (fortran_int)dimensions[0];
nrhs = (fortran_int)dimensions[1];
if (init_@lapack_func@(¶ms, n, nrhs)) {
LINEARIZE_DATA_t a_in, b_in, r_out;
init_linearize_data(&a_in, n, n, steps[1], steps[0]);
init_linearize_data(&b_in, nrhs, n, steps[3], steps[2]);
init_linearize_data(&r_out, nrhs, n, steps[5], steps[4]);
BEGIN_OUTER_LOOP_3
int not_ok;
linearize_@TYPE@_matrix(params.A, args[0], &a_in);
linearize_@TYPE@_matrix(params.B, args[1], &b_in);
not_ok =call_@lapack_func@(¶ms);
if (!not_ok) {
delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
} else {
error_occurred = 1;
nan_@TYPE@_matrix(args[2], &r_out);
}
END_OUTER_LOOP
release_@lapack_func@(¶ms);
}
set_fp_invalid_or_clear(error_occurred);
}
static void
@TYPE@_solve1(char **args, npy_intp *dimensions, npy_intp *steps,
void *NPY_UNUSED(func))
{
GESV_PARAMS_t params;
int error_occurred = get_fp_invalid_and_clear();
fortran_int n;
INIT_OUTER_LOOP_3
n = (fortran_int)dimensions[0];
if (init_@lapack_func@(¶ms, n, 1)) {
LINEARIZE_DATA_t a_in, b_in, r_out;
init_linearize_data(&a_in, n, n, steps[1], steps[0]);
init_linearize_data(&b_in, 1, n, 1, steps[2]);
init_linearize_data(&r_out, 1, n, 1, steps[3]);
BEGIN_OUTER_LOOP_3
int not_ok;
linearize_@TYPE@_matrix(params.A, args[0], &a_in);
linearize_@TYPE@_matrix(params.B, args[1], &b_in);
not_ok = call_@lapack_func@(¶ms);
if (!not_ok) {
delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
} else {
error_occurred = 1;
nan_@TYPE@_matrix(args[2], &r_out);
}
END_OUTER_LOOP
release_@lapack_func@(¶ms);
}
set_fp_invalid_or_clear(error_occurred);
}
static void
@TYPE@_inv(char **args, npy_intp *dimensions, npy_intp *steps,
void *NPY_UNUSED(func))
{
GESV_PARAMS_t params;
fortran_int n;
int error_occurred = get_fp_invalid_and_clear();
INIT_OUTER_LOOP_2
n = (fortran_int)dimensions[0];
if (init_@lapack_func@(¶ms, n, n)) {
LINEARIZE_DATA_t a_in, r_out;
init_linearize_data(&a_in, n, n, steps[1], steps[0]);
init_linearize_data(&r_out, n, n, steps[3], steps[2]);
BEGIN_OUTER_LOOP_2
int not_ok;
linearize_@TYPE@_matrix(params.A, args[0], &a_in);
identity_@TYPE@_matrix(params.B, n);
not_ok = call_@lapack_func@(¶ms);
if (!not_ok) {
delinearize_@TYPE@_matrix(args[1], params.B, &r_out);
} else {
error_occurred = 1;
nan_@TYPE@_matrix(args[1], &r_out);
}
END_OUTER_LOOP
release_@lapack_func@(¶ms);
}
set_fp_invalid_or_clear(error_occurred);
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* Cholesky decomposition */
typedef struct potr_params_struct
{
void *A;
fortran_int N;
fortran_int LDA;
char UPLO;
} POTR_PARAMS_t;
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#ftyp = fortran_real, fortran_doublereal,
fortran_complex, fortran_doublecomplex#
#lapack_func = spotrf, dpotrf, cpotrf, zpotrf#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(POTR_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->UPLO,
¶ms->N, params->A, ¶ms->LDA,
&rv);
return rv;
}
static NPY_INLINE int
init_@lapack_func@(POTR_PARAMS_t *params, char UPLO, fortran_int N)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *a;
size_t safe_N = N;
fortran_int lda = fortran_int_max(N, 1);
mem_buff = malloc(safe_N * safe_N * sizeof(@ftyp@));
if (!mem_buff) {
goto error;
}
a = mem_buff;
params->A = a;
params->N = N;
params->LDA = lda;
params->UPLO = UPLO;
return 1;
error:
free(mem_buff);
memset(params, 0, sizeof(*params));
return 0;
}
static NPY_INLINE void
release_@lapack_func@(POTR_PARAMS_t *params)
{
/* memory block base in A */
free(params->A);
memset(params, 0, sizeof(*params));
}
static void
@TYPE@_cholesky(char uplo, char **args, npy_intp *dimensions, npy_intp *steps)
{
POTR_PARAMS_t params;
int error_occurred = get_fp_invalid_and_clear();
fortran_int n;
INIT_OUTER_LOOP_2
assert(uplo == 'L');
n = (fortran_int)dimensions[0];
if (init_@lapack_func@(¶ms, uplo, n)) {
LINEARIZE_DATA_t a_in, r_out;
init_linearize_data(&a_in, n, n, steps[1], steps[0]);
init_linearize_data(&r_out, n, n, steps[3], steps[2]);
BEGIN_OUTER_LOOP_2
int not_ok;
linearize_@TYPE@_matrix(params.A, args[0], &a_in);
not_ok = call_@lapack_func@(¶ms);
if (!not_ok) {
triu_@TYPE@_matrix(params.A, params.N);
delinearize_@TYPE@_matrix(args[1], params.A, &r_out);
} else {
error_occurred = 1;
nan_@TYPE@_matrix(args[1], &r_out);
}
END_OUTER_LOOP
release_@lapack_func@(¶ms);
}
set_fp_invalid_or_clear(error_occurred);
}
static void
@TYPE@_cholesky_lo(char **args, npy_intp *dimensions, npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_cholesky('L', args, dimensions, steps);
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* eig family */
typedef struct geev_params_struct {
void *A;
void *WR; /* RWORK in complex versions, REAL W buffer for (sd)geev*/
void *WI;
void *VLR; /* REAL VL buffers for _geev where _ is s, d */
void *VRR; /* REAL VR buffers for _geev hwere _ is s, d */
void *WORK;
void *W; /* final w */
void *VL; /* final vl */
void *VR; /* final vr */
fortran_int N;
fortran_int LDA;
fortran_int LDVL;
fortran_int LDVR;
fortran_int LWORK;
char JOBVL;
char JOBVR;
} GEEV_PARAMS_t;
static NPY_INLINE void
dump_geev_params(const char *name, GEEV_PARAMS_t* params)
{
TRACE_TXT("\n%s\n"
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %p\n"\
"\t%10s: %d\n"\
"\t%10s: %d\n"\
"\t%10s: %d\n"\
"\t%10s: %d\n"\
"\t%10s: %d\n"\
"\t%10s: %c\n"\
"\t%10s: %c\n",
name,
"A", params->A,
"WR", params->WR,
"WI", params->WI,
"VLR", params->VLR,
"VRR", params->VRR,
"WORK", params->WORK,
"W", params->W,
"VL", params->VL,
"VR", params->VR,
"N", (int)params->N,
"LDA", (int)params->LDA,
"LDVL", (int)params->LDVL,
"LDVR", (int)params->LDVR,
"LWORK", (int)params->LWORK,
"JOBVL", params->JOBVL,
"JOBVR", params->JOBVR);
}
/**begin repeat
#TYPE = FLOAT, DOUBLE#
#CTYPE = CFLOAT, CDOUBLE#
#typ = float, double#
#complextyp = COMPLEX_t, DOUBLECOMPLEX_t#
#lapack_func = sgeev, dgeev#
#zero = 0.0f, 0.0#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(GEEV_PARAMS_t* params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBVL, ¶ms->JOBVR,
¶ms->N, params->A, ¶ms->LDA,
params->WR, params->WI,
params->VLR, ¶ms->LDVL,
params->VRR, ¶ms->LDVR,
params->WORK, ¶ms->LWORK,
&rv);
return rv;
}
static NPY_INLINE int
init_@lapack_func@(GEEV_PARAMS_t *params, char jobvl, char jobvr, fortran_int n)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *mem_buff2 = NULL;
npy_uint8 *a, *wr, *wi, *vlr, *vrr, *work, *w, *vl, *vr;
size_t safe_n = n;
size_t a_size = safe_n * safe_n * sizeof(@typ@);
size_t wr_size = safe_n * sizeof(@typ@);
size_t wi_size = safe_n * sizeof(@typ@);
size_t vlr_size = jobvl=='V' ? safe_n * safe_n * sizeof(@typ@) : 0;
size_t vrr_size = jobvr=='V' ? safe_n * safe_n * sizeof(@typ@) : 0;
size_t w_size = wr_size*2;
size_t vl_size = vlr_size*2;
size_t vr_size = vrr_size*2;
size_t work_count = 0;
fortran_int ld = fortran_int_max(n, 1);
/* allocate data for known sizes (all but work) */
mem_buff = malloc(a_size + wr_size + wi_size +
vlr_size + vrr_size +
w_size + vl_size + vr_size);
if (!mem_buff) {
goto error;
}
a = mem_buff;
wr = a + a_size;
wi = wr + wr_size;
vlr = wi + wi_size;
vrr = vlr + vlr_size;
w = vrr + vrr_size;
vl = w + w_size;
vr = vl + vl_size;
params->A = a;
params->WR = wr;
params->WI = wi;
params->VLR = vlr;
params->VRR = vrr;
params->W = w;
params->VL = vl;
params->VR = vr;
params->N = n;
params->LDA = ld;
params->LDVL = ld;
params->LDVR = ld;
params->JOBVL = jobvl;
params->JOBVR = jobvr;
/* Work size query */
{
@typ@ work_size_query;
params->LWORK = -1;
params->WORK = &work_size_query;
if (call_@lapack_func@(params) != 0) {
goto error;
}
work_count = (size_t)work_size_query;
}
mem_buff2 = malloc(work_count*sizeof(@typ@));
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
params->LWORK = (fortran_int)work_count;
params->WORK = work;
return 1;
error:
free(mem_buff2);
free(mem_buff);
memset(params, 0, sizeof(*params));
return 0;
}
static NPY_INLINE void
mk_@TYPE@_complex_array_from_real(@complextyp@ *c, const @typ@ *re, size_t n)
{
size_t iter;
for (iter = 0; iter < n; ++iter) {
c[iter].array[0] = re[iter];
c[iter].array[1] = @zero@;
}
}
static NPY_INLINE void
mk_@TYPE@_complex_array(@complextyp@ *c,
const @typ@ *re,
const @typ@ *im,
size_t n)
{
size_t iter;
for (iter = 0; iter < n; ++iter) {
c[iter].array[0] = re[iter];
c[iter].array[1] = im[iter];
}
}
static NPY_INLINE void
mk_@TYPE@_complex_array_conjugate_pair(@complextyp@ *c,
const @typ@ *r,
size_t n)
{
size_t iter;
for (iter = 0; iter < n; ++iter) {
@typ@ re = r[iter];
@typ@ im = r[iter+n];
c[iter].array[0] = re;
c[iter].array[1] = im;
c[iter+n].array[0] = re;
c[iter+n].array[1] = -im;
}
}
/*
* make the complex eigenvectors from the real array produced by sgeev/zgeev.
* c is the array where the results will be left.
* r is the source array of reals produced by sgeev/zgeev
* i is the eigenvalue imaginary part produced by sgeev/zgeev
* n is so that the order of the matrix is n by n
*/
static NPY_INLINE void
mk_@lapack_func@_complex_eigenvectors(@complextyp@ *c,
const @typ@ *r,
const @typ@ *i,
size_t n)
{
size_t iter = 0;
while (iter < n)
{
if (i[iter] == @zero@) {
/* eigenvalue was real, eigenvectors as well... */
mk_@TYPE@_complex_array_from_real(c, r, n);
c += n;
r += n;
iter ++;
} else {
/* eigenvalue was complex, generate a pair of eigenvectors */
mk_@TYPE@_complex_array_conjugate_pair(c, r, n);
c += 2*n;
r += 2*n;
iter += 2;
}
}
}
static NPY_INLINE void
process_@lapack_func@_results(GEEV_PARAMS_t *params)
{
/* REAL versions of geev need the results to be translated
* into complex versions. This is the way to deal with imaginary
* results. In our gufuncs we will always return complex arrays!
*/
mk_@TYPE@_complex_array(params->W, params->WR, params->WI, params->N);
/* handle the eigenvectors */
if ('V' == params->JOBVL) {
mk_@lapack_func@_complex_eigenvectors(params->VL, params->VLR,
params->WI, params->N);
}
if ('V' == params->JOBVR) {
mk_@lapack_func@_complex_eigenvectors(params->VR, params->VRR,
params->WI, params->N);
}
}
/**end repeat**/
/**begin repeat
#TYPE = CFLOAT, CDOUBLE#
#typ = COMPLEX_t, DOUBLECOMPLEX_t#
#ftyp = fortran_complex, fortran_doublecomplex#
#realtyp = float, double#
#lapack_func = cgeev, zgeev#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(GEEV_PARAMS_t* params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBVL, ¶ms->JOBVR,
¶ms->N, params->A, ¶ms->LDA,
params->W,
params->VL, ¶ms->LDVL,
params->VR, ¶ms->LDVR,
params->WORK, ¶ms->LWORK,
params->WR, /* actually RWORK */
&rv);
return rv;
}
static NPY_INLINE int
init_@lapack_func@(GEEV_PARAMS_t* params,
char jobvl,
char jobvr,
fortran_int n)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *mem_buff2 = NULL;
npy_uint8 *a, *w, *vl, *vr, *work, *rwork;
size_t safe_n = n;
size_t a_size = safe_n * safe_n * sizeof(@ftyp@);
size_t w_size = safe_n * sizeof(@ftyp@);
size_t vl_size = jobvl=='V'? safe_n * safe_n * sizeof(@ftyp@) : 0;
size_t vr_size = jobvr=='V'? safe_n * safe_n * sizeof(@ftyp@) : 0;
size_t rwork_size = 2 * safe_n * sizeof(@realtyp@);
size_t work_count = 0;
size_t total_size = a_size + w_size + vl_size + vr_size + rwork_size;
fortran_int ld = fortran_int_max(n, 1);
mem_buff = malloc(total_size);
if (!mem_buff) {
goto error;
}
a = mem_buff;
w = a + a_size;
vl = w + w_size;
vr = vl + vl_size;
rwork = vr + vr_size;
params->A = a;
params->WR = rwork;
params->WI = NULL;
params->VLR = NULL;
params->VRR = NULL;
params->VL = vl;
params->VR = vr;
params->W = w;
params->N = n;
params->LDA = ld;
params->LDVL = ld;
params->LDVR = ld;
params->JOBVL = jobvl;
params->JOBVR = jobvr;
/* Work size query */
{
@typ@ work_size_query;
params->LWORK = -1;
params->WORK = &work_size_query;
if (call_@lapack_func@(params) != 0) {
goto error;
}
work_count = (size_t) work_size_query.array[0];
/* Fix a bug in lapack 3.0.0 */
if(work_count == 0) work_count = 1;
}
mem_buff2 = malloc(work_count*sizeof(@ftyp@));
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
params->LWORK = (fortran_int)work_count;
params->WORK = work;
return 1;
error:
free(mem_buff2);
free(mem_buff);
memset(params, 0, sizeof(*params));
return 0;
}
static NPY_INLINE void
process_@lapack_func@_results(GEEV_PARAMS_t *NPY_UNUSED(params))
{
/* nothing to do here, complex versions are ready to copy out */
}
/**end repeat**/
/**begin repeat
#TYPE = FLOAT, DOUBLE, CDOUBLE#
#COMPLEXTYPE = CFLOAT, CDOUBLE, CDOUBLE#
#ftype = fortran_real, fortran_doublereal, fortran_doublecomplex#
#lapack_func = sgeev, dgeev, zgeev#
*/
static NPY_INLINE void
release_@lapack_func@(GEEV_PARAMS_t *params)
{
free(params->WORK);
free(params->A);
memset(params, 0, sizeof(*params));
}
static NPY_INLINE void
@TYPE@_eig_wrapper(char JOBVL,
char JOBVR,
char**args,
npy_intp* dimensions,
npy_intp* steps)
{
ptrdiff_t outer_steps[4];
size_t iter;
size_t outer_dim = *dimensions++;
size_t op_count = 2;
int error_occurred = get_fp_invalid_and_clear();
GEEV_PARAMS_t geev_params;
assert(JOBVL == 'N');
STACK_TRACE;
op_count += 'V'==JOBVL?1:0;
op_count += 'V'==JOBVR?1:0;
for (iter = 0; iter < op_count; ++iter) {
outer_steps[iter] = (ptrdiff_t) steps[iter];
}
steps += op_count;
if (init_@lapack_func@(&geev_params,
JOBVL, JOBVR,
(fortran_int)dimensions[0])) {
LINEARIZE_DATA_t a_in;
LINEARIZE_DATA_t w_out;
LINEARIZE_DATA_t vl_out;
LINEARIZE_DATA_t vr_out;
init_linearize_data(&a_in,
geev_params.N, geev_params.N,
steps[1], steps[0]);
steps += 2;
init_linearize_data(&w_out,
1, geev_params.N,
0, steps[0]);
steps += 1;
if ('V' == geev_params.JOBVL) {
init_linearize_data(&vl_out,
geev_params.N, geev_params.N,
steps[1], steps[0]);
steps += 2;
}
if ('V' == geev_params.JOBVR) {
init_linearize_data(&vr_out,
geev_params.N, geev_params.N,
steps[1], steps[0]);
}
for (iter = 0; iter < outer_dim; ++iter) {
int not_ok;
char **arg_iter = args;
/* copy the matrix in */
linearize_@TYPE@_matrix(geev_params.A, *arg_iter++, &a_in);
not_ok = call_@lapack_func@(&geev_params);
if (!not_ok) {
process_@lapack_func@_results(&geev_params);
delinearize_@COMPLEXTYPE@_matrix(*arg_iter++,
geev_params.W,
&w_out);
if ('V' == geev_params.JOBVL) {
delinearize_@COMPLEXTYPE@_matrix(*arg_iter++,
geev_params.VL,
&vl_out);
}
if ('V' == geev_params.JOBVR) {
delinearize_@COMPLEXTYPE@_matrix(*arg_iter++,
geev_params.VR,
&vr_out);
}
} else {
/* geev failed */
error_occurred = 1;
nan_@COMPLEXTYPE@_matrix(*arg_iter++, &w_out);
if ('V' == geev_params.JOBVL) {
nan_@COMPLEXTYPE@_matrix(*arg_iter++, &vl_out);
}
if ('V' == geev_params.JOBVR) {
nan_@COMPLEXTYPE@_matrix(*arg_iter++, &vr_out);
}
}
update_pointers((npy_uint8**)args, outer_steps, op_count);
}
release_@lapack_func@(&geev_params);
}
set_fp_invalid_or_clear(error_occurred);
}
static void
@TYPE@_eig(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_eig_wrapper('N', 'V', args, dimensions, steps);
}
static void
@TYPE@_eigvals(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_eig_wrapper('N', 'N', args, dimensions, steps);
}
/**end repeat**/
/* -------------------------------------------------------------------------- */
/* singular value decomposition */
typedef struct gessd_params_struct
{
void *A;
void *S;
void *U;
void *VT;
void *WORK;
void *RWORK;
void *IWORK;
fortran_int M;
fortran_int N;
fortran_int LDA;
fortran_int LDU;
fortran_int LDVT;
fortran_int LWORK;
char JOBZ;
} GESDD_PARAMS_t;
static NPY_INLINE void
dump_gesdd_params(const char *name,
GESDD_PARAMS_t *params)
{
TRACE_TXT("\n%s:\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18p\n"\
"%14s: %18d\n"\
"%14s: %18d\n"\
"%14s: %18d\n"\
"%14s: %18d\n"\
"%14s: %18d\n"\
"%14s: %18d\n"\
"%14s: %15c'%c'\n",
name,
"A", params->A,
"S", params->S,
"U", params->U,
"VT", params->VT,
"WORK", params->WORK,
"RWORK", params->RWORK,
"IWORK", params->IWORK,
"M", (int)params->M,
"N", (int)params->N,
"LDA", (int)params->LDA,
"LDU", (int)params->LDU,
"LDVT", (int)params->LDVT,
"LWORK", (int)params->LWORK,
"JOBZ", ' ', params->JOBZ);
}
static NPY_INLINE int
compute_urows_vtcolumns(char jobz,
fortran_int m, fortran_int n,
fortran_int *urows, fortran_int *vtcolumns)
{
fortran_int min_m_n = fortran_int_min(m, n);
switch(jobz)
{
case 'N':
*urows = 0;
*vtcolumns = 0;
break;
case 'A':
*urows = m;
*vtcolumns = n;
break;
case 'S':
{
*urows = min_m_n;
*vtcolumns = min_m_n;
}
break;
default:
return 0;
}
return 1;
}
/**begin repeat
#TYPE = FLOAT, DOUBLE#
#lapack_func = sgesdd, dgesdd#
#ftyp = fortran_real, fortran_doublereal#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(GESDD_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBZ, ¶ms->M, ¶ms->N,
params->A, ¶ms->LDA,
params->S,
params->U, ¶ms->LDU,
params->VT, ¶ms->LDVT,
params->WORK, ¶ms->LWORK,
params->IWORK,
&rv);
return rv;
}
static NPY_INLINE int
init_@lapack_func@(GESDD_PARAMS_t *params,
char jobz,
fortran_int m,
fortran_int n)
{
npy_uint8 *mem_buff = NULL;
npy_uint8 *mem_buff2 = NULL;
npy_uint8 *a, *s, *u, *vt, *work, *iwork;
size_t safe_m = m;
size_t safe_n = n;
size_t a_size = safe_m * safe_n * sizeof(@ftyp@);
fortran_int min_m_n = fortran_int_min(m, n);
size_t safe_min_m_n = min_m_n;
size_t s_size = safe_min_m_n * sizeof(@ftyp@);
fortran_int u_row_count, vt_column_count;
size_t safe_u_row_count, safe_vt_column_count;
size_t u_size, vt_size;
fortran_int work_count;
size_t work_size;
size_t iwork_size = 8 * safe_min_m_n * sizeof(fortran_int);
fortran_int ld = fortran_int_max(m, 1);
if (!compute_urows_vtcolumns(jobz, m, n, &u_row_count, &vt_column_count)) {
goto error;
}
safe_u_row_count = u_row_count;
safe_vt_column_count = vt_column_count;
u_size = safe_u_row_count * safe_m * sizeof(@ftyp@);
vt_size = safe_n * safe_vt_column_count * sizeof(@ftyp@);
mem_buff = malloc(a_size + s_size + u_size + vt_size + iwork_size);
if (!mem_buff) {
goto error;
}
a = mem_buff;
s = a + a_size;
u = s + s_size;
vt = u + u_size;
iwork = vt + vt_size;
/* fix vt_column_count so that it is a valid lapack parameter (0 is not) */
vt_column_count = fortran_int_max(1, vt_column_count);
params->M = m;
params->N = n;
params->A = a;
params->S = s;
params->U = u;
params->VT = vt;
params->RWORK = NULL;
params->IWORK = iwork;
params->M = m;
params->N = n;
params->LDA = ld;
params->LDU = ld;
params->LDVT = vt_column_count;
params->JOBZ = jobz;
/* Work size query */
{
@ftyp@ work_size_query;
params->LWORK = -1;
params->WORK = &work_size_query;
if (call_@lapack_func@(params) != 0) {
goto error;
}
work_count = (fortran_int)work_size_query;
/* Fix a bug in lapack 3.0.0 */
if(work_count == 0) work_count = 1;
work_size = (size_t)work_count * sizeof(@ftyp@);
}
mem_buff2 = malloc(work_size);
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
params->LWORK = work_count;
params->WORK = work;
return 1;
error:
TRACE_TXT("%s failed init\n", __FUNCTION__);
free(mem_buff);
free(mem_buff2);
memset(params, 0, sizeof(*params));
return 0;
}
/**end repeat**/
/**begin repeat
#TYPE = CFLOAT, CDOUBLE#
#ftyp = fortran_complex, fortran_doublecomplex#
#frealtyp = fortran_real, fortran_doublereal#
#typ = COMPLEX_t, DOUBLECOMPLEX_t#
#lapack_func = cgesdd, zgesdd#
*/
static NPY_INLINE fortran_int
call_@lapack_func@(GESDD_PARAMS_t *params)
{
fortran_int rv;
LAPACK(@lapack_func@)(¶ms->JOBZ, ¶ms->M, ¶ms->N,
params->A, ¶ms->LDA,
params->S,
params->U, ¶ms->LDU,
params->VT, ¶ms->LDVT,
params->WORK, ¶ms->LWORK,
params->RWORK,
params->IWORK,
&rv);
return rv;
}
static NPY_INLINE int
init_@lapack_func@(GESDD_PARAMS_t *params,
char jobz,
fortran_int m,
fortran_int n)
{
npy_uint8 *mem_buff = NULL, *mem_buff2 = NULL;
npy_uint8 *a,*s, *u, *vt, *work, *rwork, *iwork;
size_t a_size, s_size, u_size, vt_size, work_size, rwork_size, iwork_size;
size_t safe_u_row_count, safe_vt_column_count;
fortran_int u_row_count, vt_column_count, work_count;
size_t safe_m = m;
size_t safe_n = n;
fortran_int min_m_n = fortran_int_min(m, n);
size_t safe_min_m_n = min_m_n;
fortran_int ld = fortran_int_max(m, 1);
if (!compute_urows_vtcolumns(jobz, m, n, &u_row_count, &vt_column_count)) {
goto error;
}
safe_u_row_count = u_row_count;
safe_vt_column_count = vt_column_count;
a_size = safe_m * safe_n * sizeof(@ftyp@);
s_size = safe_min_m_n * sizeof(@frealtyp@);
u_size = safe_u_row_count * safe_m * sizeof(@ftyp@);
vt_size = safe_n * safe_vt_column_count * sizeof(@ftyp@);
rwork_size = 'N'==jobz?
(7 * safe_min_m_n) :
(5*safe_min_m_n * safe_min_m_n + 5*safe_min_m_n);
rwork_size *= sizeof(@ftyp@);
iwork_size = 8 * safe_min_m_n* sizeof(fortran_int);
mem_buff = malloc(a_size +
s_size +
u_size +
vt_size +
rwork_size +
iwork_size);
if (!mem_buff) {
goto error;
}
a = mem_buff;
s = a + a_size;
u = s + s_size;
vt = u + u_size;
rwork = vt + vt_size;
iwork = rwork + rwork_size;
/* fix vt_column_count so that it is a valid lapack parameter (0 is not) */
vt_column_count = fortran_int_max(1, vt_column_count);
params->A = a;
params->S = s;
params->U = u;
params->VT = vt;
params->RWORK = rwork;
params->IWORK = iwork;
params->M = m;
params->N = n;
params->LDA = ld;
params->LDU = ld;
params->LDVT = vt_column_count;
params->JOBZ = jobz;
/* Work size query */
{
@ftyp@ work_size_query;
params->LWORK = -1;
params->WORK = &work_size_query;
if (call_@lapack_func@(params) != 0) {
goto error;
}
work_count = (fortran_int)((@typ@*)&work_size_query)->array[0];
/* Fix a bug in lapack 3.0.0 */
if(work_count == 0) work_count = 1;
work_size = (size_t)work_count * sizeof(@ftyp@);
}
mem_buff2 = malloc(work_size);
if (!mem_buff2) {
goto error;
}
work = mem_buff2;
params->LWORK = work_count;
params->WORK = work;
return 1;
error:
TRACE_TXT("%s failed init\n", __FUNCTION__);
free(mem_buff2);
free(mem_buff);
memset(params, 0, sizeof(*params));
return 0;
}
/**end repeat**/
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
#REALTYPE = FLOAT, DOUBLE, FLOAT, DOUBLE#
#lapack_func = sgesdd, dgesdd, cgesdd, zgesdd#
*/
static NPY_INLINE void
release_@lapack_func@(GESDD_PARAMS_t* params)
{
/* A and WORK contain allocated blocks */
free(params->A);
free(params->WORK);
memset(params, 0, sizeof(*params));
}
static NPY_INLINE void
@TYPE@_svd_wrapper(char JOBZ,
char **args,
npy_intp* dimensions,
npy_intp* steps)
{
ptrdiff_t outer_steps[4];
int error_occurred = get_fp_invalid_and_clear();
size_t iter;
size_t outer_dim = *dimensions++;
size_t op_count = (JOBZ=='N')?2:4;
GESDD_PARAMS_t params;
for (iter = 0; iter < op_count; ++iter) {
outer_steps[iter] = (ptrdiff_t) steps[iter];
}
steps += op_count;
if (init_@lapack_func@(¶ms,
JOBZ,
(fortran_int)dimensions[0],
(fortran_int)dimensions[1])) {
LINEARIZE_DATA_t a_in, u_out, s_out, v_out;
init_linearize_data(&a_in, params.N, params.M, steps[1], steps[0]);
if ('N' == params.JOBZ) {
/* only the singular values are wanted */
fortran_int min_m_n = params.M < params.N? params.M : params.N;
init_linearize_data(&s_out, 1, min_m_n, 0, steps[2]);
} else {
fortran_int u_columns, v_rows;
fortran_int min_m_n = params.M < params.N? params.M : params.N;
if ('S' == params.JOBZ) {
u_columns = min_m_n;
v_rows = min_m_n;
} else {
u_columns = params.M;
v_rows = params.N;
}
init_linearize_data(&u_out,
u_columns, params.M,
steps[3], steps[2]);
init_linearize_data(&s_out,
1, min_m_n,
0, steps[4]);
init_linearize_data(&v_out,
params.N, v_rows,
steps[6], steps[5]);
}
for (iter = 0; iter < outer_dim; ++iter) {
int not_ok;
/* copy the matrix in */
linearize_@TYPE@_matrix(params.A, args[0], &a_in);
not_ok = call_@lapack_func@(¶ms);
if (!not_ok) {
if ('N' == params.JOBZ) {
delinearize_@REALTYPE@_matrix(args[1], params.S, &s_out);
} else {
delinearize_@TYPE@_matrix(args[1], params.U, &u_out);
delinearize_@REALTYPE@_matrix(args[2], params.S, &s_out);
delinearize_@TYPE@_matrix(args[3], params.VT, &v_out);
}
} else {
error_occurred = 1;
if ('N' == params.JOBZ) {
nan_@REALTYPE@_matrix(args[1], &s_out);
} else {
nan_@TYPE@_matrix(args[1], &u_out);
nan_@REALTYPE@_matrix(args[2], &s_out);
nan_@TYPE@_matrix(args[3], &v_out);
}
}
update_pointers((npy_uint8**)args, outer_steps, op_count);
}
release_@lapack_func@(¶ms);
}
set_fp_invalid_or_clear(error_occurred);
}
/**end repeat*/
/* svd gufunc entry points */
/**begin repeat
#TYPE = FLOAT, DOUBLE, CFLOAT, CDOUBLE#
*/
static void
@TYPE@_svd_N(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_svd_wrapper('N', args, dimensions, steps);
}
static void
@TYPE@_svd_S(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_svd_wrapper('S', args, dimensions, steps);
}
static void
@TYPE@_svd_A(char **args,
npy_intp *dimensions,
npy_intp *steps,
void *NPY_UNUSED(func))
{
@TYPE@_svd_wrapper('A', args, dimensions, steps);
}
/**end repeat**/
#pragma GCC diagnostic pop
/* -------------------------------------------------------------------------- */
/* gufunc registration */
static void *array_of_nulls[] = {
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL,
(void *)NULL
};
#define FUNC_ARRAY_NAME(NAME) NAME ## _funcs
#define GUFUNC_FUNC_ARRAY_REAL(NAME) \
static PyUFuncGenericFunction \
FUNC_ARRAY_NAME(NAME)[] = { \
FLOAT_ ## NAME, \
DOUBLE_ ## NAME \
}
#define GUFUNC_FUNC_ARRAY_REAL_COMPLEX(NAME) \
static PyUFuncGenericFunction \
FUNC_ARRAY_NAME(NAME)[] = { \
FLOAT_ ## NAME, \
DOUBLE_ ## NAME, \
CFLOAT_ ## NAME, \
CDOUBLE_ ## NAME \
}
/* There are problems with eig in complex single precision.
* That kernel is disabled
*/
#define GUFUNC_FUNC_ARRAY_EIG(NAME) \
static PyUFuncGenericFunction \
FUNC_ARRAY_NAME(NAME)[] = { \
FLOAT_ ## NAME, \
DOUBLE_ ## NAME, \
CDOUBLE_ ## NAME \
}
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(slogdet);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(det);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(eighlo);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(eighup);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(eigvalshlo);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(eigvalshup);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(solve);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(solve1);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(inv);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(cholesky_lo);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(svd_N);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(svd_S);
GUFUNC_FUNC_ARRAY_REAL_COMPLEX(svd_A);
GUFUNC_FUNC_ARRAY_EIG(eig);
GUFUNC_FUNC_ARRAY_EIG(eigvals);
static char equal_2_types[] = {
NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_CFLOAT,
NPY_CDOUBLE, NPY_CDOUBLE
};
static char equal_3_types[] = {
NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_CFLOAT, NPY_CFLOAT,
NPY_CDOUBLE, NPY_CDOUBLE, NPY_CDOUBLE
};
/* second result is logdet, that will always be a REAL */
static char slogdet_types[] = {
NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_CFLOAT, NPY_FLOAT,
NPY_CDOUBLE, NPY_CDOUBLE, NPY_DOUBLE
};
static char eigh_types[] = {
NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_FLOAT, NPY_CFLOAT,
NPY_CDOUBLE, NPY_DOUBLE, NPY_CDOUBLE
};
static char eighvals_types[] = {
NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_FLOAT,
NPY_CDOUBLE, NPY_DOUBLE
};
static char eig_types[] = {
NPY_FLOAT, NPY_CFLOAT, NPY_CFLOAT,
NPY_DOUBLE, NPY_CDOUBLE, NPY_CDOUBLE,
NPY_CDOUBLE, NPY_CDOUBLE, NPY_CDOUBLE
};
static char eigvals_types[] = {
NPY_FLOAT, NPY_CFLOAT,
NPY_DOUBLE, NPY_CDOUBLE,
NPY_CDOUBLE, NPY_CDOUBLE
};
static char svd_1_1_types[] = {
NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_FLOAT,
NPY_CDOUBLE, NPY_DOUBLE
};
static char svd_1_3_types[] = {
NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
NPY_CFLOAT, NPY_CFLOAT, NPY_FLOAT, NPY_CFLOAT,
NPY_CDOUBLE, NPY_CDOUBLE, NPY_DOUBLE, NPY_CDOUBLE
};
typedef struct gufunc_descriptor_struct {
char *name;
char *signature;
char *doc;
int ntypes;
int nin;
int nout;
PyUFuncGenericFunction *funcs;
char *types;
} GUFUNC_DESCRIPTOR_t;
GUFUNC_DESCRIPTOR_t gufunc_descriptors [] = {
{
"slogdet",
"(m,m)->(),()",
"slogdet on the last two dimensions and broadcast on the rest. \n"\
"Results in two arrays, one with sign and the other with log of the"\
" determinants. \n"\
" \"(m,m)->(),()\" \n",
4, 1, 2,
FUNC_ARRAY_NAME(slogdet),
slogdet_types
},
{
"det",
"(m,m)->()",
"det of the last two dimensions and broadcast on the rest. \n"\
" \"(m,m)->()\" \n",
4, 1, 1,
FUNC_ARRAY_NAME(det),
equal_2_types
},
{
"eigh_lo",
"(m,m)->(m),(m,m)",
"eigh on the last two dimension and broadcast to the rest, using"\
" lower triangle \n"\
"Results in a vector of eigenvalues and a matrix with the"\
"eigenvectors. \n"\
" \"(m,m)->(m),(m,m)\" \n",
4, 1, 2,
FUNC_ARRAY_NAME(eighlo),
eigh_types
},
{
"eigh_up",
"(m,m)->(m),(m,m)",
"eigh on the last two dimension and broadcast to the rest, using"\
" upper triangle. \n"\
"Results in a vector of eigenvalues and a matrix with the"\
" eigenvectors. \n"\
" \"(m,m)->(m),(m,m)\" \n",
4, 1, 2,
FUNC_ARRAY_NAME(eighup),
eigh_types
},
{
"eigvalsh_lo",
"(m,m)->(m)",
"eigh on the last two dimension and broadcast to the rest, using"\
" lower triangle. \n"\
"Results in a vector of eigenvalues and a matrix with the"\
"eigenvectors. \n"\
" \"(m,m)->(m)\" \n",
4, 1, 1,
FUNC_ARRAY_NAME(eigvalshlo),
eighvals_types
},
{
"eigvalsh_up",
"(m,m)->(m)",
"eigvalsh on the last two dimension and broadcast to the rest,"\
" using upper triangle. \n"\
"Results in a vector of eigenvalues and a matrix with the"\
"eigenvectors.\n"\
" \"(m,m)->(m)\" \n",
4, 1, 1,
FUNC_ARRAY_NAME(eigvalshup),
eighvals_types
},
{
"solve",
"(m,m),(m,n)->(m,n)",
"solve the system a x = b, on the last two dimensions, broadcast"\
" to the rest. \n"\
"Results in a matrices with the solutions. \n"\
" \"(m,m),(m,n)->(m,n)\" \n",
4, 2, 1,
FUNC_ARRAY_NAME(solve),
equal_3_types
},
{
"solve1",
"(m,m),(m)->(m)",
"solve the system a x = b, for b being a vector, broadcast in"\
" the outer dimensions. \n"\
"Results in vectors with the solutions. \n"\
" \"(m,m),(m)->(m)\" \n",
4, 2, 1,
FUNC_ARRAY_NAME(solve1),
equal_3_types
},
{
"inv",
"(m, m)->(m, m)",
"compute the inverse of the last two dimensions and broadcast"\
" to the rest. \n"\
"Results in the inverse matrices. \n"\
" \"(m,m)->(m,m)\" \n",
4, 1, 1,
FUNC_ARRAY_NAME(inv),
equal_2_types
},
{
"cholesky_lo",
"(m,m)->(m,m)",
"cholesky decomposition of hermitian positive-definite matrices. \n"\
"Broadcast to all outer dimensions. \n"\
" \"(m,m)->(m,m)\" \n",
4, 1, 1,
FUNC_ARRAY_NAME(cholesky_lo),
equal_2_types
},
{
"svd_m",
"(m,n)->(m)",
"svd when n>=m. ",
4, 1, 1,
FUNC_ARRAY_NAME(svd_N),
svd_1_1_types
},
{
"svd_n",
"(m,n)->(n)",
"svd when n<=m",
4, 1, 1,
FUNC_ARRAY_NAME(svd_N),
svd_1_1_types
},
{
"svd_m_s",
"(m,n)->(m,m),(m),(m,n)",
"svd when m<=n",
4, 1, 3,
FUNC_ARRAY_NAME(svd_S),
svd_1_3_types
},
{
"svd_n_s",
"(m,n)->(m,n),(n),(n,n)",
"svd when m>=n",
4, 1, 3,
FUNC_ARRAY_NAME(svd_S),
svd_1_3_types
},
{
"svd_m_f",
"(m,n)->(m,m),(m),(n,n)",
"svd when m<=n",
4, 1, 3,
FUNC_ARRAY_NAME(svd_A),
svd_1_3_types
},
{
"svd_n_f",
"(m,n)->(m,m),(n),(n,n)",
"svd when m>=n",
4, 1, 3,
FUNC_ARRAY_NAME(svd_A),
svd_1_3_types
},
{
"eig",
"(m,m)->(m),(m,m)",
"eig on the last two dimension and broadcast to the rest. \n"\
"Results in a vector with the eigenvalues and a matrix with the"\
" eigenvectors. \n"\
" \"(m,m)->(m),(m,m)\" \n",
3, 1, 2,
FUNC_ARRAY_NAME(eig),
eig_types
},
{
"eigvals",
"(m,m)->(m)",
"eigvals on the last two dimension and broadcast to the rest. \n"\
"Results in a vector of eigenvalues. \n"\
" \"(m,m)->(m),(m,m)\" \n",
3, 1, 1,
FUNC_ARRAY_NAME(eigvals),
eigvals_types
},
};
static void
addUfuncs(PyObject *dictionary) {
PyObject *f;
int i;
const int gufunc_count = sizeof(gufunc_descriptors)/
sizeof(gufunc_descriptors[0]);
for (i = 0; i < gufunc_count; i++) {
GUFUNC_DESCRIPTOR_t* d = &gufunc_descriptors[i];
f = PyUFunc_FromFuncAndDataAndSignature(d->funcs,
array_of_nulls,
d->types,
d->ntypes,
d->nin,
d->nout,
PyUFunc_None,
d->name,
d->doc,
0,
d->signature);
PyDict_SetItemString(dictionary, d->name, f);
#if 0
dump_ufunc_object((PyUFuncObject*) f);
#endif
Py_DECREF(f);
}
}
/* -------------------------------------------------------------------------- */
/* Module initialization stuff */
static PyMethodDef UMath_LinAlgMethods[] = {
{NULL, NULL, 0, NULL} /* Sentinel */
};
#if defined(NPY_PY3K)
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
UMATH_LINALG_MODULE_NAME,
NULL,
-1,
UMath_LinAlgMethods,
NULL,
NULL,
NULL,
NULL
};
#endif
#if defined(NPY_PY3K)
#define RETVAL m
PyObject *PyInit__umath_linalg(void)
#else
#define RETVAL
PyMODINIT_FUNC
init_umath_linalg(void)
#endif
{
PyObject *m;
PyObject *d;
PyObject *version;
init_constants();
#if defined(NPY_PY3K)
m = PyModule_Create(&moduledef);
#else
m = Py_InitModule(UMATH_LINALG_MODULE_NAME, UMath_LinAlgMethods);
#endif
if (m == NULL) {
return RETVAL;
}
import_array();
import_ufunc();
d = PyModule_GetDict(m);
version = PyString_FromString(umath_linalg_version_string);
PyDict_SetItemString(d, "__version__", version);
Py_DECREF(version);
/* Load the ufunc operators into the module's namespace */
addUfuncs(d);
if (PyErr_Occurred()) {
PyErr_SetString(PyExc_RuntimeError,
"cannot load _umath_linalg module.");
}
return RETVAL;
}