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 格式(以及许多其他格式)中指向下方。
我对两个方面感兴趣:
- 我想了解为什么存在限制。真的只是因为BLAS的设计者没有考虑过这样的用例吗?我是不是有人试图捕捉一个确实是功能的错误的受害者?还是限制会带来更好的性能?我很难想象后者。
- 有没有办法在不牺牲(太多)性能的情况下解决这个限制?现在我需要一些可以在 Mac 上使用 C++ 的东西,但从长远来看它应该仍然基于 BLAS,所以它可以跨平台工作。
最佳答案
我最近发现自己有效地做到了这一点:
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
我没有尝试,是否适用于其他功能还有待观察。
关于c++ - 为什么在 blas gemm 函数族中不允许非正步幅?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/44875236/