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). И тогда я прямо звоню sgemm_ в моем коде C++. (Код должен в принципе работать с любой библиотекой BLAS.)

Я не уверен, что это плохая идея или нет. Или о чем я должен позаботиться. Например, мне нужно специальное выравнивание? Или мне нужно позаботиться в многопоточной среде?

(Например, я немного искал код OpenBLAS (SGEMM специально для ядра), и, похоже, он предполагает особые требования к выравниванию (но, возможно, я ошибся).)

В основном, кажется, работает нормально. За исключением того, что в некоторых случаях (недетерминированных, может быть, в 10% случаев, для некоторых сложных тестовых случаев, я получаю nan в моем результате; и это, кажется, не происходит в нашем производственном коде).

Answer 1

У вас не должно возникнуть проблем с вызовом BLAS из C++. OpenBLAS является широко используемой библиотекой, и не ожидается, что у нее будет такая базовая проблема.

BLAS - это библиотека Fortran, поэтому вы должны иметь в виду, что она использует формат хранения Fortran для двумерных матриц. Это означает, что матрицы - это отдельные блоки памяти, в которых хранятся двумерные матрицы в главном порядке столбцов.

Вы не можете использовать двумерные динамически распределенные массивы (т.е. double **a;) так как выделенная память будет фрагментирована. Также, если вы используете двумерные статические массивы (т.е. double a[5][4]), вы должны иметь в виду, что в C/C++ порядок хранения является основным. В этом случае вы все еще можете использовать BLAS, но вы должны учитывать, что матрица транспонирована.

Я бы предложил использовать векторы с одним указателем (double *a;) и получить доступ к элементам матрицы вручную (a[i+j*m]).

OpenBLAS имеет многопоточную поддержку. При компиляции вы можете определить, будет ли он использовать потоки или OpenMP или ничего.

Что касается получаемой ошибки, я бы посоветовал проверить вашу память, поскольку такое поведение обычно основано на ошибках в памяти. В любом случае я не ожидал, что ошибка будет sgemm реализация, но так, как это вызывается.