Skip to content

Commit

Permalink
MAINT: spatial: use cython_lapack in spatial/_qhull.pyx
Browse files Browse the repository at this point in the history
  • Loading branch information
ev-br authored and tylerjereddy committed Apr 1, 2024
1 parent 729ff0f commit abb04b2
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 27 deletions.
26 changes: 8 additions & 18 deletions scipy/spatial/_qhull.pyx
Expand Up @@ -23,6 +23,8 @@ from libc.math cimport NAN
from scipy._lib.messagestream cimport MessageStream
from libc.stdio cimport FILE

from scipy.linalg.cython_lapack cimport dgetrf, dgetrs, dgecon

np.import_array()


Expand Down Expand Up @@ -170,6 +172,8 @@ cdef extern from "qhull_src/src/libqhull_r.h":
coordT* qh_sethalfspace_all(qhT *, int dim, int count, coordT* halfspaces, pointT *feasible)

cdef extern from "qhull_misc.h":
ctypedef int CBLAS_INT # actual type defined in the header file
void qhull_misc_lib_check()
int qh_new_qhull_scipy(qhT *, int dim, int numpoints, realT *points,
boolT ismalloc, char* qhull_cmd, void *outfile,
void *errfile, coordT* feaspoint) nogil
Expand Down Expand Up @@ -201,20 +205,6 @@ cdef extern from "qhull_src/src/mem_r.h":

from libc.stdlib cimport qsort

#------------------------------------------------------------------------------
# LAPACK interface
#------------------------------------------------------------------------------

cdef extern from "qhull_misc.h":
ctypedef int CBLAS_INT # actual type defined in the header file
void qhull_misc_lib_check()
void qh_dgetrf(CBLAS_INT *m, CBLAS_INT *n, double *a, CBLAS_INT *lda, CBLAS_INT *ipiv,
CBLAS_INT *info) nogil
void qh_dgetrs(char *trans, CBLAS_INT *n, CBLAS_INT *nrhs, double *a, CBLAS_INT *lda,
CBLAS_INT *ipiv, double *b, CBLAS_INT *ldb, CBLAS_INT *info) nogil
void qh_dgecon(char *norm, CBLAS_INT *n, double *a, CBLAS_INT *lda, double *anorm,
double *rcond, double *work, CBLAS_INT *iwork, CBLAS_INT *info) nogil


#------------------------------------------------------------------------------
# Qhull wrapper
Expand Down Expand Up @@ -1113,20 +1103,20 @@ def _get_barycentric_transforms(np.ndarray[np.double_t, ndim=2] points,
nrhs = ndim
lda = ndim
ldb = ndim
qh_dgetrf(&n, &n, <double*>T.data, &lda, ipiv, &info)
dgetrf(&n, &n, <double*>T.data, &lda, ipiv, &info)

# Check condition number
if info == 0:
qh_dgecon("1", &n, <double*>T.data, &lda, &anorm, &rcond,
work, iwork, &info)
dgecon("1", &n, <double*>T.data, &lda, &anorm, &rcond,
work, iwork, &info)

if rcond < rcond_limit:
# The transform seems singular
info = 1

# Compute transform
if info == 0:
qh_dgetrs("N", &n, &nrhs, <double*>T.data, &lda, ipiv,
dgetrs("N", &n, &nrhs, <double*>T.data, &lda, ipiv,
(<double*>Tinvs.data) + ndim*(ndim+1)*isimplex,
&ldb, &info)

Expand Down
11 changes: 2 additions & 9 deletions scipy/spatial/qhull_misc.h
@@ -1,19 +1,12 @@
/*
* Handle different Fortran conventions and qh_new_qhull_scipy entry point.
* Handle qh_new_qhull_scipy entry point.
*/
#ifndef QHULL_MISC_H_
#define QHULL_MISC_H_

/* for CBLAS_INT only*/
#include "npy_cblas.h"

void BLAS_FUNC(dgetrf)(CBLAS_INT*, CBLAS_INT*, double*, CBLAS_INT*, CBLAS_INT*, CBLAS_INT*);
void BLAS_FUNC(dgecon)(char*, CBLAS_INT*, double*, CBLAS_INT*, double*, double*, double*, CBLAS_INT*, CBLAS_INT*, size_t);
void BLAS_FUNC(dgetrs)(char*, CBLAS_INT*, CBLAS_INT*, double*, CBLAS_INT*, CBLAS_INT*, double*, CBLAS_INT*, CBLAS_INT*, size_t);

#define qh_dgetrf(m,n,a,lda,ipiv,info) BLAS_FUNC(dgetrf)(m,n,a,lda,ipiv,info)
#define qh_dgecon(norm,n,a,lda,anorm,rcond,work,iwork,info) BLAS_FUNC(dgecon)(norm,n,a,lda,anorm,rcond,work,iwork,info,1)
#define qh_dgetrs(trans,n,nrhs,a,lda,ipiv,b,ldb,info) BLAS_FUNC(dgetrs)(trans,n,nrhs,a,lda,ipiv,b,ldb,info,1)

#define qhull_misc_lib_check() QHULL_LIB_CHECK

#include "qhull_src/src/libqhull_r.h"
Expand Down

0 comments on commit abb04b2

Please sign in to comment.