c++ - 为什么在 blas gemm 函数族中不允许非正步幅？

Question

The netlib documentation of sgemm指出数组步长 LDA 和 LDB 必须是 >= 1，并且足够大，这样列就不会重叠。事实上，Apple 的 Accelerate/veclib 框架中的实现会检查这些条件并在违反时存在。

我真的不明白为什么会有这个限制。 BLAS 不能简单地相信我真的想要零步幅或负步幅吗？据我了解，默认情况下 Fortran 整数是有符号的，因此参数类型似乎不是原因(免责声明:我不太了解 Fortan)。

确实，存在非正数组步长的非常合理的用例:

• 零步幅:在多维数组类中，零步幅支持 numpy 风格的广播。
• 负步幅:取反步幅允许沿任何轴以相反顺序查看数组而无需复制。这可能很有用，例如翻转卷积核时，and convolutions can be implemented efficiently using gemm .或者，可以翻转图像的垂直轴，这很方便，因为存在不同的约定:轴在 postscript/pdf 中指向上方，而在 png 格式(以及许多其他格式)中指向下方。

我对两个方面感兴趣:

1. 我想了解为什么存在限制。真的只是因为BLAS的设计者没有考虑过这样的用例吗？我是不是有人试图捕捉一个确实是功能的错误的受害者？还是限制会带来更好的性能？我很难想象后者。
2. 有没有办法在不牺牲(太多)性能的情况下解决这个限制？现在我需要一些可以在 Mac 上使用 C++ 的东西，但从长远来看它应该仍然基于 BLAS，所以它可以跨平台工作。

Answer 1

我最近发现自己有效地做到了这一点:

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 的某些实现会让您默认执行此操作，而其他实现将尝试防止无限循环，作为副作用，无限循环将假定为正步幅值。

---

Intel BLAS 似乎暗示它支持负步幅，因为它提到了 abs(incx)，至少对于 AXPY:https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/blas-and-sparse-blas-routines/blas-routines/blas-level-1-routines-and-functions/cblas-axpy.html

我没有尝试，是否适用于其他功能还有待观察。