Repository: scipy/scipy
Pull request: MAINT:linalg: Expose Cython functions for generic use
Closes #16143 

Certain low-level Cython functions are useful for multiple functions, such as checking for triangular-ness, or symmetry/Hermitian structure. 

With generous help of @ev-br hopefully this PR manages the cimportable cdef'd functions. 

@rgommers I think this needs you attention since we butchered our way out for the Meson stuff.

Problem statement:
ENH: Expose Cythonized functions to other pyx files through cimport
This is more of a learning issue that I have been trying out for a few weeks now but couldn't have any progress. What I would like to have is 
 - expose a Cython API for linalg that is similar to `cython_blas`, `cython_lapack` or `cython_special`. 
 - or at least make certain Cythonized functions common to all scipy functions. 

I think I can say that I really don't understand the machinery `_generate_pyx.py` or the `custom_target` of meson yet. 

#### A particular case study: 

In `matfuncs_expm.pyx.in` file there is a function that is a fast transpose for memoryviews. What I wanted to have is to make it available through `_cythonized_linalg_utils.pyx` and `cimport` it via other pyx files. To do that it turns out that I need a pxd file for the definitions since I am constantly getting 

```
from ._cythonized_array_utils cimport swap_c_and_f_layout
^
------------------------------------------------------------

scipy\linalg\_matfuncs_expm.pyx:31:0: 'scipy\linalg\_cythonized_array_utils.pxd' not found
```

But processing pxd file is quite the black magic both in distutils and also meson. Or put better, I really don't understand it at all. And I would like to throw the towel after wasting a lot of time. 

The resulting structure I tried to have is 

```
-- __cythonized_array_utils  #  this namespace holds the cimportable 'swap_c_f_layout' function
 | ...
 | ...
 |- _matfuncs_expm.pyx.in #  this file cimports the 'swap_c_f_layout'
 |- ...
```

The optimal outcome is to expose a Cython `cython_linalg_utils` API however that turned out to be even more complicated hence the case study above is already enough for me. Maybe it can be also beneficial for others in the future to expose a Cython API in other modules.

#### What I have tried

There are many permutations I've tried but the generic attempt was to create a pxd file and change `setup.py` to

```python
    config.add_data_files('_cythonized_array_utils.pxd')
    config.add_extension('_cythonized_array_utils',
                         sources=['_cythonized_array_utils.c',
                                        '_cythonized_array_utils.pxd'])  #  also tried the "depends" etc. but
                                                                         #  documentation is awful for distutils on these
```

Hints from issue discussion:

Commit message: MAINT:linalg: Expose low-level Cython functions for generic use

Patch to review:
diff --git a/scipy/linalg/_cythonized_array_utils.pxd b/scipy/linalg/_cythonized_array_utils.pxd
new file mode 100644
index 000000000000..6263fa249272
--- /dev/null
+++ b/scipy/linalg/_cythonized_array_utils.pxd
@@ -0,0 +1,40 @@
+cimport numpy as cnp
+
+ctypedef fused lapack_t:
+    float
+    double
+    (float complex)
+    (double complex)
+
+ctypedef fused lapack_cz_t:
+    (float complex)
+    (double complex)
+
+ctypedef fused lapack_sd_t:
+    float
+    double
+
+ctypedef fused np_numeric_t:
+    cnp.int8_t
+    cnp.int16_t
+    cnp.int32_t
+    cnp.int64_t
+    cnp.uint8_t
+    cnp.uint16_t
+    cnp.uint32_t
+    cnp.uint64_t
+    cnp.float32_t
+    cnp.float64_t
+    cnp.longdouble_t
+    cnp.complex64_t
+    cnp.complex128_t
+
+ctypedef fused np_complex_numeric_t:
+    cnp.complex64_t
+    cnp.complex128_t
+
+
+cdef void swap_c_and_f_layout(lapack_t *a, lapack_t *b, int r, int c, int n) nogil
+cdef (int, int) band_check_internal_c(np_numeric_t[:, ::1]A) nogil
+cdef bint is_sym_her_real_c_internal(np_numeric_t[:, ::1]A) nogil
+cdef bint is_sym_her_complex_c_internal(np_complex_numeric_t[:, ::1]A) nogil
diff --git a/scipy/linalg/_cythonized_array_utils.pyx b/scipy/linalg/_cythonized_array_utils.pyx
index 1a074b19849a..c9beee0b3df7 100644
--- a/scipy/linalg/_cythonized_array_utils.pyx
+++ b/scipy/linalg/_cythonized_array_utils.pyx
@@ -2,29 +2,47 @@
 cimport cython
 cimport numpy as cnp
 import numpy as np
-
-ctypedef fused np_numeric_t:
-    cnp.int8_t
-    cnp.int16_t
-    cnp.int32_t
-    cnp.int64_t
-    cnp.uint8_t
-    cnp.uint16_t
-    cnp.uint32_t
-    cnp.uint64_t
-    cnp.float32_t
-    cnp.float64_t
-    cnp.longdouble_t
-    cnp.complex64_t
-    cnp.complex128_t
-
-ctypedef fused np_complex_numeric_t:
-    cnp.complex64_t
-    cnp.complex128_t
+from scipy.linalg._cythonized_array_utils cimport (
+    lapack_t,
+	np_complex_numeric_t,
+	np_numeric_t
+    )
 
 __all__ = ['bandwidth', 'issymmetric', 'ishermitian']
 
 
+# ====================== swap_c_and_f_layout : s, d, c, z ====================
+@cython.cdivision(True)
+@cython.wraparound(False)
+@cython.boundscheck(False)
+@cython.initializedcheck(False)
+cdef void swap_c_and_f_layout(lapack_t *a, lapack_t *b, int r, int c, int n) nogil:
+    """Recursive matrix transposition for square arrays"""
+    cdef int i, j, ith_row, r2, c2
+    cdef lapack_t *bb=b
+    cdef lapack_t *aa=a
+    if r < 16:
+        for j in range(c):
+            ith_row = 0
+            for i in range(r):
+            # Basically b[i*n+j] = a[j*n+i] without index math
+                bb[ith_row] = aa[i]
+                ith_row += n
+            aa += n
+            bb += 1
+    else:
+        # If tall
+        if (r > c):
+            r2 = r//2
+            swap_c_and_f_layout(a, b, r2, c, n)
+            swap_c_and_f_layout(a + r2, b+(r2)*n, r-r2, c, n)
+        else:  # Nope
+            c2 = c//2
+            swap_c_and_f_layout(a, b, r, c2, n);
+            swap_c_and_f_layout(a+(c2)*n, b+c2, r, c-c2, n)
+# ============================================================================
+
+
 @cython.embedsignature(True)
 def bandwidth(a):
     """Return the lower and upper bandwidth of a 2D numeric array.
diff --git a/scipy/linalg/_matfuncs_expm.pyx.in b/scipy/linalg/_matfuncs_expm.pyx.in
index 05f712ac93b7..41b42c87c6a1 100644
--- a/scipy/linalg/_matfuncs_expm.pyx.in
+++ b/scipy/linalg/_matfuncs_expm.pyx.in
@@ -14,24 +14,14 @@ from scipy.linalg.cython_blas cimport (sgemm, saxpy, sscal, scopy, sgemv,
                                        zgemm, zaxpy, zscal, zcopy, dgemv,
                                        csscal, zdscal, idamax)
 
+from ._cythonized_array_utils cimport (lapack_t, lapack_cz_t,
+        lapack_sd_t, swap_c_and_f_layout)
+
 __all__ = ['pick_pade_structure', 'pade_UV_calc']
 
 # See GH-14813
 cnp.import_array()
 
-ctypedef fused lapack_t:
-    float
-    double
-    (float complex)
-    (double complex)
-
-ctypedef fused lapack_cz_t:
-    (float complex)
-    (double complex)
-
-ctypedef fused lapack_sd_t:
-    float
-    double
 
 {{py:
 # Tempita variables
@@ -39,36 +29,7 @@ typenames = ['float', 'double', 'float complex', 'double complex']
 precision = ['float', 'double', 'float', 'double']
 prefixes = ['s', 'd', 'c', 'z']
 }}
-# ====================== swap_c_and_f_layout : s, d, c, z ====================
-@cython.cdivision(True)
-@cython.wraparound(False)
-@cython.boundscheck(False)
-@cython.initializedcheck(False)
-cdef void swap_c_and_f_layout(lapack_t *a, lapack_t *b, int r, int c, int n) nogil:
-    """Recursive matrix transposition for square arrays"""
-    cdef int i, j, ith_row, r2, c2
-    cdef lapack_t *bb=b
-    cdef lapack_t *aa=a
-    if r < 16:
-        for j in range(c):
-            ith_row = 0
-            for i in range(r):
-            # Basically b[i*n+j] = a[j*n+i] without index math
-                bb[ith_row] = aa[i]
-                ith_row += n
-            aa += n
-            bb += 1
-    else:
-        # If tall
-        if (r > c):
-            r2 = r//2
-            swap_c_and_f_layout(a, b, r2, c, n)
-            swap_c_and_f_layout(a + r2, b+(r2)*n, r-r2, c, n)
-        else:  # Nope
-            c2 = c//2
-            swap_c_and_f_layout(a, b, r, c2, n);
-            swap_c_and_f_layout(a+(c2)*n, b+c2, r, c-c2, n)
-# ============================================================================
+
 
 # ========================= norm1 : s, d, c, z ===============================
 # The generic abs() function is made nogil friendly in Cython 3.x (unreleased)
diff --git a/scipy/linalg/meson.build b/scipy/linalg/meson.build
index 2426523e8ed7..1271df281cc1 100644
--- a/scipy/linalg/meson.build
+++ b/scipy/linalg/meson.build
@@ -222,6 +222,30 @@ cython_lapack = py3.extension_module('cython_lapack',
   subdir: 'scipy/linalg'
 )
 
+
+py3.install_sources(
+  ['_cythonized_array_utils.pyi', '_cythonized_array_utils.pxd'],
+  pure: false,
+  subdir: 'scipy/linalg',
+)
+
+# This is needed to add pxd files into build dir ala "cp"
+# See 
+configure_file(input : '_cythonized_array_utils.pxd',
+               output : '_cythonized_array_utils.pxd',
+               copy : true)
+
+
+_cythonized_array_utils = py3.extension_module('_cythonized_array_utils',
+  linalg_init_cython_gen.process('_cythonized_array_utils.pyx'),
+  c_args: [cython_c_args, c_undefined_ok],
+  dependencies: [py3_dep],
+  include_directories: [inc_np],
+  install: true,
+  subdir: 'scipy/linalg'
+)
+
+
 _decomp_update_pyx = custom_target('_decomp_update',
   output: '_decomp_update.pyx',
   input: '_decomp_update.pyx.in',
@@ -253,14 +277,7 @@ py3.extension_module('_matfuncs_expm',
   subdir: 'scipy/linalg'
 )
 
-_cythonized_array_utils = py3.extension_module('_cythonized_array_utils',
-  linalg_init_cython_gen.process('_cythonized_array_utils.pyx'),
-  c_args: [cython_c_args, c_undefined_ok],
-  dependencies: [py3_dep],
-  include_directories: [inc_np],
-  install: true,
-  subdir: 'scipy/linalg'
-)
+
 
 python_sources = [
   '__init__.py',
@@ -310,11 +327,7 @@ py3.install_sources(
   subdir: 'scipy/linalg'
 )
 
-py3.install_sources(
-  ['_cythonized_array_utils.pyi'],
-  pure: false,
-  subdir: 'scipy/linalg'
-)
+
 
 # TODO: install (only) these pxd files. Neither py3.install_sources nor
 #       install_data accept generated data files. What
diff --git a/scipy/linalg/setup.py b/scipy/linalg/setup.py
index ec012c6f1c4d..199b70267f01 100644
--- a/scipy/linalg/setup.py
+++ b/scipy/linalg/setup.py
@@ -146,8 +146,14 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_decomp_update',
                          sources=['_decomp_update.c'])
 
+    config.add_data_files('_cythonized_array_utils.pxd')
+    
     config.add_extension('_cythonized_array_utils',
-                         sources=['_cythonized_array_utils.c'])
+                         sources=['_cythonized_array_utils.c'],
+                         depends=['_cythonized_array_utils.pyx',
+                                  '_cythonized_array_utils.pxd'],
+                                  include_dirs=['.']
+                                  )
 
     config.add_extension('_matfuncs_expm',
                          sources=['_matfuncs_expm.c'])


Code review comments:
diff_hunk: @@ -222,6 +222,30 @@ cython_lapack = py3.extension_module('cython_lapack',
   subdir: 'scipy/linalg'
 )
 
+
+py3.install_sources(
+  ['_cythonized_array_utils.pyi', '_cythonized_array_utils.pxd'],
+  pure: false,
+  subdir: 'scipy/linalg',
+)
+
+# This is needed to add pxd files into build dir ala "cp"
+# See 
path: scipy/linalg/meson.build
text: @author:
bah I forgot to add the actual link to `copy : true` warning creation

diff_hunk: @@ -146,8 +146,14 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_decomp_update',
                          sources=['_decomp_update.c'])
 
+    config.add_data_files('_cythonized_array_utils.pxd')
+    
     config.add_extension('_cythonized_array_utils',
-                         sources=['_cythonized_array_utils.c'])
+                         sources=['_cythonized_array_utils.c'],
+                         depends=['_cythonized_array_utils.pyx',
+                                  '_cythonized_array_utils.pxd'],
+                                  include_dirs=['.']
path: scipy/linalg/setup.py
text: @author:
I am not sure if this line is needed. I got varying results and I certainly don't understand anything that is going on here. But this started to give consistent results.

@user1:
Yes indeed, it's needed with `distutils`. For Meson it never is; the current source directory is automatically added.