7 #ifndef __eigenpy_ufunc_hpp__
8 #define __eigenpy_ufunc_hpp__
16 #ifdef NPY_1_19_API_VERSION
17 #define EIGENPY_NPY_CONST_UFUNC_ARG const
19 #define EIGENPY_NPY_CONST_UFUNC_ARG
24 npy_intp
const *steps) {
31 npy_intp dm = dimensions[0];
32 npy_intp dn = dimensions[1];
33 npy_intp dp = dimensions[2];
36 npy_intp is1_m = steps[0];
37 npy_intp is1_n = steps[1];
38 npy_intp is2_n = steps[2];
39 npy_intp is2_p = steps[3];
40 npy_intp os_m = steps[4];
41 npy_intp os_p = steps[5];
47 for (m = 0; m < dm; m++) {
48 for (p = 0; p < dp; p++) {
70 void *NPY_UNUSED(func)) {
75 npy_intp dN = dimensions[0];
78 npy_intp s0 = steps[0];
79 npy_intp s1 = steps[1];
80 npy_intp s2 = steps[2];
86 for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) {
87 matrix_multiply<T>(args, dimensions + 1, steps + 3);
91 #define EIGENPY_REGISTER_BINARY_OPERATOR(name, op) \
92 template <typename T1, typename T2, typename R> \
93 void binary_op_##name( \
94 char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, \
95 EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * ) { \
96 npy_intp is0 = steps[0], is1 = steps[1], os = steps[2], n = *dimensions; \
97 char *i0 = args[0], *i1 = args[1], *o = args[2]; \
99 for (k = 0; k < n; k++) { \
100 T1 &x = *static_cast<T1 *>(static_cast<void *>(i0)); \
101 T2 &y = *static_cast<T2 *>(static_cast<void *>(i1)); \
102 R &res = *static_cast<R *>(static_cast<void *>(o)); \
110 template <typename T> \
111 void binary_op_##name( \
112 char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, \
113 EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void *data) { \
114 binary_op_##name<T, T, T>(args, dimensions, steps, data); \
128 #define EIGENPY_REGISTER_UNARY_OPERATOR(name, op) \
129 template <typename T, typename R> \
130 void unary_op_##name( \
131 char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, \
132 EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void * ) { \
133 npy_intp is = steps[0], os = steps[1], n = *dimensions; \
134 char *i = args[0], *o = args[1]; \
136 for (k = 0; k < n; k++) { \
137 T &x = *static_cast<T *>(static_cast<void *>(i)); \
138 R &res = *static_cast<R *>(static_cast<void *>(o)); \
145 template <typename T> \
146 void unary_op_##name( \
147 char **args, EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *dimensions, \
148 EIGENPY_NPY_CONST_UFUNC_ARG npy_intp *steps, void *data) { \
149 unary_op_##name<T, T>(args, dimensions, steps, data); \
156 #define EIGENPY_REGISTER_BINARY_UFUNC(name, code, T1, T2, R) \
158 PyUFuncObject *ufunc = \
159 (PyUFuncObject *)PyObject_GetAttrString(numpy, #name); \
160 int _types[3] = {Register::getTypeCode<T1>(), Register::getTypeCode<T2>(), \
161 Register::getTypeCode<R>()}; \
165 if (sizeof(_types) / sizeof(int) != ufunc->nargs) { \
166 PyErr_Format(PyExc_AssertionError, \
167 "ufunc %s takes %d arguments, our loop takes %lu", #name, \
169 (unsigned long)(sizeof(_types) / sizeof(int))); \
172 if (PyUFunc_RegisterLoopForType((PyUFuncObject *)ufunc, code, \
173 internal::binary_op_##name<T1, T2, R>, \
181 #define EIGENPY_REGISTER_UNARY_UFUNC(name, code, T, R) \
183 PyUFuncObject *ufunc = \
184 (PyUFuncObject *)PyObject_GetAttrString(numpy, #name); \
185 int _types[2] = {Register::getTypeCode<T>(), Register::getTypeCode<R>()}; \
189 if (sizeof(_types) / sizeof(int) != ufunc->nargs) { \
190 PyErr_Format(PyExc_AssertionError, \
191 "ufunc %s takes %d arguments, our loop takes %lu", #name, \
193 (unsigned long)(sizeof(_types) / sizeof(int))); \
196 if (PyUFunc_RegisterLoopForType((PyUFuncObject *)ufunc, code, \
197 internal::unary_op_##name<T, R>, _types, \
205 template <
typename Scalar>
207 const int type_code = Register::getTypeCode<Scalar>();
210 #if PY_MAJOR_VERSION >= 3
211 numpy_str = PyUnicode_FromString(
"numpy");
213 numpy_str = PyString_FromString(
"numpy");
216 numpy = PyImport_Import(numpy_str);
217 Py_DECREF(numpy_str);
223 int types[3] = {type_code, type_code, type_code};
225 std::stringstream ss;
226 ss <<
"return result of multiplying two matrices of ";
227 ss << bp::type_info(
typeid(Scalar)).name();
228 PyUFuncObject *ufunc =
229 (PyUFuncObject *)PyObject_GetAttrString(numpy,
"matmul");
231 std::stringstream ss;
232 ss <<
"Impossible to define matrix_multiply for given type "
233 << bp::type_info(
typeid(Scalar)).name() << std::endl;
236 if (PyUFunc_RegisterLoopForType((PyUFuncObject *)ufunc, type_code,
237 &internal::gufunc_matrix_multiply<Scalar>,
239 std::stringstream ss;
240 ss <<
"Impossible to register matrix_multiply for given type "
241 << bp::type_info(
typeid(Scalar)).name() << std::endl;
270 #endif // __eigenpy_ufunc_hpp__