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 lib（或者可选地还有其他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，但你必须考虑矩阵是转置的。

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

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

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