The netlib documentation of sgemm states that the array strides LDA and LDB must be >= 1, and large enough so that columns don’t overlap. Indeed, the implementation in Apple’s Accelerate/veclib framework checks these conditions and exists if they are violated.
I really don’t understand why this limitation exists. Can’t BLAS simply believe me that I really want a stride of zero, or a negative stride? As far as I understand, Fortran integers are signed by default, so parameter types don’t seem to be the reaons (disclaimer: I don’t know Fortan very well).
Indeed, there exist very reasonable use cases of non-positive array strides:
- zero stride: in a multi-dimensional array class, a stride of zero enables numpy-style broadcasting.
- negative stride: negating a stride allows to view an array in reverse order along any axis without copying. This can be useful e.g. when flipping convolution kernels, and convolutions can be implemented efficiently using gemm. Alternatively, the vertical axis of an image can be flipped, which is handy since there exist different conventions: axis pointing upwards in postscript/pdf, and downwards in the png format (and many others).
I am interested in two aspects:
- I’d like to understand why the restriction exists. Is it really just because the designers of BLAS did not think about such use cases? Am I the victim of someone trying to catch a bug that’s indeed a feature? Or does the restriction result in better performance? I have a hard time imagining the latter.
- Is there a way to work around the restriction without sacrificing (too much) performance? For now I need something that works on a Mac in C++, but long term it should still be based on BLAS, so it works across platforms.
You need to sign in to view this answers
