Broken BLAS.gemm with complex rectangular matrix

Specific Domains Numerics
, ,
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. Intel Developer Zone. 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.

1 Like