| 1195 | |
| 1196 | template<typename typ, typename basetyp> |
| 1197 | static void |
| 1198 | det(char **args, |
| 1199 | npy_intp const *dimensions, |
| 1200 | npy_intp const *steps, |
| 1201 | void *NPY_UNUSED(func)) |
| 1202 | { |
| 1203 | fortran_int m; |
| 1204 | char *tmp_buff; |
| 1205 | size_t matrix_size; |
| 1206 | size_t pivot_size; |
| 1207 | size_t safe_m; |
| 1208 | /* notes: |
| 1209 | * matrix will need to be copied always, as factorization in lapack is |
| 1210 | * made inplace |
| 1211 | * matrix will need to be in column-major order, as expected by lapack |
| 1212 | * code (fortran) |
| 1213 | * always a square matrix |
| 1214 | * need to allocate memory for both, matrix_buffer and pivot buffer |
| 1215 | */ |
| 1216 | INIT_OUTER_LOOP_2 |
| 1217 | m = (fortran_int) dimensions[0]; |
| 1218 | /* avoid empty malloc (buffers likely unused) and ensure m is `size_t` */ |
| 1219 | safe_m = m != 0 ? m : 1; |
| 1220 | matrix_size = safe_m * safe_m * sizeof(typ); |
| 1221 | pivot_size = safe_m * sizeof(fortran_int); |
| 1222 | tmp_buff = (char *)malloc(matrix_size + pivot_size); |
| 1223 | |
| 1224 | if (tmp_buff) { |
| 1225 | LINEARIZE_DATA_t lin_data; |
| 1226 | typ sign; |
| 1227 | basetyp logdet; |
| 1228 | /* swapped steps to get matrix in FORTRAN order */ |
| 1229 | init_linearize_data(&lin_data, m, m, steps[1], steps[0]); |
| 1230 | |
| 1231 | BEGIN_OUTER_LOOP_2 |
| 1232 | linearize_matrix((typ*)tmp_buff, (typ*)args[0], &lin_data); |
| 1233 | slogdet_single_element(m, |
| 1234 | (typ*)tmp_buff, |
| 1235 | (fortran_int*)(tmp_buff + matrix_size), |
| 1236 | &sign, |
| 1237 | &logdet); |
| 1238 | *(typ *)args[1] = det_from_slogdet(sign, logdet); |
| 1239 | END_OUTER_LOOP |
| 1240 | |
| 1241 | free(tmp_buff); |
| 1242 | } |
| 1243 | else { |
| 1244 | /* TODO: Requires use of new ufunc API to indicate error return */ |
| 1245 | NPY_ALLOW_C_API_DEF |
| 1246 | NPY_ALLOW_C_API; |
| 1247 | PyErr_NoMemory(); |
| 1248 | NPY_DISABLE_C_API; |
| 1249 | } |
| 1250 | } |
| 1251 | |
| 1252 | |
| 1253 | /* -------------------------------------------------------------------------- */ |
nothing calls this directly
no test coverage detected