ENH: When cblas is available use it in descr->f->dot.

Importing _dotblas currently executes _dotblas.alterdot, which replaces the default descr->f->dot function with a cblas based version for float, double, complex float, and complex double data types. This PR changes the default descr->f->dot to use cblas whenever it is available. After this change, the alterdot and restoredot functions serve no purpose, so are changed to do nothing and deprecated. Note that those functions were already doing nothing when _dotblas was not available.
numpy · charris · Mar 31, 2014 · Aug 11, 2014 · Aug 14, 2014 · Aug 15, 2014
commit ace49c2ebcd42b9d41bdac11729f3c58e525480e
diff --git a/numpy/core/blasdot/_dotblas.c b/numpy/core/blasdot/_dotblas.c
@@ -24,121 +24,6 @@ static char module_doc[] =
 
 static PyArray_DotFunc *oldFunctions[NPY_NTYPES];
 
-#define MIN(a, b)   ((a) < (b) ? (a) : (b))
-
 /*
- * Convert NumPy stride to BLAS stride. Returns 0 if conversion cannot be done
- * (BLAS won't handle negative or zero strides the way we want).
- */
-static NPY_INLINE int
-blas_stride(npy_intp stride, unsigned itemsize)
-{
-    if (stride <= 0 || stride % itemsize != 0) {
-        return 0;
-    }
-    stride /= itemsize;
-
-    if (stride > INT_MAX) {
-        return 0;
-    }
-    return stride;
-}
-
-/*
- * The following functions do a "chunked" dot product using BLAS when
- * sizeof(npy_intp) > sizeof(int), because BLAS libraries can typically not
- * handle more than INT_MAX elements per call.
- *
- * The chunksize is the greatest power of two less than INT_MAX.
- */
-#if NPY_MAX_INTP > INT_MAX
-# define CHUNKSIZE  (INT_MAX / 2 + 1)
-#else
-# define CHUNKSIZE  NPY_MAX_INTP
-#endif
-
-static void
-FLOAT_dot(void *a, npy_intp stridea, void *b, npy_intp strideb, void *res,
-          npy_intp n, void *tmp)
-{
-    int na = blas_stride(stridea, sizeof(float));
-    int nb = blas_stride(strideb, sizeof(float));
-
-    if (na && nb) {
-        double r = 0.;          /* double for stability */
-        float *fa = a, *fb = b;
-
-        while (n > 0) {
-            int chunk = MIN(n, CHUNKSIZE);
-
-            r += cblas_sdot(chunk, fa, na, fb, nb);
-            fa += chunk * na;
-            fb += chunk * nb;
-            n -= chunk;
-        }
-        *((float *)res) = r;
-    }
-    else {
-        oldFunctions[NPY_FLOAT](a, stridea, b, strideb, res, n, tmp);
-    }
-}
-
-static void
-DOUBLE_dot(void *a, npy_intp stridea, void *b, npy_intp strideb, void *res,
-           npy_intp n, void *tmp)
-{
-    int na = blas_stride(stridea, sizeof(double));
-    int nb = blas_stride(strideb, sizeof(double));
-
-    if (na && nb) {
-        double r = 0.;
-        double *da = a, *db = b;
-
-        while (n > 0) {
-            int chunk = MIN(n, CHUNKSIZE);
-
-            r += cblas_ddot(chunk, da, na, db, nb);
-            da += chunk * na;
-            db += chunk * nb;
-            n -= chunk;
-        }
-        *((double *)res) = r;
-    }
-    else {
-        oldFunctions[NPY_DOUBLE](a, stridea, b, strideb, res, n, tmp);
-    }
-}
-
-static void
-CFLOAT_dot(void *a, npy_intp stridea, void *b, npy_intp strideb, void *res,
-           npy_intp n, void *tmp)
-{
-    int na = blas_stride(stridea, sizeof(npy_cfloat));
-    int nb = blas_stride(strideb, sizeof(npy_cfloat));
-
-    if (na && nb) {
-        cblas_cdotu_sub((int)n, (float *)a, na, (float *)b, nb, (float *)res);
-    }
-    else {
-        oldFunctions[NPY_CFLOAT](a, stridea, b, strideb, res, n, tmp);
-    }
-}
-
-static void
-CDOUBLE_dot(void *a, npy_intp stridea, void *b, npy_intp strideb, void *res,
-            npy_intp n, void *tmp)
-{
-    int na = blas_stride(stridea, sizeof(npy_cdouble));
-    int nb = blas_stride(strideb, sizeof(npy_cdouble));
-
-    if (na && nb) {
-        cblas_zdotu_sub((int)n, (double *)a, na, (double *)b, nb,
-                        (double *)res);
-    }
-    else {
-        oldFunctions[NPY_CDOUBLE](a, stridea, b, strideb, res, n, tmp);
-    }
-}
 
 /*
  * Helper: call appropriate BLAS dot function for typenum.
@@ -149,20 +34,8 @@ blas_dot(int typenum, npy_intp n,
          void *a, npy_intp stridea, void *b, npy_intp strideb, void *res)
 {
     PyArray_DotFunc *dot = NULL;
-    switch (typenum) {
-        case NPY_DOUBLE:
-            dot = DOUBLE_dot;
-            break;
-        case NPY_FLOAT:
-            dot = FLOAT_dot;
-            break;
-        case NPY_CDOUBLE:
-            dot = CDOUBLE_dot;
-            break;
-        case NPY_CFLOAT:
-            dot = CFLOAT_dot;
-            break;
-    }
+
+    dot = oldFunctions[typenum];
     assert(dot != NULL);
     dot(a, stridea, b, strideb, res, n, NULL);
 }
@@ -257,19 +130,15 @@ dotblas_alterdot(PyObject *NPY_UNUSED(dummy), PyObject *args)
     if (!altered) {
         descr = PyArray_DescrFromType(NPY_FLOAT);
         oldFunctions[NPY_FLOAT] = descr->f->dotfunc;
-        descr->f->dotfunc = (PyArray_DotFunc *)FLOAT_dot;
 
         descr = PyArray_DescrFromType(NPY_DOUBLE);
         oldFunctions[NPY_DOUBLE] = descr->f->dotfunc;
-        descr->f->dotfunc = (PyArray_DotFunc *)DOUBLE_dot;
 
         descr = PyArray_DescrFromType(NPY_CFLOAT);
         oldFunctions[NPY_CFLOAT] = descr->f->dotfunc;
-        descr->f->dotfunc = (PyArray_DotFunc *)CFLOAT_dot;
 
         descr = PyArray_DescrFromType(NPY_CDOUBLE);
         oldFunctions[NPY_CDOUBLE] = descr->f->dotfunc;
-        descr->f->dotfunc = (PyArray_DotFunc *)CDOUBLE_dot;
 
         altered = NPY_TRUE;
     }

diff --git a/numpy/core/bscript b/numpy/core/bscript
@@ -488,14 +488,22 @@ def pre_build(context):
                 pjoin('src', 'multiarray', 'usertypes.c')]
         else:
             sources = extension.sources
+
+        use = 'npysort npymath'
+        defines = ['_FILE_OFFSET_BITS=64',
+                   '_LARGEFILE_SOURCE=1',
+                   '_LARGEFILE64_SOURCE=1']
+
+        if bld.env.HAS_CBLAS:
+            use += ' CBLAS'
+            defines.append('HAVE_CBLAS')
+
         includes = ["src/multiarray", "src/private"]
         return context.default_builder(extension,
                                        includes=includes,
                                        source=sources,
-                                       use="npysort npymath",
-                                       defines=['_FILE_OFFSET_BITS=64',
-                                                '_LARGEFILE_SOURCE=1',
-                                                '_LARGEFILE64_SOURCE=1']
+                                       use=use,
+                                       defines=defines
                                        )
     context.register_builder("multiarray", builder_multiarray)
 

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
@@ -1078,31 +1078,33 @@ def outer(a, b, out=None):
 # try to import blas optimized dot if available
 envbak = os.environ.copy()
 try:
-    # importing this changes the dot function for basic 4 types
-    # to blas-optimized versions.
-
     # disables openblas affinity setting of the main thread that limits
     # python threads or processes to one core
     if 'OPENBLAS_MAIN_FREE' not in os.environ:
         os.environ['OPENBLAS_MAIN_FREE'] = '1'
     if 'GOTOBLAS_MAIN_FREE' not in os.environ:
         os.environ['GOTOBLAS_MAIN_FREE'] = '1'
-    from ._dotblas import dot, vdot, inner, alterdot, restoredot
+    from ._dotblas import dot, vdot, inner
 except ImportError:
     # docstrings are in add_newdocs.py
     inner = multiarray.inner
     dot = multiarray.dot
     def vdot(a, b):
         return dot(asarray(a).ravel().conj(), asarray(b).ravel())
-    def alterdot():
-            pass
-    def restoredot():
-        pass
 finally:
     os.environ.clear()
     os.environ.update(envbak)
     del envbak
 
+
+def alterdot():
+    warnings.warn("alterdot no longer does anything.", DeprecationWarning)
+
+
+def restoredot():
+    warnings.warn("restoredot no longer does anything.", DeprecationWarning)
+
+
 def tensordot(a, b, axes=2):
     """
     Compute tensor dot product along specified axes for arrays >= 1-D.

diff --git a/numpy/core/setup.py b/numpy/core/setup.py
@@ -846,15 +846,22 @@ def generate_multiarray_templated_sources(ext, build_dir):
         multiarray_src = [join('src', 'multiarray', 'multiarraymodule_onefile.c')]
         multiarray_src.append(generate_multiarray_templated_sources)
 
+    blas_info = get_info('blas_opt', 0)
+    if blas_info and  ('HAVE_CBLAS', None) in blas_info.get('define_macros', []):
+        extra_info = blas_info
+    else:
+        extra_info = {}
+
     config.add_extension('multiarray',
-                         sources = multiarray_src +
+                         sources=multiarray_src +
                                  [generate_config_h,
-                                 generate_numpyconfig_h,
-                                 generate_numpy_api,
-                                 join(codegen_dir, 'generate_numpy_api.py'),
-                                 join('*.py')],
-                         depends = deps + multiarray_deps,
-                         libraries = ['npymath', 'npysort'])
+                                  generate_numpyconfig_h,
+                                  generate_numpy_api,
+                                  join(codegen_dir, 'generate_numpy_api.py'),
+                                  join('*.py')],
+                         depends=deps + multiarray_deps,
+                         libraries=['npymath', 'npysort'],
+                         extra_info=extra_info)
 
     #######################################################################
     #                           umath module                              #