From 2495a5e77ec6ce8d4404e0c7b8622e21e87355e2 Mon Sep 17 00:00:00 2001 From: Fabian Isensee Date: Fri, 6 Nov 2026 12:22:29 +0100 Subject: [PATCH 1/3] Make the matmul, miniexpr and Menudet prefilters unable to need the GIL C-Blosc2 runs prefilters on its worker threads when nthreads > 1 (and on the calling thread without the GIL when compressions release it), but the error paths of the miniexpr and matmul prefilters raised Python exceptions: Cython took the GIL for the raise and the wrapper then called PyErr_Occurred() without it. With the caller holding the GIL, any error (e.g. an operand chunk that fails to decompress) deadlocked matmul and legacy-evaluator expressions and reductions with nthreads > 1; without it, it segfaulted. Declare miniexpr_prefilter, matmul_prefilter and their helpers `noexcept nogil` and remove every `raise` (Cython would still wrap one in an implicit `with gil`); use cdivision there and in the Menudet prefilter (the zero-division checks took the GIL as well). aux_miniexpr() records its error in its worker slot and returns -1; _store_data() and _prefilter_data() raise it (C-Blosc2 ignores the return code of a prefilter with a disposable output, i.e. of reductions). aux_matmul() frees what it allocated and returns -1. These prefilters, and the workers of parallel native graph maps, decode their inputs themselves, so operands with a codec or filter registered from Python ran that callback without the GIL (segfault; the Menudet reduction even with one thread). The routing gates (matmul fast path, _set_pref_expr(), _set_pref_matmul(), Menudet input eligibility, native graph layout) now decline such operands, so they take another path. matmul_prefilter raised "Unsupported dtype" for dtypes other than float and int; blosc2.matmul() never sent those, but it sent float16 and longdouble (kind "f"), which the kernels computed as float64: wrong results, and for float16 reads and writes past the block buffers. The fast path now takes float32/float64 only, and _set_pref_matmul() checks the dtypes before the prefilter is set. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01DiPFztdpnKMJkzrjT36BuP --- src/blosc2/blosc2_ext.pyx | 155 +++++++++++----- src/blosc2/linalg.py | 14 +- src/blosc2/native_graph.py | 3 + tests/ndarray/test_linalg.py | 83 +++++++++ tests/ndarray/test_prefilter_errors.py | 239 +++++++++++++++++++++++++ 5 files changed, 444 insertions(+), 50 deletions(-) create mode 100644 tests/ndarray/test_prefilter_errors.py diff --git a/src/blosc2/blosc2_ext.pyx b/src/blosc2/blosc2_ext.pyx index de90e254c..df68fcefe 100644 --- a/src/blosc2/blosc2_ext.pyx +++ b/src/blosc2/blosc2_ext.pyx @@ -858,6 +858,7 @@ ctypedef struct me_worker_s: void** owned_buffers int32_t* owned_capacities blosc2_context** decoders + const char* error # first error of aux_miniexpr() in this worker, see _raise_miniexpr_error() ctypedef struct me_udata: b2nd_array_t** inputs @@ -1354,6 +1355,7 @@ cdef extern from "fenv.h": int fesetenv(const fenv_t *) noexcept nogil +@cython.cdivision(True) # see aux_miniexpr() cdef int menudet_menudet_prefilter(blosc2_prefilter_params *params) noexcept nogil: cdef menudet_prefilter_state *state = params.user_data cdef menudet_prefilter_worker *worker = &state.workers[params.tid] @@ -3574,6 +3576,7 @@ cdef int64_t _store_data(blosc2_schunk *schunk, _write_op op, int64_t nchunk, da PyBuffer_Release(&buf) if rc < 0: _raise_callback_error() + _raise_miniexpr_error(schunk.storage.cparams) raise RuntimeError(error) return rc @@ -4595,6 +4598,7 @@ cdef class SChunk: finally: PyBuffer_Release(&buf) PyBuffer_Release(&chunk_buf) + _raise_miniexpr_error(self.schunk.storage.cparams) _check_rc(size, "Could not compress the data") if size == 0: raise RuntimeError("The result could not fit ") @@ -5045,11 +5049,19 @@ cdef int general_filler(blosc2_prefilter_params *params): return 0 +cdef inline int _miniexpr_fail(me_worker_s *worker, const char *error) noexcept nogil: + """Error exit of aux_miniexpr(): keep the worker's first `error` for _raise_miniexpr_error().""" + if worker.error == NULL: + worker.error = error + return -1 + # Auxiliary function for miniexpr as a prefilter # Only meant for (input and output) arrays that are blosc2.NDArray objects. +# cdivision: ZeroDivisionError checks take the GIL; all operands are sizes and indices, divisors > 0. +@cython.cdivision(True) cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, c_bool is_postfilter, uint8_t *params_output, int32_t typesize, - int32_t tid) nogil: + int32_t tid) noexcept nogil: # Declare all C variables at the beginning cdef int64_t chunk_ndim[B2ND_MAX_DIM] cdef int64_t block_ndim[B2ND_MAX_DIM] @@ -5088,12 +5100,12 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, cdef void* aux_reduc_ptr if miniexpr_handle == NULL: - raise ValueError("miniexpr: handle not assigned") + return _miniexpr_fail(worker, "miniexpr: handle not assigned") # Query valid (unpadded) items for this block rc = me_nd_valid_nitems(miniexpr_handle, nchunk, nblock, &valid_nitems) if rc != 0: - raise RuntimeError(f"miniexpr: invalid block; error code: {rc}") + return _miniexpr_fail(worker, "miniexpr: invalid block") if valid_nitems <= 0: # Nothing to compute for this block. # For reductions, keep aux_reduc neutral values untouched. @@ -5119,7 +5131,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, worker.decoders = calloc(udata.ninputs, sizeof(blosc2_context*)) if (worker.input_buffers == NULL or worker.owned_buffers == NULL or worker.owned_capacities == NULL or worker.decoders == NULL): - raise MemoryError("miniexpr: cannot allocate input buffer table") + return _miniexpr_fail(worker, "miniexpr: cannot allocate input buffer table") input_buffers = worker.input_buffers for i in range(udata.ninputs): @@ -5133,7 +5145,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, if expected_blocknitems == -1: expected_blocknitems = blocknitems elif blocknitems != expected_blocknitems: - raise ValueError("miniexpr: inconsistent block element counts across inputs") + return _miniexpr_fail(worker, "miniexpr: inconsistent block element counts across inputs") np_ndim = udata.array.ndim if np_ndim == 1 and valid_nitems == blocknitems: start_ndim[0] = nchunk * udata.array.chunkshape[0] + nblock * udata.array.blockshape[0] @@ -5145,7 +5157,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, worker.owned_buffers[i] = malloc(block_nbytes) if worker.owned_buffers[i] == NULL: worker.owned_capacities[i] = 0 - raise MemoryError("miniexpr: cannot allocate input block buffer") + return _miniexpr_fail(worker, "miniexpr: cannot allocate input block buffer") worker.owned_capacities[i] = block_nbytes input_buffers[i] = worker.owned_buffers[i] if valid_nitems < blocknitems: @@ -5206,7 +5218,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, else: input_cache = &udata.input_chunk_caches[i] if input_cache.state_lock == NULL or input_cache.ready_lock == NULL: - raise MemoryError("miniexpr: cache locks not assigned") + return _miniexpr_fail(worker, "miniexpr: cache locks not assigned") while True: PyThread_acquire_lock(input_cache.state_lock, 1) if input_cache.state == ME_CACHE_READY and input_cache.nchunk == nchunk and input_cache.data != NULL: @@ -5215,7 +5227,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, break if input_cache.state == ME_CACHE_ERROR and input_cache.nchunk == nchunk: PyThread_release_lock(input_cache.state_lock) - raise ValueError("miniexpr: error getting chunk") + return _miniexpr_fail(worker, "miniexpr: error getting chunk") if input_cache.state == ME_CACHE_LOADING: PyThread_release_lock(input_cache.state_lock) PyThread_acquire_lock(input_cache.ready_lock, 1) @@ -5232,7 +5244,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, input_cache.state = ME_CACHE_ERROR PyThread_release_lock(input_cache.state_lock) PyThread_release_lock(input_cache.ready_lock) - raise ValueError("miniexpr: error getting chunk") + return _miniexpr_fail(worker, "miniexpr: error getting chunk") if not needs_free: loaded_chunk = malloc(rc) @@ -5241,7 +5253,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, input_cache.state = ME_CACHE_ERROR PyThread_release_lock(input_cache.state_lock) PyThread_release_lock(input_cache.ready_lock) - raise MemoryError("miniexpr: cannot allocate chunk copy") + return _miniexpr_fail(worker, "miniexpr: cannot allocate chunk copy") memcpy(loaded_chunk, chunk, rc) else: loaded_chunk = chunk @@ -5258,16 +5270,16 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, break rc = blosc2_cbuffer_sizes(src, &chunk_nbytes, &chunk_cbytes, &block_nbytes) if rc < 0: - raise ValueError("miniexpr: error getting cbuffer sizes") + return _miniexpr_fail(worker, "miniexpr: error getting cbuffer sizes") if block_nbytes <= 0: - raise ValueError("miniexpr: invalid block size") + return _miniexpr_fail(worker, "miniexpr: invalid block size") if worker.owned_capacities[i] < block_nbytes: if worker.owned_buffers[i] != NULL: free(worker.owned_buffers[i]) worker.owned_buffers[i] = malloc(block_nbytes) if worker.owned_buffers[i] == NULL: worker.owned_capacities[i] = 0 - raise MemoryError("miniexpr: cannot allocate input block buffer") + return _miniexpr_fail(worker, "miniexpr: cannot allocate input block buffer") worker.owned_capacities[i] = block_nbytes input_buffers[i] = worker.owned_buffers[i] input_typesize = ndarr.sc.typesize @@ -5275,14 +5287,14 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, if expected_blocknitems == -1: expected_blocknitems = blocknitems elif blocknitems != expected_blocknitems: - raise ValueError("miniexpr: inconsistent block element counts across inputs") + return _miniexpr_fail(worker, "miniexpr: inconsistent block element counts across inputs") if worker.decoders[i] == NULL: worker.decoders[i] = blosc2_create_dctx(BLOSC2_DPARAMS_DEFAULTS) if worker.decoders[i] == NULL: - raise MemoryError("miniexpr: cannot allocate decompression context") + return _miniexpr_fail(worker, "miniexpr: cannot allocate decompression context") dctx = worker.decoders[i] if valid_nitems > blocknitems: - raise ValueError("miniexpr: valid items exceed padded block size") + return _miniexpr_fail(worker, "miniexpr: valid items exceed padded block size") # Ask in bytes: blosc2_getitem_ctx() counts in the typesize the *chunk* # records, which c-blosc2 caps to 1 above BLOSC_MAX_TYPESIZE (255), so its # unit changes silently with the data -- every block past the first once @@ -5291,7 +5303,7 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, rc = blosc2_getitem_bytes_ctx(dctx, src, chunk_cbytes, nblock * block_nbytes, block_nbytes, input_buffers[i], block_nbytes) if rc != block_nbytes: - raise ValueError("miniexpr: error decompressing the chunk") + return _miniexpr_fail(worker, "miniexpr: error decompressing the chunk") # For reduction operations, we need to track which block we're processing # The linear_block_index should be based on the same grid the output shares # with every input (raw NumPy inputs have no b2nd_array_t to read ndim from). @@ -5320,11 +5332,11 @@ cdef int aux_miniexpr(me_udata *udata, int64_t nchunk, int32_t nblock, rc = me_eval_nd(miniexpr_handle, input_buffers, udata.ninputs, aux_reduc_ptr, output_blocknitems, nchunk, nblock, udata.eval_params) if rc != 0: - raise RuntimeError(f"miniexpr: issues during evaluation; error code: {rc}") + return _miniexpr_fail(worker, "miniexpr: issues during evaluation") return 0 -cdef int matmul_block_kernel(T* A, T* B, T* C, int M, int K, int N) nogil: +cdef int matmul_block_kernel(T* A, T* B, T* C, int M, int K, int N) noexcept nogil: cdef int r, c, k cdef T a cdef int rowA, rowC, rowB @@ -5338,19 +5350,28 @@ cdef int matmul_block_kernel(T* A, T* B, T* C, int M, int K, int N) nogil: C[rowC + c] += (a * B[rowB + c]) return 0 -cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *params_output, int32_t typesize, int typecode) nogil: +cdef inline int _matmul_exit(void **input_buffers, blosc2_context *dctx, int rc) noexcept nogil: + """Free what aux_matmul() allocated and return `rc`.""" + if dctx != NULL: + blosc2_free_ctx(dctx) + free(input_buffers[0]) + free(input_buffers[1]) + return rc + +@cython.cdivision(True) # see aux_miniexpr() +cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *params_output, int32_t typesize, int typecode) noexcept nogil: # Declare all C variables at the beginning cdef b2nd_array_t* out_arr cdef b2nd_array_t* ndarr cdef c_bool first_run cdef int rc, p, q, r - cdef void** input_buffers = malloc(2 * sizeof(uint8_t*)) - cdef uint8_t** src = malloc(2 * sizeof(uint8_t*)) + cdef void* input_buffers[2] + cdef uint8_t* src[2] cdef int32_t chunk_nbytes[2] cdef int32_t chunk_cbytes[2] cdef int32_t block_nbytes[2] cdef int startA, startB, expected_blocknitems - cdef blosc2_context* dctx + cdef blosc2_context* dctx = NULL cdef int i, j, block_i, block_j, chunk_i, chunk_j, ncols, block_ncols, Bblock_ncols, Bncols, Ablock_ncols, Ancols cdef int nchunkA = 0, nchunkB = 0, nblockA = 0, nblockB = 0, offsetA = 0, offsetB = 0, offset = 0 out_arr = udata.array @@ -5359,6 +5380,7 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param cdef int coord, batch, batch_, batches = 1 cdef int out_chunk_nrows, out_chunk_ncols, out_block_nrows, out_block_ncols cdef int selected_backend = b2_get_selected_matmul_backend() + input_buffers[0] = input_buffers[1] = NULL # batches = sum(strides[i]*elcoords[i]) for i in range(ndim - 2): @@ -5394,6 +5416,8 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param memset(params_output, 0, out_arr.blocknitems * typesize) dctx = blosc2_create_dctx(BLOSC2_DPARAMS_DEFAULTS) + if dctx == NULL: + return _matmul_exit(input_buffers, dctx, -1) first_run = True while True: # chunk loop @@ -5403,10 +5427,8 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param ndim = ndarr.ndim src[i] = ndarr.sc.data[chunk_idx] rc = blosc2_cbuffer_sizes(src[i], &chunk_nbytes[i], &chunk_cbytes[i], &block_nbytes[i]) - if rc < 0: - raise ValueError("miniexpr: error getting cbuffer sizes") - if block_nbytes[i] <= 0: - raise ValueError("miniexpr: invalid block size") + if rc < 0 or block_nbytes[i] <= 0: + return _matmul_exit(input_buffers, dctx, -1) if first_run: if i == 0: q = ndarr.blockshape[ndim - 1] @@ -5431,7 +5453,7 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param block_startB = nblockB + block_j // r input_buffers[i] = malloc(block_nbytes[i]) if input_buffers[i] == NULL: - raise MemoryError("miniexpr: cannot allocate input block buffer") + return _matmul_exit(input_buffers, dctx, -1) first_run = False nblockA = block_startA @@ -5446,11 +5468,11 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param rc = blosc2_getitem_bytes_ctx(dctx, src[0], chunk_cbytes[0], startA, block_nbytes[0], input_buffers[0], block_nbytes[0]) if rc != block_nbytes[0]: - raise ValueError("matmul: error decompressing the A chunk") + return _matmul_exit(input_buffers, dctx, -1) rc = blosc2_getitem_bytes_ctx(dctx, src[1], chunk_cbytes[1], startB, block_nbytes[1], input_buffers[1], block_nbytes[1]) if rc != block_nbytes[1]: - raise ValueError("matmul: error decompressing the B chunk") + return _matmul_exit(input_buffers, dctx, -1) batch = 0 while batch < batches: batch_ = batch @@ -5526,8 +5548,7 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param else: rc = matmul_block_kernel[int64_t](input_buffers[0] + offsetA, input_buffers[1] + offsetB, params_output + offset, p, q, r) else: - with gil: - raise ValueError("Unsupported dtype") + return _matmul_exit(input_buffers, dctx, -1) batch += 1 nblockA += 1 nblockB += Bblock_ncols @@ -5539,14 +5560,7 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param break - blosc2_free_ctx(dctx) - # Free resources - for i in range(2): - free(input_buffers[i]) - free(input_buffers) - free(src) - - return 0 + return _matmul_exit(input_buffers, dctx, 0) # Aux function for prefilter and postfilter udf cdef int aux_udf(udf_udata *udata, int64_t nchunk, int32_t nblock, @@ -5621,11 +5635,16 @@ cdef int aux_udf(udf_udata *udata, int64_t nchunk, int32_t nblock, return 0 -cdef int miniexpr_prefilter(blosc2_prefilter_params *params): +# C-Blosc2 calls these two from its threads, which do not hold the GIL (_check_cparams() allows +# nthreads > 1 for them), so they and their helpers must not touch Python objects: they are +# noexcept nogil (Cython rejects Python code in them, except `raise`, which it would wrap in an +# implicit `with gil`, so there is none) and return -1 on errors. The operands they decompress +# must not use Python codecs or filters either (_reads_need_python()). +cdef int miniexpr_prefilter(blosc2_prefilter_params *params) noexcept nogil: return aux_miniexpr( params.user_data, params.nchunk, params.nblock, False, params.output, params.output_typesize, params.tid) -cdef int matmul_prefilter(blosc2_prefilter_params *params): +cdef int matmul_prefilter(blosc2_prefilter_params *params) noexcept nogil: cdef int typecode cdef mm_udata* udata = params.user_data @@ -5636,9 +5655,41 @@ cdef int matmul_prefilter(blosc2_prefilter_params *params): elif dtype_kind == 'i': typecode = 1 else: - raise ValueError("Unsupported dtype") + return -1 # _set_pref_matmul() lets no other dtype through return aux_matmul(udata, params.nchunk, params.nblock, params.output, params.output_typesize, typecode) +cdef inline c_bool _reads_need_python(blosc2_schunk *sc) noexcept: + """Whether reading `sc` runs a codec or filter registered from Python, which the miniexpr, matmul + and Menudet prefilters cannot do (they decompress their inputs on C-Blosc2's threads). Judged by + the super-chunk's cparams, not by each chunk's header: a chunk written with other cparams + (update_chunk()) is not caught, but registering a Python codec or filter somewhere does not + disable these paths for every other array (unlike _decompress_ctx()).""" + return _python_codec_ids[sc.storage.cparams.compcode] or _uses_python_filter(sc.storage.cparams) + +def reads_need_python(NDArray array): + """_reads_need_python() for the routing gates in Python (matmul, native graphs).""" + return _reads_need_python(array.array.sc) + +cdef _raise_miniexpr_error(blosc2_cparams *cparams): + """Raise the first error that aux_miniexpr() recorded in a worker slot, if any, and clear them. + Call it after each compression with miniexpr_prefilter: C-Blosc2 ignores the prefilter's return + code when the output is disposable (reductions), and only reports a generic failure otherwise.""" + cdef me_udata* udata + cdef const char* error = NULL + cdef int tid + if cparams.prefilter != miniexpr_prefilter or cparams.preparams == NULL: + return + udata = cparams.preparams.user_data + if udata == NULL: + return + for tid in range(256): + if udata.workers[tid].error != NULL: + if error == NULL: + error = udata.workers[tid].error + udata.workers[tid].error = NULL + if error != NULL: + raise RuntimeError(error.decode()) + cdef int general_udf_prefilter(blosc2_prefilter_params *params): cdef udf_udata *udata = params.user_data return aux_udf(udata, params.nchunk, params.nblock, False, params.output, params.output_typesize) @@ -7204,6 +7255,11 @@ cdef class NDArray: @_handle_locked def _set_pref_expr(self, expression, inputs, fp_accuracy, aux_reduc=None, jit=None, candidate_blocks=None): + for operand in inputs.values(): + if not isinstance(operand, np.ndarray) and _reads_need_python( + (operand.c_array).sc): + raise ValueError("miniexpr cannot read operands that use a codec or filter registered " + "from Python") # Set prefilter for miniexpr cdef blosc2_cparams* cparams = self.array.sc.storage.cparams cparams.prefilter = miniexpr_prefilter @@ -7391,12 +7447,23 @@ cdef class NDArray: @_handle_locked def _set_pref_matmul(self, inputs, fp_accuracy): cdef b2nd_array_t* operand + # matmul_prefilter() has kernels for these only, and cannot raise + if self.dtype not in (np.float32, np.float64, np.int32, np.int64) or any( + inputs[name].dtype != self.dtype for name in ("x1", "x2")): + raise TypeError( + "The fast matmul path needs operands and output of one dtype out of float32, float64, " + f"int32, int64 (native byte order), not {inputs['x1'].dtype}, {inputs['x2'].dtype} -> " + f"{self.dtype}" + ) for name in ("x1", "x2"): operand = inputs[name].c_array # aux_matmul() reads the chunks from sc.data, which frames (on disk, or in memory # with contiguous=True) do not fill if operand.sc.frame != NULL: raise ValueError("The fast matmul path does not support frame-backed operands") + if _reads_need_python(operand.sc): + raise ValueError("The fast matmul path does not support operands that use a codec or " + "filter registered from Python") # Set prefilter for miniexpr cdef blosc2_cparams* cparams = self.array.sc.storage.cparams cparams.prefilter = matmul_prefilter @@ -8130,6 +8197,8 @@ def menudet_prefilter_inputs_eligible(inputs, allow_numpy=False, output=None, ou return False if value.array.sc.storage.urlpath != NULL or value.array.sc.storage.dparams.postfilter != NULL: return False + if _reads_need_python(value.array.sc): # the prefilter decodes inputs on C-Blosc2's threads + return False if source.schunk.cparams.use_dict: return False block_bytes += int(np.prod(source.blocks)) * source.dtype.itemsize diff --git a/src/blosc2/linalg.py b/src/blosc2/linalg.py index 1dd4fd14a..3d4c16740 100644 --- a/src/blosc2/linalg.py +++ b/src/blosc2/linalg.py @@ -78,8 +78,11 @@ def _matmul_can_use_fast_path( return False # The prefilter reads the operands' chunks straight from memory, which frames (on disk, - # or in memory with contiguous=True) do not keep - if any(op.schunk.urlpath is not None or op.schunk.contiguous for op in (x1, x2)): + # or in memory with contiguous=True) do not keep, and decompresses them without the GIL + if any( + op.schunk.urlpath is not None or op.schunk.contiguous or blosc2.blosc2_ext.reads_need_python(op) + for op in (x1, x2) + ): return False # The current prefilter-backed implementation is only supported for 2-D layouts. @@ -108,11 +111,8 @@ def _matmul_can_use_fast_path( if result_blocks[:-2] != result.blocks[:-2]: return False - if x1.dtype.kind != "f": - return False - if x2.dtype.kind != "f": - return False - return x1.dtype == x2.dtype + # The prefilter has float32 and float64 kernels only (float16 and longdouble are kind "f" too) + return x1.dtype in (np.float32, np.float64) def _matmul_fast_path_context(result: blosc2.NDArray): diff --git a/src/blosc2/native_graph.py b/src/blosc2/native_graph.py index 638daee3a..0e76a76d5 100644 --- a/src/blosc2/native_graph.py +++ b/src/blosc2/native_graph.py @@ -1116,6 +1116,9 @@ def _compressed_execution_layout(plan, inputs, shape, *, reduction=True): compressed = [value for value in inputs.values() if type(value) is blosc2.NDArray] if len(shape) != 1 or not compressed or blosc2.IS_WASM or plan.info()["stages"] != 1: return None + # The workers decompress the inputs without the GIL + if any(blosc2.blosc2_ext.reads_need_python(value) for value in compressed): + return None minimum = MINIMUM_COMPRESSED_GRAPH_ITEMS if reduction else _compressed_map_minimum(inputs) if shape[0] < minimum: return None diff --git a/tests/ndarray/test_linalg.py b/tests/ndarray/test_linalg.py index a93b0b1c2..21269d044 100644 --- a/tests/ndarray/test_linalg.py +++ b/tests/ndarray/test_linalg.py @@ -437,6 +437,89 @@ def failing_set_pref_matmul(self, inputs, fp_accuracy): _toggle_miniexpr(old_flag) +# The matmul prefilter runs on C-Blosc2's threads when nthreads > 1 (none on WASM) +PREFILTER_NTHREADS = [1, pytest.param(4, marks=pytest.mark.skipif(blosc2.IS_WASM, reason="no threads"))] + + +def _int_valued_operands(dtype, nthreads, shape=(200, 200)): + # One chunk with 16 blocks per operand; integer values keep every product sum exact + rng = np.random.default_rng(0) + a = rng.integers(-2, 3, shape).astype(dtype) + b = rng.integers(-2, 3, shape).astype(dtype) + cparams = blosc2.CParams(nthreads=nthreads) + x1 = blosc2.asarray(a, chunks=shape, blocks=(50, 50), cparams=cparams) + x2 = blosc2.asarray(b, chunks=shape, blocks=(50, 50), cparams=cparams) + return a, b, x1, x2, cparams + + +# float16 and longdouble are kind "f" too, but the prefilter has float32/float64 kernels only +@pytest.mark.parametrize( + "dtype", + [ + np.float16, + pytest.param( + np.longdouble, + marks=pytest.mark.skipif(np.dtype(np.longdouble).itemsize <= 8, reason="longdouble is float64"), + ), + ], +) +def test_matmul_fast_path_skips_float_dtypes_without_kernel(monkeypatch, dtype): + old_flag = utils_mod.try_miniexpr + calls = _set_pref_matmul_call_recorder(monkeypatch) + try: + _toggle_miniexpr(True) + a, b, x1, x2, _ = _int_valued_operands(dtype, nthreads=1) + c = blosc2.matmul(x1, x2, chunks=(200, 200), blocks=(50, 50)) + + assert calls == [] + assert c.dtype == dtype + np.testing.assert_array_equal(c[:], np.matmul(a, b)) + finally: + _toggle_miniexpr(old_flag) + + +@pytest.mark.parametrize( + ("dtype1", "dtype2", "out_dtype"), + [ + (np.complex128, np.complex128, np.complex128), + (np.bool_, np.bool_, np.bool_), + (np.int16, np.int16, np.int16), + (np.float16, np.float16, np.float16), + ([("x", " +# All rights reserved. +# +# SPDX-License-Identifier: BSD-3-Clause +####################################################################### + +# Errors inside the miniexpr, matmul and Menudet prefilters. C-Blosc2 runs them on its own threads +# when nthreads > 1, and on the calling thread without the GIL when set_releasegil(True): they must +# report errors without needing the GIL, and must not run Python codecs or filters of their inputs. +# Raising there, or running such a callback, used to deadlock or crash the interpreter, +# so each case runs in its own child interpreter with a timeout: a regression then fails the cases +# it breaks instead of hanging or killing the test run (a case takes about 1 s). WASM has neither +# subprocesses nor threads: only the cases that keep the GIL run there, in-process. + +import faulthandler +import json +import math +import os +import struct +import subprocess +import sys +import warnings +import zlib + +import numpy as np +import pytest + +import blosc2 +from blosc2.utils import _toggle_miniexpr + +TIMEOUT = 30 # seconds +# The miniexpr prefilter only runs with the legacy expression evaluator (or without native Menudet +# support), with NumPy >= 2 (NumPy 1 evaluates numeric expressions with the typed evaluator) and, +# for reductions, not on WASM; the typed evaluator (default) uses the Menudet prefilter for 1-D +# operands. typed-getitem is large enough for the parallel native graph map (none on WASM), whose +# workers decode the inputs. +NO_LEGACY_MINIEXPR = pytest.mark.skipif( + blosc2.IS_WASM or np.lib.NumpyVersion(np.__version__) < "2.0.0", + reason="the legacy evaluator does not use the miniexpr prefilter here", +) +ERROR_OPS = [ + "matmul", + pytest.param("legacy-expr", marks=NO_LEGACY_MINIEXPR), + pytest.param("legacy-sum", marks=NO_LEGACY_MINIEXPR), +] +# The gate looks at the codec and at the filters: the filter on every route, the codec on two +CALLBACK_CASES = [ + *[(op, "filter") for op in ["matmul", "legacy-expr", "legacy-sum", "typed-expr", "typed-sum"]], + pytest.param( + "typed-getitem", "filter", marks=pytest.mark.skipif(blosc2.IS_WASM, reason="no native graphs") + ), + ("matmul", "codec"), + ("typed-sum", "codec"), +] +NTHREADS = [1, pytest.param(4, marks=pytest.mark.skipif(blosc2.IS_WASM, reason="no threads"))] +RELEASEGIL = [ + pytest.param(False, id="gil"), + pytest.param(True, id="releasegil", marks=pytest.mark.skipif(blosc2.IS_WASM, reason="no subprocesses")), +] +# Codec and filter ids are process-global and cannot be unregistered: ones no other test uses +CODEC_ID, FILTER_ID = 248, 249 +# Chunk header: flags byte (BLOSC_MEMCPYED bit), and the special-value bits of the blosc2 flags byte +CHUNK_FLAGS, MEMCPYED = 2, 0x2 +CHUNK_BLOSC2_FLAGS, SPECIAL_SHIFT, SPECIAL_MASK = 31, 4, 0x7 +PAST_ANY_CHUNK = 0x7FFFFFF0 + + +def _break_first_chunk(array): + """Point every block offset of chunk 0 past the chunk: its header still parses, but reading + any block fails.""" + chunk = bytearray(array.schunk.get_chunk(0)) + assert not chunk[CHUNK_FLAGS] & MEMCPYED, "a memcpyed chunk has no block offsets" + special = (chunk[CHUNK_BLOSC2_FLAGS] >> SPECIAL_SHIFT) & SPECIAL_MASK + assert special == 0, "a special-value chunk has no block offsets" + nbytes, _, blocksize = blosc2.get_cbuffer_sizes(chunk) + for nblock in range(math.ceil(nbytes / blocksize)): # the offsets follow the header + struct.pack_into("= output.size: + return 0 + output[: len(out)] = np.frombuffer(out, dtype=np.uint8) + return len(out) + + +def _decode(input, output, meta, schunk): + out = zlib.decompress(input.tobytes()) + output[: len(out)] = np.frombuffer(out, dtype=np.uint8) + return len(out) + + +def python_callback_case(op, nthreads, kind): + """`op` on operands compressed with a codec or filter (`kind`) registered from Python: whether + the result is right, and the warnings.""" + if kind == "codec": + if CODEC_ID not in blosc2.ucodecs_registry: + blosc2.register_codec("zlib_prefilter_test", CODEC_ID, _encode, _decode) + op_cparams = blosc2.CParams(nthreads=1, codec=CODEC_ID, filters=[0], filters_meta=[0]) + else: + if FILTER_ID not in blosc2.ufilters_registry: + blosc2.register_filter(FILTER_ID, _plus_one, _minus_one, "plus_one_prefilter_test") + op_cparams = blosc2.CParams(nthreads=1, filters=[FILTER_ID], filters_meta=[0]) + a, b, x1, x2 = _operands(op, cparams=op_cparams, dparams={"nthreads": 1}) + messages = [] + result = _run_op(op, x1, x2, blosc2.CParams(nthreads=nthreads), messages) + expected = a @ b if op == "matmul" else (a + b).sum() if op.endswith("sum") else a + b + return {"max_abs_diff": float(np.max(np.abs(result - expected))), "warnings": messages} + + +def _run_in_child(case, *args): + """Run `case(*args)` in a child interpreter, which must finish within TIMEOUT.""" + if blosc2.IS_WASM: + return case(*args) + env = {name: value for name, value in os.environ.items() if name != "BLOSC_NTHREADS"} + try: + proc = subprocess.run( + [sys.executable, __file__, case.__name__, json.dumps(args)], + capture_output=True, + text=True, + env=env, + check=False, + timeout=TIMEOUT, + ) + except subprocess.TimeoutExpired: # the child's faulthandler normally ends it first + pytest.fail(f"{case.__name__}{args} did not finish within {TIMEOUT} s (deadlock?)") + # A deadlock shows up as "Timeout" with the stacks of the child's threads + assert proc.returncode == 0, f"{case.__name__}{args}: exit code {proc.returncode}\n{proc.stderr[-3000:]}" + return json.loads(proc.stdout.splitlines()[-1]) + + +@pytest.mark.parametrize("releasegil", RELEASEGIL) +@pytest.mark.parametrize("nthreads", NTHREADS) +@pytest.mark.parametrize("op", ERROR_OPS) +def test_prefilter_error_raises(op, nthreads, releasegil): + outcome = _run_in_child(broken_operand_case, op, nthreads, releasegil) + assert outcome["raised"] == "RuntimeError", outcome + if op != "matmul": + # The prefilter's own error, also for the reduction, where C-Blosc2 ignores the return code + # of the prefilter (its output is disposable) + assert outcome["message"] == "miniexpr: error decompressing the chunk", outcome + assert outcome["warnings"] == [], outcome + else: + # It falls back to its chunked path (with a warning), which then reads the broken chunk + assert len(outcome["warnings"]) == 1, outcome + assert outcome["warnings"][0].startswith("RuntimeWarning: "), outcome + assert "falling back to chunked path" in outcome["warnings"][0], outcome + + +@pytest.mark.parametrize("nthreads", NTHREADS) +@pytest.mark.parametrize(("op", "kind"), CALLBACK_CASES) +def test_prefilter_skips_python_callback_operands(op, nthreads, kind): + # The prefilters would run the codec's decoder or the filter's backward function on + # C-Blosc2's threads, which do not hold the GIL + outcome = _run_in_child(python_callback_case, op, nthreads, kind) + assert outcome["max_abs_diff"] == 0, outcome + assert outcome["warnings"] == [], outcome + + +if __name__ == "__main__": + faulthandler.dump_traceback_later(TIMEOUT - 5, exit=True) + print(json.dumps(globals()[sys.argv[1]](*json.loads(sys.argv[2])))) From 2b171f59bdc66d8159f79f293f2bc1fa46db69c4 Mon Sep 17 00:00:00 2001 From: Fabian Isensee Date: Fri, 6 Nov 2026 09:15:27 +0100 Subject: [PATCH 2/3] Fix wrong matmul fast-path results for multi-chunk operands aux_matmul() read src[i] = ndarr.sc.data[chunk_idx] before its first_run branch moved nchunkA/nchunkB to the output chunk's row and column, so the first K step of every output chunk multiplied A's and B's first chunks. Only layouts with one output chunk came out right; with non-constant data (the tests used ones()/full()) a 200x200 float64 product with 100x100 chunks was off by up to 226. Compute the first chunk and block indices before the chunk loop; they depend only on the shapes. The kernel also assumes that A, B and the output share one chunk grid (output chunk rows = A's, output chunk columns = B's, same chunks along K), but _matmul_can_use_fast_path() did not check that. Such layouts gave wrong results, failed in the prefilter (and in the chunked fallback, which still runs with the prefilter set) or crashed reading past the operand's chunk list. Require the shared grid in the gate; other layouts take the chunked path. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01DiPFzsfE5i74cBSxZN13M8 --- src/blosc2/blosc2_ext.pyx | 46 ++++++++++++----------- src/blosc2/linalg.py | 8 ++++ tests/ndarray/test_linalg.py | 73 ++++++++++++++++++++++++++++++++++++ 3 files changed, 106 insertions(+), 21 deletions(-) diff --git a/src/blosc2/blosc2_ext.pyx b/src/blosc2/blosc2_ext.pyx index df68fcefe..baf0a5f89 100644 --- a/src/blosc2/blosc2_ext.pyx +++ b/src/blosc2/blosc2_ext.pyx @@ -5413,6 +5413,31 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param out_block_nrows = out_arr.blockshape[ndim - 2] out_block_ncols = out_arr.blockshape[ndim - 1] + # Move nchunkA/nchunkB to the first chunks this output chunk needs: A's chunk in the + # same row and B's chunk in the same column, both at K = 0. The chunk loop reads src[] + # from them, so this must happen before it starts. + ndarr = udata.inputs[0] + ndim = ndarr.ndim + q = ndarr.blockshape[ndim - 1] + p = ndarr.blockshape[ndim - 2] + # nchunk_ = chunks_in_row * chunk_row + chunk_col + # convert from chunk_idx to element idx chunk_i (row) + chunk_i = nchunk_ // ncols * out_chunk_nrows + nchunkA += chunk_i // ndarr.chunkshape[ndim - 2] * Ancols + # nblock_ = blocks_in_chunkrow * block_row + block_col + # convert from block_idx to element idx block_i (row) + block_i = nblock_ // block_ncols * out_block_nrows + block_startA = nblockA + block_i // p * Ablock_ncols + ndarr = udata.inputs[1] + ndim = ndarr.ndim + r = ndarr.blockshape[ndim - 1] + # convert from chunk_idx to element idx chunk_j (col) + chunk_j = nchunk_ % ncols * out_chunk_ncols + nchunkB += chunk_j // ndarr.chunkshape[ndim - 1] + # convert from block_idx to element idx block_j (col) + block_j = nblock_ % block_ncols * out_block_ncols + block_startB = nblockB + block_j // r + memset(params_output, 0, out_arr.blocknitems * typesize) dctx = blosc2_create_dctx(BLOSC2_DPARAMS_DEFAULTS) @@ -5430,27 +5455,6 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param if rc < 0 or block_nbytes[i] <= 0: return _matmul_exit(input_buffers, dctx, -1) if first_run: - if i == 0: - q = ndarr.blockshape[ndim - 1] - p = ndarr.blockshape[ndim - 2] - # nchunk_ = chunks_in_row * chunk_row + chunk_col - # convert from chunk_idx to element idx chunk_i (row) - chunk_i = nchunk_ // ncols * out_chunk_nrows - chunk_startA = nchunkA + chunk_i // ndarr.chunkshape[ndim - 2] * Ancols - nchunkA = chunk_startA - # nblock_ = blocks_in_chunkrow * block_row + block_col - # convert from block_idx to element idx block_i (row) - block_i = nblock_ // block_ncols * out_block_nrows - block_startA = nblockA + block_i // p * Ablock_ncols - else: # i = 1 - r = ndarr.blockshape[ndim - 1] - # convert from chunk_idx to element idx chunk_j (col) - chunk_j = nchunk_ % ncols * out_chunk_ncols - chunk_startB = nchunkB + chunk_j // ndarr.chunkshape[ndim - 1] - nchunkB = chunk_startB - # convert from block_idx to element idx block_j (col) - block_j = nblock_ % block_ncols * out_block_ncols - block_startB = nblockB + block_j // r input_buffers[i] = malloc(block_nbytes[i]) if input_buffers[i] == NULL: return _matmul_exit(input_buffers, dctx, -1) diff --git a/src/blosc2/linalg.py b/src/blosc2/linalg.py index 3d4c16740..7b21953db 100644 --- a/src/blosc2/linalg.py +++ b/src/blosc2/linalg.py @@ -98,6 +98,14 @@ def _matmul_can_use_fast_path( if not chunks_aligned: return False + # The prefilter maps each output chunk to A's chunk row and B's chunk column and + # walks K chunk by chunk in both, so all three must share one chunk grid + same_chunks = result.chunks[-2] == x1.chunks[-2] + same_chunks &= result.chunks[-1] == x2.chunks[-1] + same_chunks &= x1.chunks[-1] == x2.chunks[-2] + if not same_chunks: + return False + same_blocks = x2.blocks[-2] == x1.blocks[-1] same_blocks &= x2.blocks[-1] == result.blocks[-1] same_blocks &= result.blocks[-2] == x1.blocks[-2] diff --git a/tests/ndarray/test_linalg.py b/tests/ndarray/test_linalg.py index 21269d044..e0fd2de0b 100644 --- a/tests/ndarray/test_linalg.py +++ b/tests/ndarray/test_linalg.py @@ -520,6 +520,79 @@ def test_matmul_prefilter_supported_dtypes(dtype, nthreads): np.testing.assert_array_equal(out[:], np.matmul(a, b)) +def _small_int_operands(mkn, dtype): + # Small integers, so that float32 and float64 products are exact. Not constant: with + # ones()/full() all chunks hold the same values, which hides reading the wrong chunk. + m, k, n = mkn + rng = np.random.default_rng(0) + a = rng.integers(-3, 4, size=(m, k)).astype(dtype) + b = rng.integers(-3, 4, size=(k, n)).astype(dtype) + return a, b + + +@pytest.mark.parametrize( + ("mkn", "achunks", "bchunks"), + [ + ((200, 100, 100), (100, 100), (100, 100)), # 2 chunks along the output rows + ((100, 100, 200), (100, 100), (100, 100)), # 2 chunks along the output columns + ((100, 200, 100), (100, 100), (100, 100)), # 2 chunks along K + ((300, 250, 200), (100, 100), (100, 100)), # all three, partial last chunks + ((200, 300, 150), (100, 150), (150, 50)), # non-square chunks + ], +) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("nthreads", PREFILTER_NTHREADS) +def test_matmul_fast_path_multichunk(monkeypatch, mkn, achunks, bchunks, dtype, nthreads): + old_flag = utils_mod.try_miniexpr + calls = _set_pref_matmul_call_recorder(monkeypatch) + try: + _toggle_miniexpr(True) + a, b = _small_int_operands(mkn, dtype) + cparams = blosc2.CParams(nthreads=nthreads) + x1 = blosc2.asarray(a, chunks=achunks, blocks=(50, 50), cparams=cparams) + x2 = blosc2.asarray(b, chunks=bchunks, blocks=(50, 50), cparams=cparams) + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message=".*encountered in matmul", category=RuntimeWarning) + c = blosc2.matmul(x1, x2, chunks=(achunks[0], bchunks[1]), blocks=(50, 50), cparams=cparams) + expected = a @ b + + assert calls == [(expected.shape, a.shape, b.shape, np.dtype(dtype))] + np.testing.assert_array_equal(c[:], expected) + finally: + _toggle_miniexpr(old_flag) + + +@pytest.mark.parametrize( + ("achunks", "bchunks", "outchunks"), + [ + ((100, 100), (100, 100), (50, 100)), # output chunk rows differ from A's + ((100, 100), (100, 100), (100, 50)), # output chunk columns differ from B's + ((100, 100), (50, 100), (100, 100)), # A and B chunk K differently + ], +) +def test_matmul_fast_path_needs_shared_chunk_grid(monkeypatch, achunks, bchunks, outchunks): + # The prefilter maps output chunks onto A's and B's chunks one to one; other layouts + # gave wrong results, errors or crashes there, so they must take the chunked path + old_flag = utils_mod.try_miniexpr + calls = _set_pref_matmul_call_recorder(monkeypatch) + try: + _toggle_miniexpr(True) + a, b = _small_int_operands((200, 200, 200), np.float64) + x1 = blosc2.asarray(a, chunks=achunks, blocks=(50, 50)) + x2 = blosc2.asarray(b, chunks=bchunks, blocks=(50, 50)) + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message=".*encountered in matmul", category=RuntimeWarning) + c = blosc2.matmul(x1, x2, chunks=outchunks, blocks=(50, 50)) + expected = a @ b + + assert calls == [] + np.testing.assert_array_equal(c[:], expected) + finally: + _toggle_miniexpr(old_flag) + + @pytest.mark.parametrize( ("ashape", "achunks", "ablocks"), [ From 287d160c00c9eb78920ee5e24de31dc51838f864 Mon Sep 17 00:00:00 2001 From: Fabian Isensee Date: Fri, 6 Nov 2026 09:18:48 +0100 Subject: [PATCH 3/3] Keep the operands' padding out of matmul fast-path results The matmul prefilter multiplies whole blocks. When K is not a multiple of the chunk shape, the last blocks along K reach into the chunks' padding, and A's padding columns times B's padding rows were added to the result. The padding is not always zero: blosc2.full()/ones() fill it with the value, and resize() or partial writes into such arrays keep old values there. ones((170, 170)) @ full((170, 170), 3) with 200x200 chunks and 50x50 blocks was off by 90 in every element (main and 4.14.1), and NaN in the padding turned every element into NaN. Stop at the first block that starts past K (only padding follows), and zero A's padding columns and B's padding rows in a block that ends past K. Both are zeroed so that NaN or inf in one operand's padding cannot leak through 0 * NaN. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01DiPFzsfE5i74cBSxZN13M8 --- src/blosc2/blosc2_ext.pyx | 16 ++++++++++++++ tests/ndarray/test_linalg.py | 43 ++++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+) diff --git a/src/blosc2/blosc2_ext.pyx b/src/blosc2/blosc2_ext.pyx index baf0a5f89..b2ae169df 100644 --- a/src/blosc2/blosc2_ext.pyx +++ b/src/blosc2/blosc2_ext.pyx @@ -5428,6 +5428,9 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param # convert from block_idx to element idx block_i (row) block_i = nblock_ // block_ncols * out_block_nrows block_startA = nblockA + block_i // p * Ablock_ncols + # Elements of K from the current block on. The chunks along K are whole blocks (see + # _matmul_can_use_fast_path()), so each step of the block loop advances K by q. + cdef int64_t k_left = ndarr.shape[ndim - 1] ndarr = udata.inputs[1] ndim = ndarr.ndim r = ndarr.blockshape[ndim - 1] @@ -5463,6 +5466,8 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param nblockA = block_startA nblockB = block_startB while True: # block loop + if k_left <= 0: + break # only padding left (in the last chunk along K) startA = nblockA * block_nbytes[0] startB = nblockB * block_nbytes[1] # In bytes, for the reason given in aux_miniexpr(). matmul only ever @@ -5478,6 +5483,16 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param if rc != block_nbytes[1]: return _matmul_exit(input_buffers, dctx, -1) batch = 0 + if k_left < q: + # The block ends in the padding, which need not hold zeros (full() fills it, + # resize() leaves old data there): zero A's padding columns and B's padding + # rows so that their products add nothing + for i in range(block_nbytes[0] // (q * typesize)): + memset(input_buffers[0] + (i * q + k_left) * typesize, 0, + (q - k_left) * typesize) + for i in range(block_nbytes[1] // (q * r * typesize)): + memset(input_buffers[1] + (i * q + k_left) * r * typesize, 0, + (q - k_left) * r * typesize) while batch < batches: batch_ = batch offsetA = 0 @@ -5554,6 +5569,7 @@ cdef int aux_matmul(mm_udata *udata, int64_t nchunk, int32_t nblock, void *param else: return _matmul_exit(input_buffers, dctx, -1) batch += 1 + k_left -= q nblockA += 1 nblockB += Bblock_ncols if (nblockA % Ablock_ncols == 0): diff --git a/tests/ndarray/test_linalg.py b/tests/ndarray/test_linalg.py index e0fd2de0b..e3e3f83cc 100644 --- a/tests/ndarray/test_linalg.py +++ b/tests/ndarray/test_linalg.py @@ -563,6 +563,49 @@ def test_matmul_fast_path_multichunk(monkeypatch, mkn, achunks, bchunks, dtype, _toggle_miniexpr(old_flag) +def _nan_padded(data, chunks, cparams): + # Shrinking leaves the cut values in the chunks' padding + big = np.full((data.shape[0] + 60, data.shape[1] + 60), np.nan, dtype=data.dtype) + big[: data.shape[0], : data.shape[1]] = data + arr = blosc2.asarray(big, chunks=chunks, blocks=(50, 50), cparams=cparams) + arr.resize(data.shape) + return arr + + +@pytest.mark.parametrize( + ("mkn", "chunks"), + [ + ((130, 170, 90), (100, 100)), # 2 chunks along K, the last one ends inside a block + ((170, 170, 170), (200, 200)), # 1 chunk along K, ends inside a block + ((100, 100, 100), (200, 200)), # the last 2 blocks along K are only padding + ((120, 120, 120), (200, 200)), # both + ], +) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("nthreads", PREFILTER_NTHREADS) +def test_matmul_fast_path_ignores_padding(monkeypatch, mkn, chunks, dtype, nthreads): + # The prefilter multiplies whole blocks, so along K it also reads the chunks' padding, + # which need not hold zeros (full() fills it as well) + old_flag = utils_mod.try_miniexpr + calls = _set_pref_matmul_call_recorder(monkeypatch) + try: + _toggle_miniexpr(True) + a, b = _small_int_operands(mkn, dtype) + cparams = blosc2.CParams(nthreads=nthreads) + x1 = _nan_padded(a, chunks, cparams) + x2 = _nan_padded(b, chunks, cparams) + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message=".*encountered in matmul", category=RuntimeWarning) + c = blosc2.matmul(x1, x2, chunks=chunks, blocks=(50, 50), cparams=cparams) + expected = a @ b + + assert calls == [(expected.shape, a.shape, b.shape, np.dtype(dtype))] + np.testing.assert_array_equal(c[:], expected) + finally: + _toggle_miniexpr(old_flag) + + @pytest.mark.parametrize( ("achunks", "bchunks", "outchunks"), [