Hacktoberfest 2026：メンテナが10月に向けて印を付けた、オープンで初心者向けの issue。 Hacktoberfest の issue を見る

Batched GEMM: missing small-kernel guards, Fortran argument swaps, and empty-batch crash

クローズ
  • 何か問題がありますか？
#6,121 コメント 0 件 リアクション 0 件 担当者 0 名 GitHub で見る

メンテナーはふだん 1 日以内に返信

まだ誰も着手していません。

評価

難易度
4/5
見積もり時間
3〜5日
初心者へのやさしさ
32/100
issue の種類
バグ
明瞭さ
おおむね明確
活発さ
活発
技術スタック
c, fortran
領域
backend

調査の方向性

Four separate defects are listed across interface/gemm_batch.c, interface/gemm_batch_strided.c and kernel/generic/gemm_small_matrix_kernel_nn.c. Start by comparing the small-kernel dispatch in those interfaces with the alpha and permit checks in the ordinary GEMM path, then run the reproducer's check_a through check_d against a patched library. Done means all four checks return PASS, and the reproducer's exit status is 0.

索引モデルが issue の本文から書いたものです。

説明

On develop 33746423, the batched GEMM interfaces differ from the corresponding ordinary GEMM operations in the four cases below. The original witnesses were reproduced on an x86-64 HASWELL build, natively on Apple M5, and under QEMU on arm64. The numerical cases use one matrix per batch; the zero-dimension case uses a two-member strided batch.

A. The alpha-zero guard is missing from batched small-matrix dispatch

The guard added to ordinary GEMM by #6083/#6085 is still missing from cblas_?gemm_batch and cblas_?gemm_batch_strided. With alpha=0, beta=1, C initialized to 1 and NaNs in A, both double-precision batch forms produce NaNs; the ordinary cblas_dgemm control leaves C unchanged.

gemm_batch.c:351 and gemm_batch_strided.c:364 select small kernels using only M·N·K <= 100³, without ordinary GEMM's alpha/permit checks. This shared dispatch affects the real and complex variants.

B. The Fortran interfaces swap operands or transpose flags

For column-major A = [1 3; 2 4] and B = [5 7; 6 8]:

Call Ordinary GEMM control, in memory order Batched result
dgemm_batch_, N,N 23 34 31 46 19 22 43 50 (B·A)
dgemm_batch_strided_, N,T 26 38 30 44 17 39 23 53 (Aᵀ·B)

The grouped Fortran branch exchanges TRANSA/TRANSB and passes b_array as A and a_array as B, without exchanging M/N. The strided Fortran branch keeps A/B in place but exchanges their transpose flags.

C. Zero-sized strided members leave uninitialized entries to dispatch

cblas_dgemm_batch_strided with M=0, N=K=2 and count=2 segfaults in the original witness. There are no C entries to update.

The zero-dimension branch skips initialization with continue, but the final dispatch still passes the original count. N=0 takes the same branch. The grouped interface instead maintains a compacted count.

D. SBGEMM batch reaches a small kernel that multiplies BF16 encodings

For BF16 inputs containing the same A and B, alpha=1 and beta=0, ordinary cblas_sbgemm returns 23 34 31 46. The strided batch returns approximately 5.41581e+08 5.44760e+08 5.43674e+08 5.46865e+08.

This shares the missing permit check in A. Ordinary SBGEMM's generic small-matrix permit is zero, but the batch interfaces bypass it and reach the generic multiplication with IFLOAT equal to uint16_t BF16 storage.

Reproducer

Linux, LP64 OpenBLAS (32-bit BLAS integers), with BF16 support enabled:

gcc -O2 -Wall -Wextra -Werror repro.c -ldl -lm -o repro
./repro /absolute/path/to/libopenblas.so

The checker below uses four-element matrix strides and checks each result against its exact expected value. Its exit status is 0 only when all four checks pass; mismatches or signals return 1, and missing symbols, failed controls or an incompatible integer configuration return 2. It sets one OpenBLAS thread.

The platform observations above come from the original witnesses. This revised checker adds explicit verdicts and conventional positive strides; it has not yet been rerun against the patched OpenBLAS library.

repro.c 
/* Exact expected-result checks for batched GEMM on LP64 OpenBLAS builds.
 * Derived from the review packet's repro.c. A failed or missing control, a
 * fault, and an unchanged wrong output are failures, never evidence of a fix.
 * Usage: repro LIBRARY [CORETYPE]
 * All checks run with OPENBLAS_NUM_THREADS=1. This does not test parallel use.
 */
#define _GNU_SOURCE
#include <dlfcn.h>
#include <errno.h>
#include <math.h>
#include <signal.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/resource.h>
#include <sys/wait.h>
#include <unistd.h>

enum { Col = 102, N_ = 111, T_ = 112 };
enum { PASS = 0, FAIL = 1, ERROR = 2 };
typedef uint16_t bf16;
static void *lib;
static void *need(const char *name) {
    void *p = dlsym(lib, name);
    if (!p) { printf("MISSING %s\n", name); exit(ERROR); }
    return p;
}
typedef void (*dgemm_c)(int, int, int, int, int, int, double, const double *, int, const double *, int, double, double *, int);
typedef void (*dstr_c)(int, int, int, int, int, int, double, const double *, int, int, const double *, int, int, double, double *, int, int, int);
typedef void (*dbat_c)(int, const int *, const int *, const int *, const int *, const int *, const double *, const double **, const int *, const double **, const int *, const double *, double **, const int *, int, const int *);
typedef void (*dbat_f)(const char *, const char *, const int *, const int *, const int *, const double *, const double **, const int *, const double **, const int *, const double *, double **, const int *, const int *, const int *);
typedef void (*dstr_f)(const char *, const char *, const int *, const int *, const int *, const double *, const double *, const int *, const int *, const double *, const int *, const int *, const double *, double *, const int *, const int *, const int *);
typedef void (*sbgemm_c)(int, int, int, int, int, int, float, const bf16 *, int, const bf16 *, int, float, float *, int);
typedef void (*sbstr_c)(int, int, int, int, int, int, float, const bf16 *, int, int, const bf16 *, int, int, float, float *, int, int, int);

static void show(const char *label, const double *c) { printf("  %-34s C = %g %g %g %g\n", label, c[0], c[1], c[2], c[3]); }
static int same4(const double *x, const double *y) { return !memcmp(x, y, 4 * sizeof(double)); }
static int control_ok(const char *label, const double *c, const double *want) {
    if (same4(c, want)) return 1;
    printf("  CONTROL-WRONG %s: expected %g %g %g %g\n", label, want[0], want[1], want[2], want[3]);
    return 0;
}
static int check_a(void) {
    double a[4] = {NAN, NAN, NAN, NAN}, b[4] = {1, 1, 1, 1};
    double c0[4] = {1, 1, 1, 1}, c1[4] = {1, 1, 1, 1}, c2[4] = {1, 1, 1, 1};
    ((dgemm_c)need("cblas_dgemm"))(Col, N_, N_, 2, 2, 2, 0.0, a, 2, b, 2, 1.0, c0, 2);
    ((dstr_c)need("cblas_dgemm_batch_strided"))(Col, N_, N_, 2, 2, 2, 0.0, a, 2, 4, b, 2, 4, 1.0, c1, 2, 4, 1);
    int ta = N_, tb = N_, m = 2, n = 2, k = 2, ld = 2, gs = 1; double al = 0, be = 1;
    const double *ap[1] = {a}, *bp[1] = {b}; double *cp[1] = {c2};
    ((dbat_c)need("cblas_dgemm_batch"))(Col, &ta, &tb, &m, &n, &k, &al, ap, &ld, bp, &ld, &be, cp, &ld, 1, &gs);
    printf("A alpha=0, beta=1, A=NaN, C=1 (expected C unchanged)\n");
    show("cblas_dgemm (control)", c0); show("cblas_dgemm_batch_strided", c1); show("cblas_dgemm_batch", c2);
    const double ones[4] = {1, 1, 1, 1};
    if (!control_ok("cblas_dgemm", c0, ones)) return ERROR;
    return same4(c1, ones) && same4(c2, ones) ? PASS : FAIL;
}
static int check_b(void) {
    double a[4] = {1, 2, 3, 4}, b[4] = {5, 6, 7, 8}, c0[4] = {0}, c1[4] = {-1, -1, -1, -1}, c2[4] = {0}, c3[4] = {-1, -1, -1, -1};
    int m = 2, n = 2, k = 2, ld = 2, one = 1, s4 = 4; double al = 1, be = 0;
    dgemm_c g = (dgemm_c)need("cblas_dgemm");
    g(Col, N_, N_, 2, 2, 2, 1.0, a, 2, b, 2, 0.0, c0, 2);
    const double *ap[1] = {a}, *bp[1] = {b}; double *cp[1] = {c1};
    ((dbat_f)need("dgemm_batch_"))("N", "N", &m, &n, &k, &al, ap, &ld, bp, &ld, &be, cp, &ld, &one, &one);
    g(Col, N_, T_, 2, 2, 2, 1.0, a, 2, b, 2, 0.0, c2, 2);
    ((dstr_f)need("dgemm_batch_strided_"))("N", "T", &m, &n, &k, &al, a, &ld, &s4, b, &ld, &s4, &be, c3, &ld, &s4, &one);
    printf("B Fortran interfaces, A = [1 3; 2 4], B = [5 7; 6 8]\n");
    show("cblas_dgemm N,N (control)", c0); show("dgemm_batch_ N,N", c1);
    show("cblas_dgemm N,T (control)", c2); show("dgemm_batch_strided_ N,T", c3);
    const double nn[4] = {23, 34, 31, 46}, nt[4] = {26, 38, 30, 44};
    int ok0 = control_ok("cblas_dgemm N,N", c0, nn), ok2 = control_ok("cblas_dgemm N,T", c2, nt);
    if (!ok0 || !ok2) return ERROR;
    return same4(c1, nn) && same4(c3, nt) ? PASS : FAIL;
}
static int check_c(void) {
    double a[8] = {1, 1, 1, 1, 1, 1, 1, 1}, b[8] = {1, 1, 1, 1, 1, 1, 1, 1};
    double c0[8] = {5, 5, 5, 5, 5, 5, 5, 5};
    double c1[8] = {5, 5, 5, 5, 5, 5, 5, 5};
    printf("C cblas_dgemm_batch_strided M=0, N=K=2, count=2 (expected: return, no write)\n");
    ((dgemm_c)need("cblas_dgemm"))(Col, N_, N_, 0, 2, 2, 1.0, a, 2, b, 2, 1.0, c0, 2);
    for (int i = 0; i < 8; i++) if (c0[i] != 5) {
        printf("  CONTROL-WRONG cblas_dgemm M=0 at C[%d]\n", i); return ERROR;
    }
    ((dstr_c)need("cblas_dgemm_batch_strided"))(Col, N_, N_, 0, 2, 2, 1.0, a, 2, 4, b, 2, 4, 1.0, c1, 2, 4, 2);
    for (int i = 0; i < 8; i++) if (c1[i] != 5) {
        printf("  WRONG cblas_dgemm_batch_strided M=0: C[%d]=%g, want 5\n", i, c1[i]); return FAIL;
    }
    printf("  returned, all 8 C elements unchanged\n");
    return PASS;
}
static bf16 tobf(float f) { uint32_t u; memcpy(&u, &f, 4); return (bf16)(u >> 16); }
static int check_d(void) {
    bf16 a[4], b[4]; float c0[4] = {0}, c1[4] = {-1, -1, -1, -1};
    for (int i = 0; i < 4; i++) { a[i] = tobf((float)(i + 1)); b[i] = tobf((float)(i + 5)); }
    ((sbgemm_c)need("cblas_sbgemm"))(Col, N_, N_, 2, 2, 2, 1.0f, a, 2, b, 2, 0.0f, c0, 2);
    ((sbstr_c)need("cblas_sbgemm_batch_strided"))(Col, N_, N_, 2, 2, 2, 1.0f, a, 2, 4, b, 2, 4, 0.0f, c1, 2, 4, 1);
    printf("D SBGEMM, BF16 A = [1 3; 2 4], B = [5 7; 6 8], alpha=1, beta=0\n");
    printf("  %-34s C = %g %g %g %g\n", "cblas_sbgemm (control)", c0[0], c0[1], c0[2], c0[3]);
    printf("  %-34s C = %g %g %g %g\n", "cblas_sbgemm_batch_strided", c1[0], c1[1], c1[2], c1[3]);
    const float want[4] = {23, 34, 31, 46};
    if (memcmp(c0, want, sizeof want)) { printf("  CONTROL-WRONG cblas_sbgemm\n"); return ERROR; }
    return !memcmp(c1, want, sizeof want) ? PASS : FAIL;
}
int main(int argc, char **argv) {
    if (argc < 2) { fprintf(stderr, "usage: %s LIBRARY [CORETYPE]\n", argv[0]); return ERROR; }
    struct rlimit no_core = {0, 0}; setrlimit(RLIMIT_CORE, &no_core);
    setvbuf(stdout, NULL, _IOLBF, 0);
    if (argc > 2 && setenv("OPENBLAS_CORETYPE", argv[2], 1)) { perror("setenv"); return ERROR; }
    if (setenv("OPENBLAS_NUM_THREADS", "1", 1)) { perror("setenv"); return ERROR; }
    lib = dlopen(argv[1], RTLD_NOW | RTLD_LOCAL);
    if (!lib) { fprintf(stderr, "dlopen: %s\n", dlerror()); return ERROR; }
    const char *(*get_config)(void) = (const char *(*)(void))dlsym(lib, "openblas_get_config");
    if (get_config) {
        const char *config = get_config();
        printf("LIBRARY %s\nCONFIG %s\n", argv[1], config ? config : "unavailable");
        if (config && strstr(config, "USE64BITINT")) {
            fprintf(stderr, "ERROR: this reproducer requires LP64 OpenBLAS (32-bit BLAS integers).\n");
            return ERROR;
        }
    } else {
        printf("LIBRARY %s\nCONFIG unavailable; requires LP64 OpenBLAS (32-bit BLAS integers)\n", argv[1]);
    }
    int (*checks[4])(void) = {check_a, check_b, check_c, check_d};
    int passed = 0, failed = 0, errors = 0;
    for (int i = 0; i < 4; i++) {
        fflush(stdout);
        pid_t pid = fork();
        if (pid < 0) { perror("fork"); printf("CHECK %c ERROR fork\n", 'A' + i); errors++; continue; }
        if (pid == 0) { alarm(20); int rc = checks[i](); fflush(stdout); _exit(rc); }
        int st = 0; pid_t got;
        do { got = waitpid(pid, &st, 0); } while (got < 0 && errno == EINTR);
        if (got != pid) { perror("waitpid"); printf("CHECK %c ERROR waitpid\n", 'A' + i); errors++; }
        else if (WIFSIGNALED(st)) { printf("CHECK %c FAIL signal=%d\n", 'A' + i, WTERMSIG(st)); failed++; }
        else if (!WIFEXITED(st)) { printf("CHECK %c ERROR child status=%d\n", 'A' + i, st); errors++; }
        else if (WEXITSTATUS(st) == PASS) { printf("CHECK %c PASS\n", 'A' + i); passed++; }
        else if (WEXITSTATUS(st) == FAIL) { printf("CHECK %c FAIL\n", 'A' + i); failed++; }
        else { printf("CHECK %c ERROR child-exit=%d\n", 'A' + i, WEXITSTATUS(st)); errors++; }
    }
    printf("SUMMARY: %s passed=%d/4 failed=%d errors=%d core=%s threads=1\n",
           errors ? "ERROR" : failed ? "FAIL" : "PASS", passed, failed, errors, argc > 2 ? argv[2] : "dispatch");
    return errors ? ERROR : failed ? FAIL : PASS;
}
Candidate fix

The two-file diff below reuses ordinary GEMM's alpha/permit condition, corrects the Fortran argument mappings, and returns for zero-sized strided members. I am offering it as a candidate fix; I can prepare a PR with regression tests if useful.

batch-gemm.patch 
diff --git a/interface/gemm_batch.c b/interface/gemm_batch.c
index 7a8789e..e36b861 100644
--- a/interface/gemm_batch.c
+++ b/interface/gemm_batch.c
@@ -298,19 +298,19 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE *  transa_array, enum CB
       group_ldb = ldb_array[i];
       group_ldc = ldc_array[i];
       
-      if (transb_array[i] == 'N')   group_transa = 0;
-      if (transb_array[i] == 'T')   group_transa = 1;
+      if (transa_array[i] == 'N')   group_transa = 0;
+      if (transa_array[i] == 'T')   group_transa = 1;
 #ifndef COMPLEX
-      if (transb_array[i] == 'C')   group_transa = 1;
+      if (transa_array[i] == 'C')   group_transa = 1;
 #else
-      if (transb_array[i] == 'C')   group_transa = 3;
+      if (transa_array[i] == 'C')   group_transa = 3;
 #endif
-      if (transa_array[i] == 'N')   group_transb = 0;
-      if (transa_array[i] == 'T')   group_transb = 1;
+      if (transb_array[i] == 'N')   group_transb = 0;
+      if (transb_array[i] == 'T')   group_transb = 1;
 #ifndef COMPLEX
-      if (transa_array[i] == 'C')   group_transb = 1;
+      if (transb_array[i] == 'C')   group_transb = 1;
 #else
-      if (transa_array[i] == 'C')   group_transb = 3;
+      if (transb_array[i] == 'C')   group_transb = 3;
 #endif
 #endif
 
@@ -348,7 +348,16 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE *  transa_array, enum CB
 #endif
 
 #ifdef SMALL_MATRIX_OPT
-    if (MNK <= 100.0*100.0*100.0){
+#if !defined(COMPLEX)
+    if (*(FLOAT *)(group_alpha) != ZERO &&
+        GEMM_SMALL_MATRIX_PERMIT(group_transa, group_transb, group_m, group_n, group_k,
+                                 *(FLOAT *)(group_alpha), *(FLOAT *)(group_beta))){
+#else
+    if (((FLOAT *)(group_alpha))[0] != ZERO &&
+        GEMM_SMALL_MATRIX_PERMIT(group_transa, group_transb, group_m, group_n, group_k,
+                                 ((FLOAT *)(group_alpha))[0], ((FLOAT *)(group_alpha))[1],
+                                 ((FLOAT *)(group_beta))[0], ((FLOAT *)(group_beta))[1])){
+#endif
       group_routine=NULL;
 #if !defined(COMPLEX)
       if(*(FLOAT *)(group_beta) == 0.0){
@@ -388,16 +397,15 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE *  transa_array, enum CB
       args_array[count].beta=group_beta;
 
 #if defined(CBLAS)
-      if (order == CblasColMajor) {      
-	args_array[count].a=(a_array[matrix_idx+j]);
-	args_array[count].b=(b_array[matrix_idx+j]);
-      }else if(order == CblasRowMajor){
-#endif
+      if (order == CblasRowMajor) {
 	args_array[count].a=(b_array[matrix_idx+j]);
 	args_array[count].b=(a_array[matrix_idx+j]);
-#if defined(CBLAS)
-      }
+      } else
 #endif
+      {
+	args_array[count].a=(a_array[matrix_idx+j]);
+	args_array[count].b=(b_array[matrix_idx+j]);
+      }
       
       args_array[count].c=(c_array[matrix_idx+j]);
       
diff --git a/interface/gemm_batch_strided.c b/interface/gemm_batch_strided.c
index f711b77..841933e 100644
--- a/interface/gemm_batch_strided.c
+++ b/interface/gemm_batch_strided.c
@@ -311,19 +311,19 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE transa, enum CBLAS_TRANS
       group_ldb = ldb;
       group_ldc = ldc;
       
-      if (tb == 'N')   group_transa = 0;
-      if (tb == 'T')   group_transa = 1;
+      if (ta == 'N')   group_transa = 0;
+      if (ta == 'T')   group_transa = 1;
 #ifndef COMPLEX
-      if (tb == 'C')   group_transa = 1;
+      if (ta == 'C')   group_transa = 1;
 #else
-      if (tb == 'C')   group_transa = 3;
+      if (ta == 'C')   group_transa = 3;
 #endif
-      if (ta == 'N')   group_transb = 0;
-      if (ta == 'T')   group_transb = 1;
+      if (tb == 'N')   group_transb = 0;
+      if (tb == 'T')   group_transb = 1;
 #ifndef COMPLEX
-      if (ta == 'C')   group_transb = 1;
+      if (tb == 'C')   group_transb = 1;
 #else
-      if (ta == 'C')   group_transb = 3;
+      if (tb == 'C')   group_transb = 3;
 #endif
 #endif
 
@@ -352,7 +352,8 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE transa, enum CBLAS_TRANS
       return;
     }
 
-    if (group_m == 0 || group_n == 0) continue;
+    /* Every member shares M and N: an empty member means an empty batch. */
+    if (group_m == 0 || group_n == 0) { free(args_array); return; }
 
     group_mode=mode;
 
@@ -361,7 +362,16 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE transa, enum CBLAS_TRANS
 #endif
 
 #ifdef SMALL_MATRIX_OPT
-    if (MNK <= 100.0*100.0*100.0){
+#if !defined(COMPLEX)
+    if (*(FLOAT *)(&alpha) != ZERO &&
+        GEMM_SMALL_MATRIX_PERMIT(group_transa, group_transb, group_m, group_n, group_k,
+                                 *(FLOAT *)(&alpha), *(FLOAT *)(&beta))){
+#else
+    if (((FLOAT *)(alpha))[0] != ZERO &&
+        GEMM_SMALL_MATRIX_PERMIT(group_transa, group_transb, group_m, group_n, group_k,
+                                 ((FLOAT *)(alpha))[0], ((FLOAT *)(alpha))[1],
+                                 ((FLOAT *)(beta))[0], ((FLOAT *)(beta))[1])){
+#endif
       group_routine=NULL;
 #if !defined(COMPLEX)
       if(beta == 0.0){

Found through AI-assisted code review and differential testing (Claude); the report and reproducer were also reviewed with GPT-6 Astra.

主要言語
C
スター
7.6k
フォーク
1.7k
平均マージ
1日 11時間
マージ済み PR（30日）
55

環境構築

このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

OpenMathLib/OpenBLAS のほかの issue

OpenMathLib/OpenBLAS の issue をすべて見る

似ている issue

C の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。