Question about REAL_GEMM_R macro in potrf_L_single.c

[HecaiYuan]

Environment

• Arch: loongarch64
• CPU Model: Loongson‑3C6000/D
• Total CPU(s): 128 (2 sockets × 32 cores × 2‑threads‑per‑core), 4 NUMA nodes
• CPU flags: cpucfg lam ual fpu lsx lasx crc32 complex crypto ptw lvz lbt_x86 lbt_arm lbt_mips
• Build command:

make TARGET=LA464 DYNAMIC_ARCH=1 USE_THREAD=1 USE_OPENMP=1 DEBUG=1 \
     FC=gfortran CC=gcc \
     'COMMON_OPT=-O2 -g -fPIC -fopenmp -pthread' \
     'FCOMMON_OPT=-O2 -g -fPIC -fopenmp -pthread -frecursive'

make TARGET=LA464 DYNAMIC_ARCH=1 USE_THREAD=1 USE_OPENMP=1 DEBUG=1 \
FC=gfortran CC=gcc \
COMMON_OPT="-O2 -g -fPIC -fopenmp -pthread" \
FCOMMON_OPT="-O2 -g -fPIC -fopenmp -pthread -frecursive" \
lapack-test

1.
When running `lapack‑test`, double‑precision LAPACK test cases crash with segmentation fault inside LoongArch64 dgemm assembly kernel:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0  0x7ffff24ed3b3 in ???
#1  0x7ffff24ec0ff in ???
#2  0x7ffffcfc096f in ???
#3  0x1201d8b20 in ???
	at ../kernel/loongarch64/dgemm_kernel_16x6.S:423
Testing DOUBLE PRECISION Linear‑Equation‑routines‑LIN/xlintstd < dtest.in > dtest.out ---- TESTING ... FAILED

2.
Additionally I observed **infinite loop in cpotrf** under this build configuration.

[HecaiYuan]

File lapack/potrf/potrf_L_single.c:

#define GEMM_PQ  MAX(GEMM_P, GEMM_Q)
//leave some space for GEMM_ALIGN in sb2
#define REAL_GEMM_R (GEMM_R - 2*GEMM_PQ)

Concrete calculation with LA464 GEMM parameters

Target LA464 uses these parameters:

TABLE_NAME.dgemm_p = 112;
TABLE_NAME.dgemm_q = 289;
TABLE_NAME.dgemm_r = 353;

TABLE_NAME.cgemm_p = 128;
TABLE_NAME.cgemm_q = 256;
TABLE_NAME.cgemm_r = 512;

Calculation:

1. dpotrf (double): GEMM_PQ = MAX(112,289)=289 REAL_GEMM_R = 353 - 2*289 = -225
2. cpotrf (complex): GEMM_PQ = MAX(128,256)=256 REAL_GEMM_R = 512 - 2*256 = 0

These values feed into this for‑loop:
for(js = j + bk + min_j; js < n; js += REAL_GEMM_R){
....
}

- REAL_GEMM_R = -225: Negative loop step, the variable js decreases continuously, causing out-of-bounds memory access and SIGSEGV in dgemm assembly kernel.
- REAL_GEMM_R = 0: Zero loop step,js never changes, resulting in infinite loop in cpotrf.

[HecaiYuan]

Hi @martin-frbg,

Verification

After adjusting architecture GEMM parameters to satisfy GEMM_R > 2 * MAX(GEMM_P, GEMM_Q), both the SIGSEGV in dpotrf and infinite‑loop in cpotrf disappear completely.

I have two further questions:

1. What is the original design rationale for computing REAL_GEMM_R as GEMM_R - 2*GEMM_PQ?
2. Are there explicit size constraints for GEMM_R, GEMM_P, GEMM_Q for this potrf code path?

[HecaiYuan]

For reference, you can reproduce the SIGSEGV bug with the following C program

// dpotrf_N50.c
// gcc ./dpotrf_N50.c -o dpotrf_N50 ./libopenblas.a -I./lapack-netlib/LAPACKE/include/ -lm -lpthread -lgfortran -lgomp -g
#include <stdio.h>
#include <stdlib.h>
#include "lapacke.h"

int main() {
    int N = 50;
    double *A = calloc(N*N, sizeof(double));
    for (int j = 0; j < N; j++) {
        A[j*N + j] = 100.0;
        for (int i = j+1; i < N; i++)
            A[i*N + j] = 0.1;
    }
    int info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', N, A, N);
    printf("dpotrf info=%d\n", info);
    free(A);
    return 0;
}

[martin-frbg]

After adjusting architecture GEMM parameters to satisfy GEMM_R > 2 * MAX(GEMM_P, GEMM_Q), both the SIGSEGV in dpotrf and infinite‑loop in cpotrf disappear completely.
Thanks. Looks like this code path needs a compile-time check that the result is positive and non-zero for every target...

1. What is the original design rationale for computing `REAL_GEMM_R` as `GEMM_R - 2*GEMM_PQ`?

Honestly I doubt that anybody still knows - this is basically original GotoBLAS code, written about 20-25 years ago and abandoned by Kazushige Goto in 2010 when he took up commercial employment. libgoto was made available as occasional source releases with no documentation beyond the handful of files named GotoBLAS0* that we still keep. If he used any version control system at all back then, its history has been lost.

From what little is in our git history, the original actually had REAL_GEMM_R as just GEMM_R - GEMM_PQ , and Zhang Xianyi must have encountered an array overrun on some unknown platform in 2013 shortly after he revived GotoBLAS as OpenBLAS (one of several such attempts at the time, but the one that caught on). His solution was to introduce the factor of 2 "to leave room for alignment" (though the alignment formula in interface/lapack/potrf.c matches the one used in interface/gemm.c), which somehow appears to have worked for everybody until now. I took over maintainership from him when he wanted to refocus on RISCV hardware design in 2017, and he hasn't been very active since.

2. Are there explicit size constraints for `GEMM_R`, `GEMM_P`, `GEMM_Q` for this potrf code path?

None documented, perhaps none intentional. I'd think it is possible that the modification of the formula specifically for L_single only happened to mask a problem elsewhere in the code.

[martin-frbg]

That said, (a) this special table for Loongson cpu data used in the DYNAMIC_ARCH case (kernel/setparam-ref.c, from #4593) has curiously small numbers for GEMM_R compared to what's in param.h - and I finally realize I already "fixed" #5602 by increasing the R value in #5611 earlier this year, apparently without noticing the worse danger in the potrf_L_single branch . I wonder if removing dgemm_r from that special table altogether to use the larger default would have a noticeable impact on performance

[HecaiYuan]

That said, (a) this special table for Loongson cpu data used in the DYNAMIC_ARCH case (kernel/setparam-ref.c, from #4593) has curiously small numbers for GEMM_R compared to what's in param.h - and I finally realize I already "fixed" #5602 by increasing the R value in #5611 earlier this year, apparently without noticing the worse danger in the potrf_L_single branch . I wonder if removing dgemm_r from that special table altogether to use the larger default would have a noticeable impact on performance

I’ll run benchmarks on Loongson to find a balanced dgemm_r value.

[martin-frbg]

I do not think there is a theoretical justification for that factor of 2, so the alternative would be to remove it from the calculation and try to find the target configuration that fails under those conditions. At first glance, only the outdated DEC Alpha cpus stand out as having extremely small P and Q that might have required the doubling.