Julia: Lost elements when .+= is used on sparse matrices

[falbarelli]

Hello,
I use Julia to work in Quantum Information and I think I found a bug involving sparse matrices and the .+= operator.
This is the first time I report a bug, so I hope I'm following the right etiquette.

The code below is a minimal example of the problem.
I use matrices defined via the QuantumOptics.jl package because these are the matrices that made me find the problem. However, these matrices are not part of the problem, unluckily I was not able to reproduce the bug on more simple and random matrices.

using QuantumOptics
function Jxj(j::Int, n::Int)
    b = SpinBasis(1//2)
    return full(tensor(vcat([identityoperator(b) for i = 1:j-1], [1/2*sigmax(b)], [identityoperator(b) for i = j + 1:n])...)).data
end
function ghz_state(n::Int)
    b = SpinBasis(1//2)
    return (tensor([spinup(b) for i in 1:n]...).data + tensor([spindown(b) for i in 1:n]...).data)/sqrt(2)
end

Nj=3
Lj = [Jxj(j,Nj) for j = 1:Nj]
Lj_sparse = [sparse(Jxj(j,Nj)) for j = 1:Nj]

rho = ghz_state(Nj) * ghz_state(Nj)'
rho_sparse=sparse(rho)

sparse_result1=spzeros(Complex{Float64},2^Nj,2^Nj)
sparse_result1 .+= Lj_sparse[1]*rho_sparse*Lj[1]
sparse_result1 .+= Lj_sparse[2]*rho_sparse*Lj[2]
#same behavior if I use: 
# sparse_result1 .= Lj_sparse[1]*rho_sparse*Lj[1] .+ sparse_result1
# sparse_result1 .= Lj_sparse[2]*rho_sparse*Lj[2] .+ sparse_result1

sparse_result2=spzeros(Complex{Float64},2^Nj,2^Nj)
sparse_result2.=Lj_sparse[1]*rho_sparse*Lj_sparse[1].+Lj_sparse[2]*rho_sparse*Lj_sparse[2]

sparse_result2-sparse_result1

The matrix sparse_result1 has two missing elements and thus is not equal to sparse_result2. The problem gets worse if the operator .+= is used again (i.e. in a for loop) and I keep loosing pieces of the matrix.
Moreover, with Nj=2 (i.e. 4 by 4 matrices, instead of 8 by 8) the bug is not present.
As written in the code I get the same bug also if I don't use .+= but I explicitly write the sum.

Theversioninfo() of my laptop:

Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-7360U CPU @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Prescott)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

but the same bug occurs with JuliaPro on a Linux workstation:

Julia Version 0.6.2
Commit d386e40* (2017-12-13 18:08 UTC)
Platform Info:
  OS: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2609 v3 @ 1.90GHz
  WORD_SIZE: 64
  BLAS: libmkl_rt
  LAPACK: libmkl_rt
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

[yuyichao]

Dup of https://github.com/JuliaLang/julia/issues/21693

[nalimilan]

Are you sure that's the same problem? I don't see where an input aliases the output here.

[yuyichao]

a .+= b is a .= a .+ b

[nalimilan]

Indeed. I couldn't believe we don't support that very standard operation!

[falbarelli]

Thanks everyone for the feedback.
For the moment i circumvented the problem by creating a temporary variable... instead of a .+= b I need to do c.= a .+ b; a .= c.
However, I believe that it is quite dangerous that this behaviour is not at least signalled by some warning/error message

[StefanKarpinski]

Agreed, that's why we have an open issue about it 👍