diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/clacon2.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/clacon2.c
index 107bb64..660b42f 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/clacon2.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/clacon2.c
@@ -106,6 +106,11 @@ clacon2_(int *n, complex *v, complex *x, float *est, int *kase, int isave[3])
     extern float smach(char *);
     extern int icmax1_slu(int *, complex *, int *);
     extern double scsum1_slu(int *, complex *, int *);
+#ifdef _CRAY
+    extern int CCOPY(int *, complex *, int *, complex *, int *);
+#else
+    extern int ccopy_(int *, complex *, int *, complex *, int *);
+#endif
 
     safmin = smach("Safe minimum");  /* lamch_("Safe minimum"); */
     if ( *kase == 0 ) {
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/dmach.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/dmach.c
index 73beacb..cafdf1c 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/dmach.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/dmach.c
@@ -11,6 +11,7 @@ at the top-level directory.
 #include <float.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 
 double dmach(char *cmach)
 {
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_cdrop_row.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_cdrop_row.c
index 4987548..89a0f68 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_cdrop_row.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_cdrop_row.c
@@ -28,6 +28,7 @@ extern void caxpy_(int *, complex *, complex [], int *, complex [], int *);
 extern void ccopy_(int *, complex [], int *, complex [], int *);
 extern float scasum_(int *, complex *, int *);
 extern float scnrm2_(int *, complex *, int *);
+extern int scopy_(int *, float *, int *, float *, int *);
 extern double dnrm2_(int *, double [], int *);
 extern int icamax_(int *, complex [], int *);
 
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_zdrop_row.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_zdrop_row.c
index f434dd9..a8dda39 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_zdrop_row.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/ilu_zdrop_row.c
@@ -23,6 +23,7 @@ at the top-level directory.
 #include <stdlib.h>
 #include "slu_zdefs.h"
 
+extern void dcopy_(int *, double [], int *, double [], int *);
 extern void zswap_(int *, doublecomplex [], int *, doublecomplex [], int *);
 extern void zaxpy_(int *, doublecomplex *, doublecomplex [], int *, doublecomplex [], int *);
 extern void zcopy_(int *, doublecomplex [], int *, doublecomplex [], int *);
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/slacon2.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/slacon2.c
index 7c93341..9fcbac9 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/slacon2.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/slacon2.c
@@ -99,10 +99,12 @@ slacon2_(int *n, float *v, float *x, int *isgn, float *est, int *kase, int isave
     int i;
     float temp;
 #ifdef _CRAY
+    extern int IDAMAX(int *, float *, int *);
     extern int ISAMAX(int *, float *, int *);
     extern float SASUM(int *, float *, int *);
     extern int SCOPY(int *, float *, int *, float *, int *);
 #else
+    extern int idamax_(int *, double *, int *);
     extern int isamax_(int *, float *, int *);
     extern float sasum_(int *, float *, int *);
     extern int scopy_(int *, float *, int *, float *, int *);
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/smach.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/smach.c
index fff6c5f..0b69991 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/smach.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/smach.c
@@ -11,6 +11,7 @@ at the top-level directory.
 #include <float.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 
 float smach(char *cmach)
 {
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/sp_ienv.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/sp_ienv.c
index 855d901..c9d25cd 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/sp_ienv.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/sp_ienv.c
@@ -25,6 +25,8 @@ at the top-level directory.
  */
 #include "slu_Cnames.h"
 
+extern int     input_error(char *, int *);
+
 /*! \brief
 
  <pre>
diff --git a/scipy/sparse/linalg/dsolve/SuperLU/SRC/zlacon2.c b/scipy/sparse/linalg/dsolve/SuperLU/SRC/zlacon2.c
index b43c619..ed5f2b7 100644
--- a/scipy/sparse/linalg/dsolve/SuperLU/SRC/zlacon2.c
+++ b/scipy/sparse/linalg/dsolve/SuperLU/SRC/zlacon2.c
@@ -106,6 +106,11 @@ zlacon2_(int *n, doublecomplex *v, doublecomplex *x, double *est, int *kase, int
     extern double dmach(char *);
     extern int izmax1_slu(int *, doublecomplex *, int *);
     extern double dzsum1_slu(int *, doublecomplex *, int *);
+#ifdef _CRAY
+    extern int CCOPY(int *, doublecomplex *, int *, doublecomplex *, int *);
+#else
+    extern int zcopy_(int *, doublecomplex *, int *, doublecomplex *, int *);
+#endif
 
     safmin = dmach("Safe minimum");  /* lamch_("Safe minimum"); */
     if ( *kase == 0 ) {
diff --git a/scipy/cluster/_vq.pyx b/scipy/cluster/_vq.pyx
index 9db0e35..4dd7ef2 100644
--- a/scipy/cluster/_vq.pyx
+++ b/scipy/cluster/_vq.pyx
@@ -12,7 +12,7 @@ from __future__ import absolute_import
 cimport cython
 import numpy as np
 cimport numpy as np
-from .cluster_blas cimport f_dgemm, f_sgemm
+from scipy.linalg.cython_blas cimport dgemm, sgemm
 
 from libc.math cimport sqrt
 
@@ -52,11 +52,11 @@ cdef inline void cal_M(int nobs, int ncodes, int nfeat, vq_type *obs,
     # Call BLAS functions with Fortran ABI
     # Note that BLAS Fortran ABI uses column-major order
     if vq_type is float32_t:
-        f_sgemm("T", "N", &ncodes, &nobs, &nfeat,
-                &alpha, code_book, &nfeat, obs, &nfeat, &beta, M, &ncodes)
+        sgemm("T", "N", &ncodes, &nobs, &nfeat,
+              &alpha, code_book, &nfeat, obs, &nfeat, &beta, M, &ncodes)
     else:
-        f_dgemm("T", "N", &ncodes, &nobs, &nfeat,
-                &alpha, code_book, &nfeat, obs, &nfeat, &beta, M, &ncodes)
+        dgemm("T", "N", &ncodes, &nobs, &nfeat,
+              &alpha, code_book, &nfeat, obs, &nfeat, &beta, M, &ncodes)
 
 
 cdef int _vq(vq_type *obs, vq_type *code_book,
diff --git a/scipy/interpolate/_ppoly.pyx b/scipy/interpolate/_ppoly.pyx
index 308af66..e23a706 100644
--- a/scipy/interpolate/_ppoly.pyx
+++ b/scipy/interpolate/_ppoly.pyx
@@ -14,18 +14,14 @@ cimport cython
 cimport libc.stdlib
 cimport libc.math
 
+from scipy.linalg.cython_lapack cimport dgeev
+
 ctypedef double complex double_complex
 
 ctypedef fused double_or_complex:
     double
     double complex
 
-cdef extern from "blas_defs.h":
-    void c_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)
-
 cdef extern from "numpy/npy_math.h":
     double nan "NPY_NAN"
 
@@ -930,8 +926,8 @@ cdef int croots_poly1(double[:,:,::1] c, double y, int ci, int cj,
 
     # Compute companion matrix eigenvalues
     info = 0
-    c_dgeev("N", "N", &order, a, &order, <double*>wr, <double*>wi,
-            NULL, &order, NULL, &order, work, &lwork, &info)
+    dgeev("N", "N", &order, a, &order, <double*>wr, <double*>wi,
+          NULL, &order, NULL, &order, work, &lwork, &info)
     if info != 0:
         # Failure
         return -2
diff --git a/scipy/spatial/qhull.pyx b/scipy/spatial/qhull.pyx
index cfc5d34..d823ae4 100644
--- a/scipy/spatial/qhull.pyx
+++ b/scipy/spatial/qhull.pyx
@@ -21,6 +21,8 @@ from . cimport setlist
 from libc cimport stdlib
 from scipy._lib.messagestream cimport MessageStream
 
+from scipy.linalg.cython_lapack cimport dgetrf, dgecon, dgetrs
+
 from numpy.compat import asbytes
 import os
 import sys
@@ -206,12 +208,6 @@ from libc.stdlib cimport qsort
 
 cdef extern from "qhull_misc.h":
     void qhull_misc_lib_check()
-    void qh_dgetrf(int *m, int *n, double *a, int *lda, int *ipiv,
-                   int *info) nogil
-    void qh_dgetrs(char *trans, int *n, int *nrhs, double *a, int *lda,
-                   int *ipiv, double *b, int *ldb, int *info) nogil
-    void qh_dgecon(char *norm, int *n, double *a, int *lda, double *anorm,
-                   double *rcond, double *work, int *iwork, int *info) nogil
 
 
 #------------------------------------------------------------------------------
@@ -1116,11 +1112,11 @@ 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,
+                dgecon("1", &n, <double*>T.data, &lda, &anorm, &rcond,
                           work, iwork, &info)
 
                 if rcond < rcond_limit:
@@ -1129,7 +1125,7 @@ def _get_barycentric_transforms(np.ndarray[np.double_t, ndim=2] points,
 
             # 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)