OpenBLAS sgemm, особые требования?иногда я получаю нан - PullRequest
Новая
       58

OpenBLAS sgemm, особые требования?иногда я получаю нан

0 голосов
/ 13 февраля 2019

У меня где-то есть это объявление в моем заголовке: 

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 в моем результате; и, похоже, этого не происходит в нашем производственном коде).

1 Ответ

0 голосов
/ 13 февраля 2019

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

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

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

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

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

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

...