Broken BLAS.gemm with complex rectangular matrix

question
package
blas

#1

Hi,
I ran the following short code on two different versions of Julia and got results that are broken in different ways.

import Base.LinAlg.BLAS.gemm

A = convert(Array{Complex64,2}, [-6im -5im; -4 3; 2im im])
println(size(A))

println("plain *")
println(transpose(A)*A)
println(ctranspose(A)*A)

println("blas gemm")
println(gemm('t','n',A,A))
println(gemm('c','n',A,A))

output from Julia 0.4.5 (provided as a Linux Mint package)

(3,2)
plain *
Complex{Float32}[-24.0f0 + 0.0f0im -44.0f0 + 0.0f0im
                 -44.0f0 + 0.0f0im -17.0f0 + 0.0f0im]
Complex{Float32}[56.0f0 + 0.0f0im 20.0f0 + 0.0f0im
                 20.0f0 + 0.0f0im 35.0f0 + 0.0f0im]
blas gemm
Complex{Float32}[-24.0f0 + 0.0f0im -44.0f0 + 0.0f0im -5.835328f9 - 0.19085431f0im
                 -44.0f0 + 0.0f0im -17.0f0 + 0.0f0im 4.376496f9 - 0.095427155f0im]
Complex{Float32}[56.0f0 + 0.0f0im 20.0f0 + 0.0f0im -5.835328f9 + 0.19085431f0im
                 20.0f0 + 0.0f0im 35.0f0 + 0.0f0im 4.376496f9 + 0.095427155f0im]

output from nightly build 0.6.0-2a19c36283

(3, 2)
plain *
Complex{Float32}[-24.0+0.0im -44.0+0.0im; -44.0+0.0im -17.0+0.0im]
Complex{Float32}[56.0+0.0im 20.0+0.0im; 20.0+0.0im 35.0+0.0im]
blas gemm
ERROR: LoadError: DimensionMismatch("A has size (2,3), B has size (2,3), C has size (2, 3)")
Stacktrace:
 [1] gemm!(::Char, ::Char, ::Complex{Float32}, ::Array{Complex{Float32},2}, ::Array{Complex{Float32},2}, ::Complex{Float32}, ::Array{Complex{Float32},2}) at ./linalg/blas.jl:1028
 [2] gemm(::Char, ::Char, ::Complex{Float32}, ::Array{Complex{Float32},2}, ::Array{Complex{Float32},2}) at ./linalg/blas.jl:1042
 [3] gemm(::Char, ::Char, ::Array{Complex{Float32},2}, ::Array{Complex{Float32},2}) at ./linalg/blas.jl:1045
 [4] include_from_node1(::String) at ./loading.jl:539
 [5] include(::String) at ./sysimg.jl:14
 [6] process_options(::Base.JLOptions) at ./client.jl:305
 [7] _start() at ./client.jl:371
while loading /home/me/src/jl/test-blas.jl, in expression starting on line 14

gemm works ok with the square matrices that I’ve tried so far.
I think I reduced the problem to a point where all possible external causes have been ruled out.

Should I go ahead and open a GitHub issue?

Thanks.
David


#2

Use upper case flags ‘C’, ‘T’, ‘N’. gemm isn’t exported, so documentation is incomplete.


#3

It should possibly be doing more input validation if that doesn’t harm performance.


#4

Arguments are not allowed to alias in BLAS calls. That is a Fortran rule that allows for better compiler optimizations. See e.g. https://software.intel.com/en-us/blogs/2009/07/10/doctor-fortran-in-ive-come-here-for-an-argument-side-2. The functions in the BLAS module are considered expert functions and should generally not be used directly. What is the reason you call them?


#5

Thanks for your replies.
Using the upper case ‘C’, ‘T’ worked.

I was trying to do quantum mechanics with BLAS/LAPACK in the space of complex, highly rectangular m×n-matrices (n≫m~10).


#6

Notice that * calls BLAS.