Почему неположительные шаги запрещены в семействе функций Blas Gemm?

Question

Почему неположительные шаги запрещены в семействе функций Blas Gemm?

Netlib документация sgemm утверждает, что массив шагает LDA а также LDB должно быть >= 1и достаточно большой, чтобы столбцы не перекрывались. Действительно, реализация в Apple Accelerate / veclib Framework проверяет эти условия и существует, если они нарушаются.

Я действительно не понимаю, почему существует это ограничение. Разве BLAS не может просто поверить мне, что я действительно хочу шаг в ноль или шаг в отрицательном направлении? Насколько я понимаю, целые числа Фортрана подписаны по умолчанию, поэтому типы параметров, кажется, не являются причиной (отказ от ответственности: я не очень хорошо знаю Фортан).

Действительно, существуют очень разумные случаи использования неположительных шагов массива:

нулевой шаг: в классе многомерных массивов нулевой шаг позволяет осуществлять вещание в стиле numpy.
отрицательный шаг: отрицание шага позволяет просматривать массив в обратном порядке по любой оси без копирования. Это может быть полезно, например, при переворачивании ядер свертки, и свертки могут быть эффективно реализованы с использованием гемм. В качестве альтернативы можно перевернуть вертикальную ось изображения, что удобно, поскольку существуют разные соглашения: ось направлена вверх в postscript / pdf и вниз в формате png (и многие другие).

Меня интересуют два аспекта:

Я хотел бы понять, почему существует ограничение. Неужели дизайнеры BLAS не задумывались о таких случаях использования? Я жертва кого-то, пытающегося поймать ошибку, которая действительно особенность? Или ограничение приводит к лучшей производительности? Мне трудно представить себе последнее.
Есть ли способ обойти ограничение, не жертвуя (слишком много) производительностью? На данный момент мне нужно что-то, что работает на Mac в C ++, но в долгосрочной перспективе оно должно по-прежнему основываться на BLAS, поэтому оно работает на разных платформах.

1

arrays blas c++stride

Решение

Другие решения

Других решений пока нет …

Источник

Accepted Answer

Я недавно обнаружил, что эффективно делаю это:

double s[4] = {10., 2., 3.1, 4.1};
dscal_(4, 3., s, -1);
assert( s[1] == 2.*3. );

dscal_ является самой простой функцией BLAS, умножающей массив на скаляр, ее сигнатура:

void sscal(int, double, double*, int); // doesn't seem to return anything useful

В моем конкретном дистрибутиве BLAS (поставляется с Fedora 28) это выдало тихую ошибку, так как функция, похоже, ничего не делает.
Более того, dscal_ похоже, даже не возвращает код ошибки, поэтому я не могу отловить эту ошибку без функции обтекания (у моих массивов есть положительные или отрицательные результаты во время выполнения, но я не могу контролировать значение во всех случаях).

Все эти случаи потерпели неудачу:

double s[4] = {10., 2., 3.1, 4.1};
dscal_(4, 3., s, -1); // negative incx does nothing
dscal_(4, 3., &s[4], -1); // negative incx does nothing
dscal_(-4, 3., &s[4], 1); // negative n does nothing
dscal_(-4, 3., &s[4], -1); // negative n or incx does nothing
assert( s[1] == 2. );

Это говорит мне о том, что, хотя, вероятно, нигде не задокументированы шаг (incx) должен быть положительным, (а также размер).
К счастью, для многих функций BLAS вызов может быть преобразован в положительные результаты.
Мне нужна была функция-обертка, чтобы отладить это любым способом, поэтому написал следующую функцию-обертку:

void signed_dscal(int n, double alpha, double* x, int incx){
int prod = incx*n;
if(prod > 0) dscal(abs(n), alpha, x, abs(incx));
else         dscal(abs(n), alpha, x + prod, abs(incx));
}

Таким образом, я могу позвонить signed_dscal с положительными или отрицательными шагами и размерами.

int main(){
{
double d[4] = {10., 2., 3.1, 4.1};
signed_dscal(4, 3., d, 1);
assert( d[1] == 6. );
}
{
double d[4] = {10., 2., 3.1, 4.1};
signed_dscal(4, 3., &d[4], -1);
assert( d[1] == 6. );
}
{
double d[4] = {10., 2., 3.1, 4.1};
signed_dscal(-4, 3., &d[4], 1);
assert( d[1] == 6. );
}
{
double d[4] = {10., 2., 3.1, 4.1};
signed_dscal(-4, 3., d, -1);
assert( d[1] == 6. );
}
return 0;
}

(Обратите внимание, что incx=0 до сих пор дело нельзя исправить.)

Я до сих пор не понимаю, в чем логика этого. Возможно, некоторые реализации BLAS позволят вам сделать это по умолчанию, а другие попытаются защитить от бесконечных циклов, которые в качестве побочного эффекта примут положительные значения шага.

0