From c000e129b1521e7b9b89e178d011d25582fdbb08 Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Fri, 24 Mar 2017 21:02:37 +0000 Subject: MAINT: Always use parameter structs to call lapack Getting the arguments to lapack calls in the correct order in two places is harder than doing it in one. --- numpy/linalg/umath_linalg.c.src | 294 +++++++++++++++++++++------------------- 1 file changed, 153 insertions(+), 141 deletions(-) (limited to 'numpy/linalg/umath_linalg.c.src') diff --git a/numpy/linalg/umath_linalg.c.src b/numpy/linalg/umath_linalg.c.src index 2073ce3d5..aadb697ea 100644 --- a/numpy/linalg/umath_linalg.c.src +++ b/numpy/linalg/umath_linalg.c.src @@ -1280,11 +1280,8 @@ init_@lapack_func@(EIGH_PARAMS_t* params, char JOBZ, char UPLO, { npy_uint8 *mem_buff = NULL; npy_uint8 *mem_buff2 = NULL; - @typ@ query_work_size; - fortran_int query_iwork_size; - fortran_int lwork = -1; - fortran_int liwork = -1; - fortran_int info; + 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@); @@ -1296,18 +1293,33 @@ init_@lapack_func@(EIGH_PARAMS_t* params, char JOBZ, char UPLO, goto error; a = mem_buff; w = mem_buff + safe_N * safe_N * sizeof(@typ@); - LAPACK(@lapack_func@)(&JOBZ, &UPLO, &N, - (@ftyp@*)a, &lda, (@ftyp@*)w, - &query_work_size, &lwork, - &query_iwork_size, &liwork, - &info); - if (info != 0) - goto error; + 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; + } - work = mem_buff; - 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; @@ -1315,18 +1327,10 @@ init_@lapack_func@(EIGH_PARAMS_t* params, char JOBZ, char UPLO, work = mem_buff2; iwork = mem_buff2 + lwork*sizeof(@typ@); - params->A = a; - params->W = w; - params->WORK = work; - params->RWORK = NULL; /* unused */ - params->IWORK = iwork; - params->N = N; params->LWORK = lwork; - params->LRWORK = 0; /* unused */ + params->WORK = work; params->LIWORK = liwork; - params->JOBZ = JOBZ; - params->UPLO = UPLO; - params->LDA = lda; + params->IWORK = iwork; return 1; @@ -1349,7 +1353,6 @@ init_@lapack_func@(EIGH_PARAMS_t* params, char JOBZ, char UPLO, #fbasetyp=fortran_real,fortran_doublereal# #lapack_func=cheevd,zheevd# */ - static inline fortran_int call_@lapack_func@(EIGH_PARAMS_t *params) { @@ -1375,14 +1378,10 @@ init_@lapack_func@(EIGH_PARAMS_t *params, { npy_uint8 *mem_buff = NULL; npy_uint8 *mem_buff2 = NULL; - @ftyp@ query_work_size; - @fbasetyp@ query_rwork_size; - fortran_int query_iwork_size; - fortran_int lwork = -1; - fortran_int lrwork = -1; - fortran_int liwork = -1; + fortran_int lwork; + fortran_int lrwork; + fortran_int liwork; npy_uint8 *a, *w, *work, *rwork, *iwork; - fortran_int info; size_t safe_N = N; fortran_int lda = fortran_int_max(N, 1); @@ -1393,40 +1392,50 @@ init_@lapack_func@(EIGH_PARAMS_t *params, a = mem_buff; w = mem_buff + safe_N * safe_N * sizeof(@typ@); - LAPACK(@lapack_func@)(&JOBZ, &UPLO, &N, - (@ftyp@*)a, &lda, (@fbasetyp@*)w, - &query_work_size, &lwork, - &query_rwork_size, &lrwork, - &query_iwork_size, &liwork, - &info); - if (info != 0) - goto error; + 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; + 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->A = a; - params->W = w; params->WORK = work; params->RWORK = rwork; params->IWORK = iwork; - params->N = N; params->LWORK = lwork; params->LRWORK = lrwork; params->LIWORK = liwork; - params->JOBZ = JOBZ; - params->UPLO = UPLO; - params->LDA = lda; return 1; @@ -1984,9 +1993,6 @@ init_@lapack_func@(GEEV_PARAMS_t *params, char jobvl, char jobvr, fortran_int n) size_t vl_size = vlr_size*2; size_t vr_size = vrr_size*2; size_t work_count = 0; - @typ@ work_size_query; - fortran_int do_size_query = -1; - fortran_int rv; fortran_int ld = fortran_int_max(n, 1); /* allocate data for known sizes (all but work) */ @@ -2004,27 +2010,12 @@ init_@lapack_func@(GEEV_PARAMS_t *params, char jobvl, char jobvr, fortran_int n) w = vrr + vrr_size; vl = w + w_size; vr = vl + vl_size; - LAPACK(@lapack_func@)(&jobvl, &jobvr, &n, - (void *)a, &ld, (void *)wr, (void *)wi, - (void *)vl, &ld, (void *)vr, &ld, - &work_size_query, &do_size_query, - &rv); - - if (0 != rv) - 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->A = a; params->WR = wr; params->WI = wi; params->VLR = vlr; params->VRR = vrr; - params->WORK = work; params->W = w; params->VL = vl; params->VR = vr; @@ -2032,10 +2023,30 @@ init_@lapack_func@(GEEV_PARAMS_t *params, char jobvl, char jobvr, fortran_int n) params->LDA = ld; params->LDVL = ld; params->LDVR = ld; - params->LWORK = (fortran_int)work_count; 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); @@ -2180,9 +2191,6 @@ init_@lapack_func@(GEEV_PARAMS_t* params, 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; - @typ@ work_size_query; - fortran_int do_size_query = -1; - fortran_int rv; size_t total_size = a_size + w_size + vl_size + vr_size + rwork_size; fortran_int ld = fortran_int_max(n, 1); @@ -2195,24 +2203,6 @@ init_@lapack_func@(GEEV_PARAMS_t* params, vl = w + w_size; vr = vl + vl_size; rwork = vr + vr_size; - LAPACK(@lapack_func@)(&jobvl, &jobvr, &n, - (void *)a, &ld, (void *)w, - (void *)vl, &ld, (void *)vr, &ld, - (void *)&work_size_query, &do_size_query, - (void *)rwork, - &rv); - if (0 != rv) - 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->A = a; params->WR = rwork; @@ -2221,16 +2211,38 @@ init_@lapack_func@(GEEV_PARAMS_t* params, params->VRR = NULL; params->VL = vl; params->VR = vr; - params->WORK = work; params->W = w; params->N = n; params->LDA = ld; params->LDVL = ld; params->LDVR = ld; - params->LWORK = (fortran_int)work_count; 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); @@ -2240,6 +2252,7 @@ init_@lapack_func@(GEEV_PARAMS_t* params, return 0; } + static inline void process_@lapack_func@_results(GEEV_PARAMS_t *NPY_UNUSED(params)) { @@ -2537,19 +2550,33 @@ init_@lapack_func@(GESDD_PARAMS_t *params, iwork = vt + vt_size; /* fix vt_column_count so that it is a valid lapack parameter (0 is not) */ - vt_column_count = vt_column_count < 1? 1 : vt_column_count; + 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 */ { - /* compute optimal work size */ @ftyp@ work_size_query; - fortran_int do_query = -1; - fortran_int rv; - LAPACK(@lapack_func@)(&jobz, &m, &n, - (void*)a, &ld, (void*)s, (void*)u, &ld, - (void*)vt, &vt_column_count, - &work_size_query, &do_query, - (void*)iwork, &rv); - if (0!=rv) + + 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; @@ -2562,22 +2589,8 @@ init_@lapack_func@(GESDD_PARAMS_t *params, work = mem_buff2; - params->M = m; - params->N = n; - params->A = a; - params->S = s; - params->U = u; - params->VT = vt; - params->WORK = work; - params->RWORK = NULL; - params->IWORK = iwork; - params->M = m; - params->N = n; - params->LDA = ld; - params->LDU = ld; - params->LDVT = vt_column_count; params->LWORK = work_count; - params->JOBZ = jobz; + params->WORK = work; return 1; error: @@ -2665,20 +2678,31 @@ init_@lapack_func@(GESDD_PARAMS_t *params, iwork = rwork + rwork_size; /* fix vt_column_count so that it is a valid lapack parameter (0 is not) */ - vt_column_count = vt_column_count < 1? 1 : vt_column_count; + 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 */ { - /* compute optimal work size */ @ftyp@ work_size_query; - fortran_int do_query = -1; - fortran_int rv; - LAPACK(@lapack_func@)(&jobz, &m, &n, - (void*)a, &ld, (void*)s, (void*)u, &ld, - (void*)vt, &vt_column_count, - &work_size_query, &do_query, - (void*)rwork, - (void*)iwork, &rv); - if (0!=rv) + + 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; @@ -2691,20 +2715,8 @@ init_@lapack_func@(GESDD_PARAMS_t *params, work = mem_buff2; - params->A = a; - params->S = s; - params->U = u; - params->VT = vt; - params->WORK = work; - params->RWORK = rwork; - params->IWORK = iwork; - params->M = m; - params->N = n; - params->LDA = ld; - params->LDU = ld; - params->LDVT = vt_column_count; params->LWORK = work_count; - params->JOBZ = jobz; + params->WORK = work; return 1; error: -- cgit v1.2.1