OpenBLAS #83: simple GEMM test with magic(5)

zc_gemm_magic5.c

// https://github.com/xianyi/OpenBLAS/issues/83
//
// gcc -O0 -g3 -ggdb3 -std=gnu99 -Iinstall\include zc_gemm_magic5.c -o zc_gemm_magic5 -l:install\lib\libopenblas.lib

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <openblas/cblas.h>

#if 1
  typedef float float_type;
  #define cblas_gemm cblas_sgemm
#else
  typedef double float_type;
  #define cblas_gemm cblas_dgemm
#endif

const float_type magic5[] = {
  17.0, 24.0,  1.0,  8.0, 15.0,
  23.0,  5.0,  7.0, 14.0, 16.0,
   4.0,  6.0, 13.0, 20.0, 22.0,
  10.0, 12.0, 19.0, 21.0,  3.0,
  11.0, 18.0, 25.0,  2.0,  9.0
};
const int N = 5;

// expected = magic5 * magic5
const float_type expected[] = {
  1090.0,  900.0,  725.0,  690.0,  820.0,
   850.0, 1075.0,  815.0,  720.0,  765.0,
   700.0,  840.0, 1145.0,  840.0,  700.0,
   765.0,  720.0,  815.0, 1075.0,  850.0,
   820.0,  690.0,  725.0,  900.0, 1090.0
};

int main() {
  const int magic5size = N * N * sizeof(*magic5);
  float_type* A = malloc(magic5size);
  float_type* B = malloc(magic5size);
  float_type* C = malloc(magic5size);
  memcpy(A, magic5, magic5size);
  memcpy(B, magic5, magic5size);

  cblas_gemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N, N, N, 1.0, A, N, B, N, 0.0, C, N);

  printf("sizeof(float) = %ld\n\n", (long) sizeof(float_type));

  printf("C = [\n");
  for (int i = 0; i < N; ++i) {
    printf("  ");
    for (int j = 0; j < N; ++j) {
      printf("%f, ", C[i*N + j]);
    }
    printf("\n");
  }
  printf("]\n\n");

  printf("C ./ expected = [\n");
  for (int i = 0; i < N; ++i) {
    printf("  ");
    for (int j = 0; j < N; ++j) {
      printf("%f, ", C[i*N + j] / expected[i*N + j]);
    }
    printf("\n");
  }
  printf("]\n");

  free(A);
  free(B);
  free(C);
  return 0;
}