Bug in BLAS (and LAPACK) when using strided matrices

Previous Topic Next Topic
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view

Bug in BLAS (and LAPACK) when using strided matrices

Paulo Jabardo
I was explaining to a student how matrices are passed to BLAS and how they should be organized in memory when I ran the simple test below. I expected to get an error but what I got was a wrong answer.

Matrix arguments in BLAS (and LAPACK) accept StridedMatrices. But this is wrong as is shown in the following example:

using Base.LinAlg.BLAS

A = reshape(linspace(1,100,100), 10, 10)
y = [linspace(1,3,3);]

A1 = A[1:3:7, 2:3:8]
A2 = sub(A, 1:3:7, 2:3:8)
A3 = copy(A2)

z1 = gemv('N', 1.0, A1, y)
z2 = gemv('N', 1.0, A2, y)
z3 = gemv('N', 1.0, A3, y)

hcat(z1, z2, z3)

which returns

3x3 Array{Float64,2}:
 306.0  306.0  306.0
 324.0  312.0  324.0
 342.0  318.0  342.0

Notice the second column.

I think the solution is to specify a new matrix type, perhaps BlasMatrix that allows only UnitRange indexing in subarrays. This problem occurs in both version 0.4.3 and 0.5

I will file an issue right now.