1

sgemm 的 netlib 文档指出数组跨步LDA并且LDB必须是>= 1，并且足够大，以便列不会重叠。事实上，Apple 的 Accelerate/veclib 框架中的实现会检查这些条件，如果它们被违反，它们就会存在。

我真的不明白为什么存在这种限制。BLAS 不能简单地相信我真的想要零步幅或负步幅吗？据我了解，Fortran 整数默认是有符号的，所以参数类型似乎不是原因（免责声明：我不太了解 Fortan）。

实际上，存在非常合理的非正数组步幅用例：

  • 零步幅：在多维数组类中，零步幅启用 numpy 样式的广播。
  • 负步幅：否定步幅允许沿任何轴以相反顺序查看数组，而无需复制。这可能很有用，例如在翻转卷积核时，可以使用 gemm 有效地实现卷积。或者，可以翻转图像的垂直轴，这很方便，因为存在不同的约定：轴在 postscript/pdf 中向上，在 png 格式（以及许多其他格式）中向下。

我对两个方面感兴趣：

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

1 回答 1

0

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

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

于 2018-12-17T07:36:59.730 回答