Fix: CBLAS order handling in gemm_batch_strided for column-major layout

[ajz34]

Hi devs!

This PR will fix calling gemm_batch_strided in CBLAS with column-major layout.
Current code will not set batched-group stride (group_lda/b/c) if column-major (where they were setted if row-major, refer to this code). Consequently, the program will raise error like this:

 ** On entry to DGEMM_BATCH_STRIDED  parameter number  8 had an illegal value

This is due to the group_lda/b/c initialized to zero (but not set proper value later), causing this code returns error code 8.

The fix is strightford: initialize group_lda/b/c like what row-major's code does.

---

This was found when I'm developing rust's numpy-like toolkit RSTSR, when trying to use batched gemm for broadcasted matrix-multiplication implementation with OpenBLAS backend.
This bug was originally found in AI code agent session (Claude Code with model GLM-5.3).

[ajz34]

This comment presents the reproduction of MWE.

Note that this MWE must use a slightly large matrix (to bypass SMALL_MATRIX_OPT, where small matrix batched gemm is probably not correct, as presented in #5980).

The C code to reproduce the MWE.

/*
 * MWE: cblas_dgemm_batch_strided rejects every CblasColMajor call.
 * Notes, root cause and expected output: see bug1_strided_colmajor_large.md
 *
 * Build: gcc bug1_strided_colmajor_large.c -o bug1_large \
 *            -I<openblas>/build/generated -I<openblas>/build \
 *            -Wl,-rpath,<openblas>/build/lib -L<openblas>/build/lib -lopenblas
 *
 * Co-authored-by: Claude Code <noreply@anthropic.com>
 * Co-authored-by: glm-5.3 <service@zhipuai.cn>
 */
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>

/* M*N*K > 100^3 keeps the small-kernel dispatch path out of the picture */
#define M 101
#define N 103
#define K 97
#define BATCH 3
#define SA (M * K) /* packed col-major: stride = matrix footprint */
#define SB (K * N)
#define SC (M * N)

int main(void)
{
    double *a = malloc(sizeof(double) * BATCH * SA);
    double *b = malloc(sizeof(double) * BATCH * SB);
    double *c = calloc(BATCH * SC, sizeof(double));
    double *ref = malloc(sizeof(double) * BATCH * SC);
    if (!a || !b || !c || !ref)
        return 2;

    /* distinct integer data per batch, exact in double precision */
    for (int i = 0; i < BATCH * SA; i++)
        a[i] = (i * 7 + 3) % 19;
    for (int i = 0; i < BATCH * SB; i++)
        b[i] = (i * 5 + 1) % 23;

    cblas_dgemm_batch_strided(CblasColMajor, CblasNoTrans, CblasNoTrans,
                              M, N, K, 1.0,
                              a, M, SA, b, K, SB, 0.0, c, M, SC, BATCH);

    long long checksum = 0;
    int bad = 0;
    for (int t = 0; t < BATCH; t++) {
        for (int j = 0; j < N; j++)
            for (int i = 0; i < M; i++) {
                double s = 0;
                for (int l = 0; l < K; l++)
                    s += a[t * SA + l * M + i] * b[t * SB + j * K + l];
                ref[t * SC + j * M + i] = s;
                checksum += (long long)s;
                if (c[t * SC + j * M + i] != s)
                    bad++;
            }
        printf("batch %d: c[0]=%.1f\n", t, c[t * SC]);
    }
    printf("checksum=%lld  %s\n", checksum, bad ? "FAIL" : "OK");
    return bad != 0;
}

The output from current code (commit e0cabe9)

 ** On entry to DGEMM_BATCH_STRIDED  parameter number  8 had an illegal value
batch 0: c[0]=0.0
batch 1: c[0]=0.0
batch 2: c[0]=0.0
checksum=299684775  FAIL

The expected output (this PR)

batch 0: c[0]=9868.0
batch 1: c[0]=9577.0
batch 2: c[0]=9304.0
checksum=299684775  OK

The Python equilvant.

"""
Rebuilds the same packed col-major batch buffers, runs the batched product
with stacked matmul, and prints the same c[0]s and checksum. NumPy is
row-major, so a col-major MxK matrix in a flat buffer is read as a row-major
(K,M) matrix: reshape then swap the last two axes.
"""
import numpy as np

M, N, K, BATCH = 101, 103, 97, 3
SA, SB, SC = M * K, K * N, M * N

a = ((np.arange(BATCH * SA) * 7 + 3) % 19).astype(np.float64)
b = ((np.arange(BATCH * SB) * 5 + 1) % 23).astype(np.float64)

A = a.reshape(BATCH, K, M).swapaxes(-1, -2)  # (B, M, K)
B = b.reshape(BATCH, N, K).swapaxes(-1, -2)  # (B, K, N)

C = A @ B  # batched GEMM

# flatten each (M, N) result back to col-major order
c = np.ascontiguousarray(C.swapaxes(-1, -2)).reshape(-1)

for t in range(BATCH):
    print(f"batch {t}: c[0]={c[t * SC]:.1f}")
print(f"checksum={c.sum():.0f}")

[martin-frbg]

Thank you. Not sure how that happened, but I guess I should have made sure I tested both branches