Multiply two complex matrices fails (some time)

[wenxgwen]

Here is the julia session and the code:

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" to list help topics
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.2.0-rc2-prerelease+4331
 _/ |\__'_|_|_|\__'_|  |  Commit 902b446 2013-10-29 23:03:55 UTC
|__/                   |  i686-linux-gnu

julia> reload("test")
ERROR: InexactError()
 in setindex! at array.jl:413
 in matmul2x2 at linalg/matmul.jl:594
 in generic_matmatmul at linalg/matmul.jl:412
 in generic_matmatmul at linalg/matmul.jl:398
 in * at linalg/matmul.jl:391
 in include at boot.jl:238
 in include_from_node1 at loading.jl:114
 in reload_path at loading.jl:140
 in reload at loading.jl:73
at /home/wen/Dropbox/prog/julia/STK/test.jl:27

julia> S
2x2 Array{Complex{T<:Real},2}:
 1.0+1.46958e-15im  1.0+1.95943e-15im
 1.0+1.95943e-15im  1.0+9.79965e-15im

julia> T
2x2 Array{Complex{T<:Real},2}:
 1.0-7.34788e-16im               0+0im
             0+0im  -1.0+4.89983e-15im

julia> S*T
ERROR: InexactError()
 in setindex! at array.jl:413
 in matmul2x2 at linalg/matmul.jl:594
 in generic_matmatmul at linalg/matmul.jl:412
 in generic_matmatmul at linalg/matmul.jl:398
 in * at linalg/matmul.jl:391

julia> 

The code test.jl :

dim=2

K=[1 0; 0 2]
V=[1 1; 1 0]

# The following works:  S*T does not fail
# S=zeros(Complex{Float64},dim,dim)
# T=zeros(Complex{Float64},dim,dim)

# The following does not work: S*T fails
S=zeros(Complex,dim,dim)
T=zeros(Complex,dim,dim)

i=0
for i1 = 1:dim
  i+=1
  l1=[1; i1]
  T[i,i] =exp(1im*(pi*(V*l1)'*K*(V*l1))[1,1])
  j=0
for i2 = 1:dim
  j+=1
  l2=[1; i2]
  S[i,j]=exp(-2im*pi*((V*l2)'*K*(V*l1))[1,1])
end
end

S*T

Xiao-Gang

[pao: formatting; GitHub uses triple-backticks for code blocks]

[ViralBShah]

I get the same thing.

@JeffBezanson @StefanKarpinski We should take a look and see if this is relevant for 0.2

[StefanKarpinski]

If this has never come up before, I have a hard time believing that it's high enough priority to go in 0.2.

[ViralBShah]

If @JeffBezanson can take a quick look, that would be nice. I agree that it is not unsettling enough to hold the release.

[timholy]

For one thing, even if it does get fixed you should never do the version that currently fails.

[timholy]

@wenxgwen, if my last comment seems at all mysterious to you, see http://docs.julialang.org/en/latest/manual/faq/#how-do-abstract-or-ambiguous-fields-in-types-interact-with-the-compiler

[Keno]

I'm pretty sure the difference is in the type of the output array. Generic matmul does C = Array(typeof(x*y+x*y), mA, nB) while matmul2x2 does Array(promote_type(T,S), 2, 2). The former causes problems for obvious reasons (x=one(Complex) may not be the right type)

[StefanKarpinski]

For one thing, even if it does get fixed you should never do the version that currently fails.

"never" seems a bit strong. I'm sure there's situations where you don't want that extra type parameter.

[jiahao]

The issues with Array{Complex{T}, 2} are much more extensive than what is reported in this issue. #4807

[jiahao]

I'm quite disturbed by this bug. Creating Complex instead of Complex64/128 variables is exactly the kind of thing that a new user would try, have break, and give up in frustration.

[wenxgwen]

@jiahao
Indeed, new user like me will use Complex instead of Complex64/128. Most end user have no clue what is 64 or 128. If Complex64 may be used why not complex13, or complex8, or Complex1000.

@timholy
Thank you very much for the link. In practice most people will use Complex, and I do hope Julia can generate fast code for Complex.

[jiahao]

In practice most people will use Complex, and I do hope Julia can generate fast code for Complex.

@wenxgwen It is not possible to generate efficient code with just Complex for the reasons explained in the link provided by @timholy. Most users will need to choose a precision in Complex64/128, just as in Fortran you need to choose complex*8 or complex*16. It is not hard to educate users about this choice. Complex{T} is more general; e.g. you can have Complex{BigFloat} to do arbitrary precision floating point.

The real issue here is that arrays of Complex{T} are just broken. When fixed, it still won’t be as fast as Complex64/128, and that’s a feature of Julia that probably won’t change.

[ivarne]

Maybe a rename is in order? Like: Complex-> ComplexNumber, and alias ComplexNumber{128} as Complex?

[timholy]

@wenxgwen, I'll echo what @jiahao said, what you're asking for is simply not possible. It's not a language thing, it's a CPU architecture and assembly instruction thing. You have to compile it down to something that uses the hardware to do the computations (there's a reason CPUs contain FPUs...), and you can only do that when the types are known concretely. @StefanKarpinski, perhaps "never" was too strong, but certainly this should be an extreme corner case.

I'm mostly in favor of @ivarne's suggestion, it would definitely circumvent the new user problem. One downside is that it introduces a slight asymmetry between the naming schemes for real and complex (there's no "alias" for Float64 currently).

[ivarne]

To be consistent with Float64, we might disregard Complex altogether. That would also make sure we break all existing Complex code, and not leave some old code working, but with different results.

We already have ComplexPair as an alias for Complex, and that is a better name. Maybe the alias and type should be switched so that Complex can be deprecated? (if we could deprecate an type alias? feature request?)

[wenxgwen]

I'm also in favor of @ivarne's suggestion. For example we can use ComplexType and RealType to mean generic abstract type and use Complex and Real to mean Complex(128) and Float{64}. This might be more user friendly.

[andreasnoack]

What about making zero(AbstractType) and one(AbstractType) return an error?

[ivarne]

Error for abstract types are a very good idea. Same for AbstractType[]. There are some cases where a default is obvious, but in the general case it is not, and I can't see any reason why we should not encourage use of abstract typed variables.

I made an issue for this in #4808