c - OpenBLAS sgemm，特殊要求？有时我得到南

Question

我的标题中某处有此声明:

typedef float real;
typedef int integer;
extern "C" {
extern int sgemm_(char *transa, char *transb,
  integer *m, integer *n, integer *k,
  const real *alpha,
  const real *a, integer *lda,
  const real *b, integer *ldb,
  const real *beta,
  real *c, integer *ldc);
}

然后我链接到 OpenBLAS 库(或者也可以选择其他 BLAS 库，例如 MKL)。然后我在我的 C++ 代码中直接调用 sgemm_。 (该代码原则上应该适用于任何 BLAS 库。)

我不确定这是否是个坏主意。或者我应该注意什么。例如。我需要特别对齐吗？还是我需要在多线程环境中小心？

(例如，我仔细查看了 OpenBLAS 代码(特别是 SGEMM 内核)，看起来它假定了特殊的对齐要求(但也许我弄错了)。)

它似乎大部分都运行良好。除了在某些情况下(不确定的，可能是 10% 的情况，对于一些复杂的测试用例，我在结果中得到 nan ；而且它似乎不会发生在我们的生产代码中)。

Answer 1

从 C++ 调用 BLAS 应该没有任何问题。 OpenBLAS 是一个广泛使用的库，预计不会有这样的基本问题。

BLAS 基本上是一个 Fortran 库，因此您必须记住，它对二维矩阵使用 Fortran 存储格式。这意味着矩阵是按列主要顺序存储二维矩阵的单个内存块。

您不能使用二维动态分配的数组(即 double **a;)，因为分配的内存将是碎片化的。此外，如果您使用二维静态数组(即 double a[5][4])，您应该记住，在 C/C++ 中，存储顺序是行优先的。在这种情况下，您仍然可以使用 BLAS，但您必须考虑矩阵已转置。

我建议使用单指针 vector (double *a;) 并手动访问矩阵元素 (a[i+j*m])。

OpenBLAS 支持多线程。在编译时，您可以定义它是否使用线程或 OpenMP 或什么都不使用。

至于你得到的错误，我建议检查你的内存，因为这种行为通常是基于内存错误。在任何情况下，我都不希望错误出现在 sgemm 实现上，而是出现在调用它的方式上。