Support StridedArrays with strides not supported by BLAS routines

[samuela]

I'm getting "Intel MKL ERROR: Parameter 10 was incorrect on entry to DGEMM ." printed to stdout with no other indication of what went wrong.

[dmbates]

More context will be needed. DGEMM performs multiplication of dense matrices. Does that help you to narrow down where the error is occurring? Also, it will help to have the output of

versioninfo()

[jiahao]

Please post a minimal example of code that triggered this problem and the result of versioninfo(). Otherwise it is impossible to diagnose this issue other than to point you to the signature of dgemm, where parameter 10 is LDB, the leading dimension of B in a A*B matrix product.

[samuela]

Yeah, I'm trying to track it down right now. I'll post as I figure it out.

[samuela]

julia> versioninfo()
Julia Version 0.3.4-pre+84
Commit b296937 (2014-12-22 13:18 UTC)
Platform Info:
  System: Linux (x86_64-redhat-linux)
  CPU: Intel(R) Xeon(R) CPU E5-2670 v2 @ 2.50GHz
  WORD_SIZE: 64
  BLAS: libmkl_rt
  LAPACK: libmkl_rt
  LIBM: libimf
  LLVM: libLLVM-3.3

[samuela]

Alright, this reproduces exactly the error:

using ArrayViews

function moving_window(arr, n)
    shp = (size(arr, 1) * n, size(arr, 2) - n + 1)
    strided_view(arr, shp, ArrayViews.contrank(arr), strides(arr))
end

x = rand(5, 5)

rand(10, 4) * moving_window(x, 2)

[samuela]

Oops, that should be rand(4, 10) instead of rand(10, 4).

[jiahao]

Thanks. I can reproduce the bug. Parameter 10 is calculated in linalg/blas.jl:gemm! as max(1,stride(B,2)) = 5, which is strange since B is a 10x4 matrix. The stride is probably being calculated incorrectly.

[samuela]

Well, I'm specifying a custom stride in moving_window which likely causes the confusion.

[samuela]

Can DGEMM even handle a matrix laid out in this way? I'm not familiar with BLAS.

[jiahao]

I'm not sure exactly what your code is supposed to do, but if you set

julia> x=reshape(1.:25.,5,5)
5x5 Array{Float64,2}:
 1.0   6.0  11.0  16.0  21.0
 2.0   7.0  12.0  17.0  22.0
 3.0   8.0  13.0  18.0  23.0
 4.0   9.0  14.0  19.0  24.0
 5.0  10.0  15.0  20.0  25.0

then the B in gemm! is

[1.0 6.0 11.0 16.0
 2.0 7.0 12.0 17.0
 3.0 8.0 13.0 18.0
 4.0 9.0 14.0 19.0
 5.0 10.0 15.0 20.0
 6.0 11.0 16.0 21.0
 7.0 12.0 17.0 22.0
 8.0 13.0 18.0 23.0
 9.0 14.0 19.0 24.0
 10.0 15.0 20.0 25.0]

Is that the expected matrix?

[samuela]

Yup, exactly

[jiahao]

If the desired stride is indeed 5, then BLAS cannot handle this computation. The documentation linked to above states that the parameter LDB to DGEMM must be at least the common dimension (10 in this case).

I'm not really sure what to do here. Maybe we can insert a runtime check in LinAlg.BLAS.gemm! if the strides are what BLAS expects and fall back to generic_matmul! otherwise?