trmm: Inf in B becomes NaN through the unreferenced triangle of A (0*Inf in diagonal micro-panels) #967

Reproducer

Take the 2×2 identity matrix $A$ and the vector $B$:

Mathematically, we would expect

With dtrmm('L', 'L', 'N', 'N', alpha = 1), which computes B := A*B with
the lower triangle of A, Reference BLAS produces [1; Inf], whereas BLIS
returns [NaN; Inf]. A(1,2) lies in the strictly upper triangle and is
not referenced for uplo = 'L'. We reproduced this on an Apple M5, AMD
Zen 4, Intel Ice Lake and Intel Alder Lake; the precise results are
documented below. Whether the BLAS interface requires the Reference BLAS
behavior for IEEE-754 special values is still an open question, see
Reference-LAPACK discussion #1469.

The minimal C reproducer:

#include <math.h>
#include <stdio.h>

void dtrmm_(const char *side, const char *uplo, const char *transa,
            const char *diag, const int *m, const int *n, const double *alpha,
            const double *A, const int *lda, double *B, const int *ldb);

int main(void)
{
    int m = 2, n = 1;
    double alpha = 1;
    double A[] = { 1, 0,     /* column 1 */
                   0, 1 };   /* column 2 */
    double B[] = { 1, INFINITY };

    printf("A = [%g %g; %g %g]\n", A[0], A[2], A[1], A[3]);
    printf("B = [%g; %g]\n", B[0], B[1]);
    dtrmm_("L", "L", "N", "N", &m, &n, &alpha, A, &m, B, &m);
    printf("dtrmm('L', 'L', 'N', 'N'): B = [%g; %g]\n", B[0], B[1]);
}

./configure --prefix=$PREFIX auto && make -j && make install   # BLIS
cc tiny.c -o tiny -L$PREFIX/lib -lblis -Wl,-rpath,$PREFIX/lib && ./tiny

Reference BLAS (LAPACK 3.12.1)             BLIS master (a8037cfc)
A = [1 0; 0 1]                             A = [1 0; 0 1]
B = [1; inf]                               B = [1; inf]
dtrmm('L', 'L', 'N', 'N'): B = [1; inf]    dtrmm('L', 'L', 'N', 'N'): B = [nan; inf]

BLIS evaluates A(1,2)*B(2) = 0*Inf = NaN, although A(1,2) is not
referenced. strmm behaves the same way.

Results on x86-64

The same 2×2 program (tiny.c above, unchanged), linked statically or
directly against each library, sequential. BLIS master was built with
./configure auto on each node; the configuration reported is
bli_arch_string(bli_arch_query_id()) at run time.

Compilers: GCC 14.2.0 (EPYC, Xeon), GCC 12.2.0 (Core i5). On the
EPYC 9454, BLIS master selects zen3, since it has no zen4
configuration.

OpenBLAS also returns [NaN; Inf] for this example (0.3.21 on the
Core i5-12500, 0.3.34 on the EPYC 9454, the Xeon 8358 and the Apple M5).

Summary

?trmm computes 0*Inf = NaN for elements of the triangle of A that is
not referenced, so an Inf in B turns other entries of the result into
NaN. The reference BLAS multiplies only by elements of the referenced
triangle. BLIS does not read the unreferenced triangle (NaN stored there
does not leak into the result); the zeros come from packing (see Cause
below).

Background: Reference-LAPACK/lapack#1469

Systematic test

A test over all combinations of side, uplo, transa and diag, s and d,
sizes 1 to 100 of the triangular dimension (the other dimension 3),
A = identity, one Inf in B at every position of the triangular
dimension (17216 calls). The unreferenced triangle of A holds 0 in one
run and NaN in another. A result entry counts as wrong if it is NaN where
a loop over the referenced elements of op(A) gives no NaN.

NaN that comes from explicitly stored zeros inside the referenced
triangle is not counted. The reference BLAS itself computes 0*Inf for
those when side = 'L' (it skips zeros of B, not of A), but skips zeros
of A when side = 'R'. So that case is implementation-defined, and this
issue is only about the unreferenced triangle.

All 16 combinations of side, uplo, transa and diag are affected, in both
precisions. The results are identical whether the unreferenced triangle
holds 0 or NaN.

The spurious NaN stay within one register block of the Inf:

dtrmm('L', 'U', 'N', 'N'), m = 24, n = 1, A = identity, Inf in row p (0-based)

           reference BLAS              BLIS firestorm             BLIS generic
p = 0:     -                           1 2 3 4 5                  1 2 3 4 5 6 7
p = 1:     0                           0 2 3 4 5                  0 2 3 4 5 6 7
p = 6:     0 1 2 3 4 5                 0 1 2 3 4 5 7 8 9 10 11    0 1 2 3 4 5 7

(NaN rows listed. The NaN rows of the reference BLAS come from the
explicitly stored zeros above the diagonal, which are part of the
referenced triangle for uplo = 'U'.)

Cause

For micro-panels that intersect the diagonal, the triangular operand is
packed with explicit zeros in its unstored part:

• ref_kernels/1m/bli_packm_cxc_diag_ref.c zeroes the whole diagonal block
  ("start by zeroing out the whole block") and then writes the stored
  triangle;
• frame/1m/packm/bli_packm_struc_cxk.c writes zeros (f_setv) for the
  unstored parts p10/p12 of a triangular matrix.

The trmm macro-kernels (frame/3/trmm/bli_trmm_{ll,lu,rl,ru}_ker_var2.c)
then call the regular gemm micro-kernel over the full diagonal block
(k_a1011 = min(diagoffa_i + MR, k) in bli_trmm_ll_ker_var2.c,
k_b0111 in bli_trmm_ru_ker_var2.c), which evaluates
0 * B(k,j) for these zeros.

On the M5, dtrmm('L', 'L', ...) runs bli_trmm_ru_ker_var2 (the
firestorm micro-kernel bli_dgemm_armv8a_asm_8x6r prefers rows, so the
operation is transposed), and the NaN appear in blocks of NR = 6 rows; with
generic they appear in blocks of 8 rows.

A fix would have to avoid multiplying the zeros of the unstored part of the
diagonal blocks with B, for example with a micro-kernel (or an edge path)
for the diagonal blocks that only accumulates over the stored triangle.

Environment

• Apple M5 (MacBook Air), macOS 26.6 (25G72)
• Homebrew clang 23.1.4
• BLIS master a8037cf, ./configure auto
  (firestorm) and ./configure generic, default options
• Reference BLAS: LAPACK 3.12.1 (Homebrew lapack 3.12.1_1)