Add MKL BLAS64 support and multi-threaded PARDISO in wheels (#192)

* Add MKL BLAS64 support and multi-threaded PARDISO in wheels

Previously the MKL backend only supported LP64 (32-bit BLAS integers)
and shipped sequential MKL in pip wheels, giving single-threaded
PARDISO. This adds ILP64 support and ships multi-threaded MKL.

meson.build:
- Select LP64 or ILP64 MKL component libraries based on use_blas64
- Prefer iomp (Intel OpenMP) variant for multi-threaded MKL/PARDISO,
  fall back to sequential if iomp pkg-config is unavailable
- Windows uses sequential MKL due to a broken iomp .pc file in conda
- Add mkl_rt to MKL backend dependencies so MKL_Set_Interface_Layer()
  is available at runtime (symbol not exported by component libraries)

pyproject.toml:
- Add compiler pkgconfig/lib paths to cibuildwheel environment so
  mkl-dynamic-lp64-iomp is discoverable (its .pc file lives under
  the compiler runtime, not the MKL directory)
- Result: Linux x86_64 wheels now bundle libiomp5 + libmkl_intel_thread
  instead of libmkl_sequential, enabling multi-threaded PARDISO
- Install intel-openmp on Windows for future iomp support

build.yml:
- Add use_blas64 matrix dimension to build_mkl CI job
- Simplify build step to always pass link_mkl, conditionally add
  use_blas64

Depends on cvxgrp/scs#383 (updates scs_source submodule).

* Fail fast on unsupported MKL configurations

* Update scs_source MKL initialization order

* Stop enforcing MKL threading layer at runtime

* Document MKL wheel and BLAS64 contracts

* Update scs_source for CMake BLAS64 fix

* Fix pixi MKL dependency set in CI

* Update scs_source for expanded MKL CMake matrix

* Update scs_source for clearer MKL init error

* Update scs_source for MKL mismatch CI

* Update scs_source for verbose MKL mismatch test

* Update scs_source for mismatch test ABI fix

* Update scs_source to merged SCS master

* Simplify Meson mkl_rt lookup

Author: bodono

.github/workflows/build.yml

@@ -58,10 +58,12 @@ jobs:
         os: [ ubuntu-latest ]
         python-version: [ 3.9, "3.10", "3.11", "3.12", "3.13", "3.14"]
         link_mkl: [true]
+        use_blas64: [true, false]
 
     env:
       PYTHON_VERSION: ${{ matrix.python-version }}
       LINK_MKL: ${{ matrix.link_mkl }}
+      USE_BLAS64: ${{ matrix.use_blas64 }}
 
     steps:
       - uses: actions/checkout@v6
@@ -90,12 +92,11 @@ jobs:
       - name: Build
         run: |
           pixi r python -c "import numpy as np;print(np.show_config())"
-          if [[ "$LINK_MKL" == "true" ]]; then
-
-            pixi r python -m pip install --verbose --no-build-isolation -Csetup-args=-Dlink_mkl=true .
-          else
-            pixi r python -m pip install --verbose --no-build-isolation .
+          MESON_ARGS="-Csetup-args=-Dlink_mkl=true"
+          if [[ "$USE_BLAS64" == "true" ]]; then
+            MESON_ARGS="$MESON_ARGS -Csetup-args=-Duse_blas64=true"
           fi
+          pixi r python -m pip install --verbose --no-build-isolation $MESON_ARGS .
       - name: Test
         run: |
           pixi r pytest
@@ -198,7 +199,7 @@ jobs:
 
       - name: Install MKL from conda on Windows
         if: runner.os == 'Windows'
-        run: conda install -y -c https://software.repos.intel.com/python/conda/ -c conda-forge mkl mkl-devel pkgconfig
+        run: conda install -y -c https://software.repos.intel.com/python/conda/ -c conda-forge mkl mkl-devel intel-openmp pkgconfig
 
       - name: Build wheels
         uses: pypa/cibuildwheel@v3.4.1

README.md

@@ -48,6 +48,9 @@ backends can be enabled with build-time flags:
 # MKL Pardiso direct solver
 pip install . -Csetup-args=-Dlink_mkl=true
 
+# Use 64-bit BLAS/LAPACK integers (ILP64 / BLAS64)
+pip install . -Csetup-args=-Duse_blas64=true
+
 # GPU direct solver (cuDSS)
 pip install . -Csetup-args=-Dlink_cudss=true -Csetup-args=-Dint32=true
 
@@ -58,6 +61,11 @@ pip install . -Csetup-args=-Duse_lapack=true
 pip install . -Csetup-args=-Duse_spectral_cones=true
 ```
 
+Notes:
+- Linux x86_64 wheels are built and tested against threaded MKL, and CI asserts a `libiomp5` dependency on the packaged `_scs_mkl` extension. Windows currently falls back to sequential MKL because Intel's conda `pkg-config` metadata for the threaded variant is still broken.
+- `BLAS64` is a general SCS build mode for ILP64 BLAS/LAPACK libraries, not an MKL-only feature.
+- For the MKL Pardiso backend specifically, `BLAS64` must be paired with 64-bit SCS integers (`DLONG` / `int32=false`), and SCS now fails early if another library in the process has already fixed MKL to an incompatible LP64/ILP64 interface layer.
+
 ## Usage
 
 ```python

meson.build

@@ -22,19 +22,55 @@ print(incdir)
 
 # Get BLAS
 blas_deps = []
+mkl_pkg_name = ''
 if get_option('link_mkl')
-    # Prefer component libraries (lp64, sequential) so that wheel-repair
-    # tools (auditwheel / delvewheel) can detect and bundle them.
-    blas_deps = [dependency('mkl-dynamic-lp64-seq', required : false)]
+    if get_option('use_blas64') and get_option('int32')
+        error('MKL BLAS64 requires 64-bit SCS integers. Re-run Meson with -Dint32=false.')
+    endif
+    # Link against MKL component libraries. The integer width must match:
+    #   use_blas64=false (default) -> LP64  (32-bit BLAS integers)
+    #   use_blas64=true            -> ILP64 (64-bit BLAS integers)
+    #
+    # Note: PARDISO integer width is independent — it is controlled by the
+    # pardiso vs pardiso_64 entry point (selected via DLONG), not by the
+    # interface layer. See linsys/mkl/direct/private.c for details.
+    # Prefer Intel OpenMP (iomp) for multi-threaded MKL, fall back to
+    # sequential if iomp is not available.
+    #
+    # TODO: On Windows the conda MKL iomp .pc file (mkl-dynamic-*-iomp)
+    # has broken Cflags (raw path without -I prefix). Revisit once Intel
+    # ships a fixed pkg-config file; until then Windows gets sequential MKL.
+    if get_option('use_blas64')
+        _mkl_iomp = 'mkl-dynamic-ilp64-iomp'
+        _mkl_seq = 'mkl-dynamic-ilp64-seq'
+    else
+        _mkl_iomp = 'mkl-dynamic-lp64-iomp'
+        _mkl_seq = 'mkl-dynamic-lp64-seq'
+    endif
+    if host_machine.system() == 'windows'
+        # On Windows the conda MKL iomp .pc file has broken Cflags (raw path
+        # without -I prefix). Use sequential MKL until Intel ships a fix.
+        blas_deps = [dependency(_mkl_seq, required : false)]
+        if blas_deps[0].found()
+            mkl_pkg_name = _mkl_seq
+        endif
+    else
+        blas_deps = [dependency(_mkl_iomp, required : false)]
+        if not blas_deps[0].found()
+            blas_deps = [dependency(_mkl_seq, required : false)]
+            if blas_deps[0].found()
+                mkl_pkg_name = _mkl_seq
+            endif
+        else
+            mkl_pkg_name = _mkl_iomp
+        endif
+    endif
     if not blas_deps[0].found()
         blas_deps = [cc.find_library('mkl_rt', required : false)]
     endif
     if not blas_deps[0].found()
         blas_deps = [dependency('mkl-sdl', required : false)]
     endif
-    if not blas_deps[0].found()
-      blas_deps = [ dependency('mkl-dynamic-ilp64-iomp', required: false) ]
-    endif
 else
     if host_machine.system() == 'darwin'
         blas_deps = [dependency('Accelerate')]
@@ -276,16 +312,35 @@ if get_option('use_gpu')
 endif
 
 if get_option('link_mkl')
+  # The MKL backend calls MKL_Set_Interface_Layer() for a runtime interface
+  # check. This symbol is only exported by mkl_rt (the single dynamic library),
+  # not by component libraries (mkl_intel_lp64 etc.) that pkg-config may link.
+  # Add mkl_rt explicitly so the symbol is always available.
+  mkl_rt_dep = dependency('mkl-sdl', required : false)
+  if not mkl_rt_dep.found()
+    mkl_rt_dep = cc.find_library('mkl_rt', required : false)
+  endif
+  if not mkl_rt_dep.found() and mkl_pkg_name != ''
+    # Search next to the pkg-config-provided MKL component libraries.
+    _mkl_libdir = blas_deps[0].get_variable(pkgconfig: 'libdir', default_value: '')
+    if _mkl_libdir != ''
+      mkl_rt_dep = cc.find_library('mkl_rt', dirs: [_mkl_libdir], required : false)
+    endif
+  endif
+  if not mkl_rt_dep.found()
+    error('link_mkl=true requires mkl_rt (or mkl-sdl). _scs_mkl calls MKL_Set_Interface_Layer() at startup and cannot be built safely without it.')
+  endif
+  _mkl_deps = _deps + [mkl_rt_dep]
   py.extension_module(
     '_scs_mkl',
     'scs/scspy.c',
     'scs_source/linsys/mkl/direct/private.c',
     common_linsys_sources,
     scs_core_sources,
     spectral_cone_sources,
-    c_args: common_c_args + ['-DPY_MKL'],
+    c_args: common_c_args + ['-DPY_MKL', '-DSCS_MKL=1'],
     include_directories: common_includes + ['scs_source/linsys/mkl/direct'],
-    dependencies: _deps,
+    dependencies: _mkl_deps,
     install_dir: scs_dir,
     install: true,
   )

pyproject.toml

@@ -75,7 +75,7 @@ before-all = [
   "rpm --import https://yum.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB",
   "dnf install -y intel-oneapi-mkl-devel",
 ]
-environment = {PKG_CONFIG_PATH = "/opt/intel/oneapi/mkl/latest/lib/pkgconfig:$PKG_CONFIG_PATH", LD_LIBRARY_PATH = "/opt/intel/oneapi/mkl/latest/lib/intel64:/opt/intel/oneapi/mkl/latest/lib:$LD_LIBRARY_PATH"}
+environment = {PKG_CONFIG_PATH = "/opt/intel/oneapi/mkl/latest/lib/pkgconfig:/opt/intel/oneapi/compiler/latest/lib/pkgconfig:$PKG_CONFIG_PATH", LD_LIBRARY_PATH = "/opt/intel/oneapi/mkl/latest/lib/intel64:/opt/intel/oneapi/mkl/latest/lib:/opt/intel/oneapi/compiler/latest/lib:$LD_LIBRARY_PATH"}
 config-settings = {setup-args = "-Dlink_mkl=true"}
 
 [[tool.cibuildwheel.overrides]]

scs_source

@@ -1,4 +1,6 @@
 from __future__ import print_function, division
+import os
+import subprocess
 import sys
 import platform
 import scs
@@ -41,6 +43,15 @@
 params = {"verbose": True, "eps_abs": 1e-7, "eps_rel": 1e-7, "eps_infeas": 1e-7}
 
 
+@pytest.mark.skipif(
+    not (sys.platform == "linux" and platform.machine() == "x86_64"),
+    reason="Threaded MKL linkage is only checked on Linux x86_64",
+)
+def test_mkl_module_links_intel_openmp():
+    out = subprocess.check_output(["ldd", os.fspath(_scs_mkl.__file__)], text=True)
+    assert "libiomp5" in out
+
+
 def test_solve_feasible():
     data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.MKL, **params)