ENH: Adopt new macOS Accelerate BLAS/LAPACK Interfaces, including ILP64

[Developer-Ecosystem-Engineering]

macOS 13.3 shipped with an updated Accelerate framework that provides BLAS / LAPACK. The new version is aligned with Netlib's v3.9.1 and also supports ILP64. The changes here adopt those new interfaces when available.

• New interfaces are used when ACCELERATE_NEW_LAPACK is defined.
• ILP64 interfaces are used when both ACCELERATE_NEW_LAPACK and ACCELERATE_LAPACK_ILP64 are defined.

macOS 13.3 now ships with 3 different sets of BLAS / LAPACK interfaces:

• LP64 / LAPACK v3.2.1 - legacy interfaces kept for compatibility
• LP64 / LAPACK v3.9.1 - new interfaces
• ILP64 / LAPACK v3.9.1 - new interfaces with ILP64 support

For LP64, we want to support building against macOS 13.3+ SDK, but having it work on pre-13.3 systems. To that end, we create wrappers for each API that do a runtime check on which set of API is available and should be used.

ILP64 is only supported on macOS 13.3+ and does not use additional wrappers.

We've included support for both distutils and Meson builds.

All tests pass on Apple silicon and Intel based Macs.

Benchmarks
ILP64 Accelerate vs OpenBLAS

   before           after         ratio
 [73f0cf4f]       [d1572653]
 <openblas-ilp64>       <accelerate-ilp64>
          n/a              n/a      n/a  bench_linalg.Linalg.time_op('det', 'float16')
          n/a              n/a      n/a  bench_linalg.Linalg.time_op('pinv', 'float16')
          n/a              n/a      n/a  bench_linalg.Linalg.time_op('svd', 'float16')
       failed           failed      n/a  bench_linalg.LinalgSmallArrays.time_det_small_array
  3.96±0.1μs       5.04±0.4μs     1.27  bench_linalg.Linalg.time_op('norm', 'float32')
  1.43±0.04ms         1.43±0ms     1.00  bench_linalg.Einsum.time_einsum_outer(<class 'numpy.float32'>)
   12.7±0.4μs       12.7±0.3μs     1.00  bench_linalg.Einsum.time_einsum_sum_mul2(<class 'numpy.float32'>)
   24.1±0.8μs      24.1±0.04μs     1.00  bench_linalg.Linalg.time_op('norm', 'float16')
   9.48±0.2ms       9.48±0.3ms     1.00  bench_linalg.Einsum.time_einsum_outer(<class 'numpy.float64'>)
     609±20μs          609±2μs     1.00  bench_linalg.Einsum.time_einsum_noncon_outer(<class 'numpy.float32'>)
     64.9±2μs      64.7±0.07μs     1.00  bench_linalg.Einsum.time_einsum_contig_outstride0(<class 'numpy.float64'>)
  1.24±0.03ms      1.24±0.01ms     1.00  bench_linalg.Einsum.time_einsum_noncon_outer(<class 'numpy.float64'>)
      102±3μs        102±0.2μs     1.00  bench_linalg.Einsum.time_einsum_contig_contig(<class 'numpy.float64'>)
   21.9±0.8μs      21.8±0.02μs     1.00  bench_linalg.Einsum.time_einsum_multiply(<class 'numpy.float64'>)
   22.8±0.2ms       22.7±0.3ms     0.99  bench_linalg.Eindot.time_einsum_ijk_jil_kl
   13.3±0.4μs      13.3±0.02μs     0.99  bench_linalg.Einsum.time_einsum_sum_mul2(<class 'numpy.float64'>)
   9.56±0.3μs       9.49±0.2μs     0.99  bench_linalg.Einsum.time_einsum_noncon_contig_contig(<class 'numpy.float64'>)
   7.31±0.2μs      7.26±0.08μs     0.99  bench_linalg.Einsum.time_einsum_noncon_contig_outstride0(<class 'numpy.float32'>)
   5.60±0.2ms      5.55±0.02ms     0.99  bench_linalg.Eindot.time_einsum_ij_jk_a_b
     37.1±1μs       36.7±0.1μs     0.99  bench_linalg.Einsum.time_einsum_contig_outstride0(<class 'numpy.float32'>)
   13.5±0.4μs      13.4±0.05μs     0.99  bench_linalg.Einsum.time_einsum_sum_mul(<class 'numpy.float64'>)
  1.03±0.03μs         1.02±0μs     0.99  bench_linalg.LinalgSmallArrays.time_norm_small_array
     51.6±2μs      51.0±0.09μs     0.99  bench_linalg.Einsum.time_einsum_contig_contig(<class 'numpy.float32'>)
   15.2±0.5μs      15.0±0.04μs     0.99  bench_linalg.Einsum.time_einsum_noncon_sum_mul2(<class 'numpy.float64'>)
   13.9±0.4μs      13.7±0.02μs     0.99  bench_linalg.Einsum.time_einsum_noncon_sum_mul2(<class 'numpy.float32'>)
     415±10μs        409±0.4μs     0.99  bench_linalg.Eindot.time_einsum_i_ij_j
   9.29±0.3μs      9.01±0.03μs     0.97  bench_linalg.Einsum.time_einsum_noncon_mul(<class 'numpy.float64'>)
   18.2±0.6μs      17.6±0.04μs     0.97  bench_linalg.Einsum.time_einsum_multiply(<class 'numpy.float32'>)
     509±40μs         492±10μs     0.97  bench_linalg.Einsum.time_einsum_mul(<class 'numpy.float64'>)
   9.63±0.3μs      9.28±0.09μs     0.96  bench_linalg.Einsum.time_einsum_noncon_contig_contig(<class 'numpy.float32'>)
   9.08±0.2μs      8.73±0.02μs     0.96  bench_linalg.Einsum.time_einsum_noncon_mul(<class 'numpy.float32'>)
   15.6±0.5μs      15.0±0.04μs     0.96  bench_linalg.Einsum.time_einsum_noncon_sum_mul(<class 'numpy.float64'>)
   7.74±0.2μs      7.39±0.04μs     0.95  bench_linalg.Einsum.time_einsum_noncon_contig_outstride0(<class 'numpy.float64'>)
   18.6±0.6μs      17.7±0.03μs     0.95  bench_linalg.Einsum.time_einsum_noncon_multiply(<class 'numpy.float32'>)
   14.5±0.4μs      13.7±0.03μs     0.95  bench_linalg.Einsum.time_einsum_noncon_sum_mul(<class 'numpy.float32'>)
   13.3±0.6μs       12.5±0.3μs     0.94  bench_linalg.Einsum.time_einsum_sum_mul(<class 'numpy.float32'>)
   23.5±0.5μs      21.9±0.05μs     0.93  bench_linalg.Einsum.time_einsum_noncon_multiply(<class 'numpy.float64'>)
     264±20μs          243±4μs     0.92  bench_linalg.Einsum.time_einsum_mul(<class 'numpy.float32'>)
    177±50μs        132±0.6μs     0.75  bench_linalg.Eindot.time_dot_trans_at_a
  10.7±0.3μs      7.13±0.01μs     0.67  bench_linalg.Linalg.time_op('norm', 'int16')
   97.5±2μs       64.7±0.1μs     0.66  bench_linalg.Eindot.time_matmul_trans_a_at
  8.87±0.3μs         5.76±0μs     0.65  bench_linalg.Linalg.time_op('norm', 'longfloat')
  8.90±0.3μs      5.77±0.01μs     0.65  bench_linalg.Linalg.time_op('norm', 'float64')
  8.48±0.3μs      5.40±0.01μs     0.64  bench_linalg.Linalg.time_op('norm', 'int64')
    106±2μs         66.5±8μs     0.63  bench_linalg.Eindot.time_inner_trans_a_a
 8.25±0.3μs         5.16±0μs     0.62  bench_linalg.Linalg.time_op('norm', 'int32')
     103±5ms       64.6±0.5ms     0.62  bench_import.Import.time_linalg
     106±3μs       66.0±0.1μs     0.62  bench_linalg.Eindot.time_dot_trans_a_at
    202±20μs        124±0.6μs     0.61  bench_linalg.Eindot.time_matmul_trans_at_a
   31.5±10μs      19.3±0.02μs     0.61  bench_linalg.Eindot.time_dot_d_dot_b_c
   32.4±20μs      19.7±0.03μs     0.61  bench_linalg.Eindot.time_matmul_d_matmul_b_c
    5.05±1ms      3.06±0.09ms     0.61  bench_linalg.Linalg.time_op('svd', 'complex128')
  5.35±0.9ms      3.09±0.09ms     0.58  bench_linalg.Linalg.time_op('svd', 'complex64')
    6.37±3ms       3.27±0.1ms     0.51  bench_linalg.Linalg.time_op('pinv', 'complex128')
    7.26±8ms       3.24±0.1ms     0.45  bench_linalg.Linalg.time_op('pinv', 'complex64')
   519±100μs        219±0.8μs     0.42  bench_linalg.Linalg.time_op('det', 'complex64')
  31.3±0.9μs       12.8±0.1μs     0.41  bench_linalg.Linalg.time_op('norm', 'complex128')
  2.44±0.7ms          924±1μs     0.38  bench_linalg.Linalg.time_op('pinv', 'float64')
  29.9±0.8μs      10.8±0.01μs     0.36  bench_linalg.Linalg.time_op('norm', 'complex64')
  2.56±0.5ms          924±1μs     0.36  bench_linalg.Linalg.time_op('pinv', 'float32')
  2.63±0.5ms        924±0.6μs     0.35  bench_linalg.Linalg.time_op('pinv', 'int64')
  2.68±0.7ms         927±10μs     0.35  bench_linalg.Linalg.time_op('pinv', 'int32')
  2.68±0.5ms         927±10μs     0.35  bench_linalg.Linalg.time_op('pinv', 'int16')
  2.93±0.6ms          925±2μs     0.32  bench_linalg.Linalg.time_op('pinv', 'longfloat')
   809±500μs        215±0.2μs     0.27  bench_linalg.Linalg.time_op('det', 'complex128')
  3.67±0.9ms         895±20μs     0.24  bench_linalg.Eindot.time_tensordot_a_b_axes_1_0_0_1
   489±100μs         114±20μs     0.23  bench_linalg.Eindot.time_inner_trans_a_ac
  3.64±0.7ms        777±0.3μs     0.21  bench_linalg.Lstsq.time_numpy_linalg_lstsq_a__b_float64
    755±90μs         157±10μs     0.21  bench_linalg.Eindot.time_dot_a_b
    4.63±1ms          899±9μs     0.19  bench_linalg.Linalg.time_op('svd', 'longfloat')
    5.19±1ms         922±10μs     0.18  bench_linalg.Linalg.time_op('svd', 'float64')
   599±200μs         89.4±2μs     0.15  bench_linalg.Eindot.time_matmul_trans_atc_a
   956±200μs         140±10μs     0.15  bench_linalg.Eindot.time_matmul_a_b
    6.45±3ms         903±10μs     0.14  bench_linalg.Linalg.time_op('svd', 'float32')
    6.42±3ms        896±0.7μs     0.14  bench_linalg.Linalg.time_op('svd', 'int32')
    6.47±4ms          902±5μs     0.14  bench_linalg.Linalg.time_op('svd', 'int64')
    6.52±1ms          899±2μs     0.14  bench_linalg.Linalg.time_op('svd', 'int16')
   799±300μs          109±2μs     0.14  bench_linalg.Eindot.time_dot_trans_atc_a
   502±100μs       65.0±0.2μs     0.13  bench_linalg.Eindot.time_dot_trans_a_atc
   542±300μs      64.2±0.05μs     0.12  bench_linalg.Eindot.time_matmul_trans_a_atc
   458±300μs      41.6±0.09μs     0.09  bench_linalg.Linalg.time_op('det', 'int32')
   471±100μs      41.9±0.03μs     0.09  bench_linalg.Linalg.time_op('det', 'float32')
   510±100μs      43.6±0.06μs     0.09  bench_linalg.Linalg.time_op('det', 'int16')
   478±200μs      39.6±0.05μs     0.08  bench_linalg.Linalg.time_op('det', 'longfloat')
   599±200μs      39.6±0.09μs     0.07  bench_linalg.Linalg.time_op('det', 'float64')
   758±300μs       41.6±0.1μs     0.05  bench_linalg.Linalg.time_op('det', 'int64')

[Developer-Ecosystem-Engineering]

Reviewing failures

[rgommers]

Thanks @Developer-Ecosystem-Engineering, support for the new Accelerate is something that macOS users will get excited about I think. The benchmarks look great.

There's one CI failure for the Pyodide/Emscripten build which is real:

wasm-ld: error: function signature mismatch: dgeqrf_
>>> defined as (i32, i32, i32, i32, i32, i32, i32, i32) -> void in build/temp.emscripten_3_1_32_wasm32-3.11/numpy/linalg/umath_linalg.o
>>> defined as (i32, i32, i32, i32, i32, i32, i32, i32) -> i32 in build/temp.emscripten_3_1_32_wasm32-3.11/numpy/linalg/lapack_lite/f2c_d_lapack.o

due to this PR changing signatures from int to void (in this case: extern "C" fortran_int FNAME(dgeqrf) to void BLAS_FUNC(dgeqrf)). Was that int to void change actually needed? Netlib uses lapack_int, and I see that clapack.h in the Accelerate framework uses plain int.

[seberg]

I am probably just ignorant, but can't we reuse the existing env variable rather than adding ACCELERATE_LAPACK_ILP64? Also, do we really need to include thousands of lines of header files rather than making it in a sense "just another" blas/lapack?

[Developer-Ecosystem-Engineering]

Thanks @Developer-Ecosystem-Engineering, support for the new Accelerate is something that macOS users will get excited about I think. The benchmarks look great.

There's one CI failure for the Pyodide/Emscripten build which is real:

wasm-ld: error: function signature mismatch: dgeqrf_
>>> defined as (i32, i32, i32, i32, i32, i32, i32, i32) -> void in build/temp.emscripten_3_1_32_wasm32-3.11/numpy/linalg/umath_linalg.o
>>> defined as (i32, i32, i32, i32, i32, i32, i32, i32) -> i32 in build/temp.emscripten_3_1_32_wasm32-3.11/numpy/linalg/lapack_lite/f2c_d_lapack.o

due to this PR changing signatures from int to void (in this case: extern "C" fortran_int FNAME(dgeqrf) to void BLAS_FUNC(dgeqrf)). Was that int to void change actually needed? Netlib uses lapack_int, and I see that clapack.h in the Accelerate framework uses plain int.

Yeah, we swapped it to fortran_int which seemed more appropriate

edit: oh boy, I've made it worse! Let's discuss the env variable bit and come back with another pass to ensure we are on the same page.

[rgommers]

I am probably just ignorant, but can't we reuse the existing env variable rather than adding ACCELERATE_LAPACK_ILP64?

It's not an environment variable, it's a define that is not used in our code base but required for Accelerate to expose the new interface. From https://developer.apple.com/documentation/accelerate/blas:

"Apple provides the BLAS and LAPACK libraries under the Accelerate framework to be in line with LAPACK 3.9.1. These new interfaces provide additional functionality, as well as a new ILP64 interface. To use the new interfaces, define ACCELERATE_NEW_LAPACK before including the Accelerate or vecLib headers. For ILP64 interfaces, also define ACCELERATE_LAPACK_ILP64. For Swift projects, specify ACCELERATE_NEW_LAPACK=1 and ACCELERATE_LAPACK_ILP64=1 as preprocessor macros in Xcode build settings."

Also, do we really need to include thousands of lines of header files rather than making it in a sense "just another" blas/lapack?

The problem is that these interfaces are new in macOS 13.3, so if we treat it like a separate library that means we'd have to double the amount of wheels we ship for macOS - one set for <13.3 and one for >=13.3. I'd much rather carry a bunch of autogenerated shims than adding to our wheel building load.

That reminds me: @Developer-Ecosystem-Engineering it would be good to add the script(s) used to generate the shim headers to this PR. That way we can regenerate them if needed in the future.

[mattip]

Is the intent that we stop shipping OpenBLAS in the macOS wheels?

[rgommers]

Is the intent that we stop shipping OpenBLAS in the macOS wheels?

To me, yes. That should be separate from this PR, but if Accelerate is significantly faster, passes all tests and results in much smaller wheels, I don't see why we shouldn't do that.

[seberg]

Well, the code is easy enough since its generated, but it still seems a bit strange to effectively ship a a blas "implementation" (which just dispatches, sure). So I am wondering:

• Isn't this exclusively needed for distributable wheels? I.e. there is already no problem to compile against Accelerate (or its a minor thing to enable). Just the wheel will not be portable if I create them using the 13.3 ABI.
• If/since I presume this is only really about wheels, there are a few more questions:
  • SciPy will need the same wrapper and so will others maybe, could this be not in the main codebase, but rather a (maybe header only) blas/lapack implementation like any other?
  • The dispatching code seems (I may be wrong) to fall back to the legacy Accelerate. Do we build wheels for any MacOS versions where our tests will not pass due to the old Accelerate being buggy?
  • A bit related: We are using 64bit interfaces currently with OpenBLAS. Even without bugs maybe that is actually a reason to use OpenBLAS for the wheels when 64bit interfaces are not available?

[mattip]

Compilation is failing:

INFO: compiling C++ sources
INFO: C compiler: g++ -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Werror=vla -Werror=nonnull -Werror=pointer-arith -Wlogical-op -Wno-sign-compare -Werror=undef -fPIC

INFO: compile options: '-DHAVE_CBLAS -Inumpy/core/include -Ibuild/src.linux-x86_64-3.9/numpy/core/include/numpy -Ibuild/src.linux-x86_64-3.9/numpy/distutils/include -Inumpy/core/src/common -Inumpy/core/src -Inumpy/core -Inumpy/core/src/npymath -Inumpy/core/src/multiarray -Inumpy/core/src/umath -Inumpy/core/src/npysort -Inumpy/core/src/_simd -I/home/runner/work/numpy/numpy/builds/venv/include -I/opt/hostedtoolcache/Python/3.9.17/x64/include/python3.9 -Ibuild/src.linux-x86_64-3.9/numpy/core/src/common -Ibuild/src.linux-x86_64-3.9/numpy/core/src/npymath -c'
extra options: '-fno-threadsafe-statics -D__STDC_VERSION__=0 -fno-exceptions -fno-rtti'
INFO: g++: numpy/linalg/umath_linalg.cpp
numpy/linalg/umath_linalg.cpp:301:1: error: expected unqualified-id before ‘{’ token
  301 | {
      | ^

[Developer-Ecosystem-Engineering]

Well, the code is easy enough since its generated, but it still seems a bit strange to effectively ship a a blas "implementation" (which just dispatches, sure). So I am wondering:

• Isn't this exclusively needed for distributable wheels? I.e. there is already no problem to compile against Accelerate (or its a minor thing to enable). Just the wheel will not be portable if I create them using the 13.3 ABI.

• If/since I presume this is only really about wheels, there are a few more questions:

  • SciPy will need the same wrapper and so will others maybe, could this be not in the main codebase, but rather a (maybe header only) blas/lapack implementation like any other?
  • The dispatching code seems (I may be wrong) to fall back to the legacy Accelerate. Do we build wheels for any MacOS versions where our tests will not pass due to the old Accelerate being buggy?
  • A bit related: We are using 64bit interfaces currently with OpenBLAS. Even without bugs maybe that is actually a reason to use OpenBLAS for the wheels when 64bit interfaces are not available?

• Fixed the print failures
• Provided a variant without the dispatch as there are open questions for the project to decide what might work best for its needs

[charris]

@Developer-Ecosystem-Engineering This needs a rebase.

[mattip]

We no longer require the header shims?
The ppc64le failure on travis is unrelated.

[rgommers]

Thanks for the discussion on this PR and in today's community meeting everyone. Looks like we're good to get this merged without the dispatching shims.

I've fixed the merge conflicts and added a CI job, which is now possible (still in beta, and the default XCode version was too old - but it was still pretty simple to get it to work).

It all LGTM, modulo one comment about the cblas.h check not being quite right. I'll revisit tomorrow.

The dispatching code seems (I may be wrong) to fall back to the legacy Accelerate. Do we build wheels for any MacOS versions where our tests will not pass due to the old Accelerate being buggy?

I haven't checked, but we were not forbidding use of Accelerate previously, and I don't think we've had bug reports in a long time.

A bit related: We are using 64bit interfaces currently with OpenBLAS. Even without bugs maybe that is actually a reason to use OpenBLAS for the wheels when 64bit interfaces are not available?

Indeed - I think that'd be a regression for macOS <13.3, and we can't really do that. And it seems like everyone was fine with extra wheels for >=13.3, so let's go that way (EDIT: should be macOS >=14.0, because Python packaging isn't smart enough to be able to target minor macOS versions)). This PR now has a pleasantly small diff.

[rgommers]

pleasantly fast:)

[rgommers]

One thing I noticed when adding the CI job is that xcrun -sdk macosx --show-sdk-version will show the SDK version which is controlled by XCode version - and that may not match the macOS version. The default GitHub Actions macos-13 image right now has macOS 13.5 but the SDK for 13.1, so the new Accelerate will not be picked up.

I believe that this is only a potential issue at build time, not at runtime. The installed Accelerate version should be tied to macOS version, not SDK version. If that wouldn't be the case, we'd have a problem - because we can only control wheel selection by macOS version. @Developer-Ecosystem-Engineering just to make sure: do I have this right? There is no way to get macOS 13.3 or 14.0 where the Accelerate shared library with the new symbols included will be missing, correct?

[andyfaff]

Typically those images have multiple sdk installed. So if you want to choose a different one you can.

[rgommers]

Typically those images have multiple sdk installed. So if you want to choose a different one you can.

Yes, that's what the CI job added in this PR does. I asked because I want to make doubly sure that we can ship *_macosx_14_0_arm64.whl's. I'm pretty sure that's right, but it's hard to verify.

[rgommers]

This all looks happy now, and the merge isn't conditional on the last question - that's for distributing wheels for 2.0, which we'll have to revisit in Oct/Nov when macOS 14.0 is released. So in it goes. Thanks a lot @Developer-Ecosystem-Engineering!

[rgommers]

@charris I've marked this for backporting, it'd be good to include these changes in 1.26.0rc1.

[mattip]

Is the idea that by releasing this in the RC1 we could get feedback about whether accelerate is included in various versions of MacOS? Or should we wait with the backport until we get a clear answer?

[rgommers]

We're not building wheels with Accelerate for 1.26.0, and as I said above the question I had is only relevant for that. The purpose of backporting is simply to complete build system support. It's actually easier to build against Accelerate than against OpenBLAS for users on up-to-date macOS versions. Also, Accelerate is actually more robust than OpenBLAS - I see local floating point error related issues with OpenBLAS on arm64, while Accelerate passes the full test suite.

[rgommers]

There is no way to get macOS 13.3 or 14.0 where the Accelerate shared library with the new symbols included will be missing, correct?

I got confirmation that this is correct, Accelerate is guaranteed to always be there for macOS >=13.3. We're all good here.

[charris]

@rgommers The backport adds a new file macos.yml that I assume is part of the CI cleanup in main. Are there other associated files that need to be backported/modified? I assume that the wheel builds are still done by cirrus, does anything need to be added for that?

tools/ci/cirrus_macosx_arm64.yml seems to have changed quite a bit and was renamed in main to cirrus_arm.yml.

[charris]

I'm thinking that both .cirrus.star and cirrus_arm.yml should be taken from main.

[rgommers]

Oh, the macos.yml can be added as is, or left out - whatever you prefer. They are new CI jobs that pass on main and will also pass on 1.26.x if added. Leaving them out won't do any harm though.

There is no impact on Cirrus CI config from this PR. However, yes indeed you should be able to copy .cirrus.star and cirrus_arm.yml directly, and that seems useful to me.

[charris]

macos.yml can be added as is, or left out

I kept it, but had to update environment.yml. Decided to keep the Cython>=3.0 there, didn't see any harm in that.