arm64_bug.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "cblas.h"


int main() {
    int found_error;
    int k;
    int i;
    int n = 1000;
    int m = 10;
    while (m <= 1000) {
        k = 2;
        while (k <= 512) {
            found_error = 0;
            printf("dgemm m=%d, n=%d, k=%d\n", m, n, k);
            double *a = (double*) malloc(sizeof(double) * m * k);
            double *b = (double*) malloc(sizeof(double) * n * k);
            double *c = (double*) malloc(sizeof(double) * m * n);
            for (i=0; i<m*k; i++) {
                a[i] = 1.0;
            }
            for (i=0; i<n*k; i++) {
                b[i] = 1.0;
            }
            for (i=0; i<m*n; i++) {
                c[i] = 0.0;
            }
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1., a, k, b, k, 0.0, c, n);
            for (i=0; i<m*n; i++) {
                if (isnan(c[i]) | c[i] < k - 1e-12 | c[i] > k + 1e-12) {
                    printf("c[%d]=%f\n", i, c[i]);
                    found_error = 1;
                }
            }
            if (!found_error) {
                printf("dgemm ok\n");
            }
            free(a);
            free(b);
            free(c);
            k *= 2;
        }
        m *= 2;
    }
    return 0;
}