Fail "sometimes" in unitary transformation

[pochoi]

I was using Expokit to compute some unitary transformations. However, I found out that the output is wrong "sometime".

using Expokit

badk = zeros(Int64, 0)
badx = Array{Array{Float64,1}}(0)
for k in 1:1000
    x = expmv(-0.1, sparse(9:-1:1, 1:9, -4:4), ones(9))
    if norm(x) > 4
        push!(badk, k)
        push!(badx, x)
    end
end

length(badk)

sparse(9:-1:1, 1:9, -4:4) is just a 9x9 skew-symmetric matrix, then its exponential is a unitary/orthonormal matrix. Therefore norm(x) should be approximately norm(ones(9)) = 3.

Depending on you luck, you will get non-empty badk above. It happens on my Macbook as well as my office linux computing server.

[acroy]

I can sort of confirm your observation using Julia v0.5. But it happens very rarely - I can increase the number of trials to 1.000.000 and the first run never failed so far. Only when I re-run the loop I get 1-2 elements in badk. My guess is that we need to set hm to zero inside expmv like we discussed for phimv in #18.

Try to replace x = expmv(-0.1, S, ones(9)) by x = phimv(-0.1, S, zeros(9), ones(9)) where S=sparse(9:-1:1, 1:9, -4:4) as a quick fix.

[acroy]

@pochoi : Please reopen if the problem persists. The fix is merged into master.

[pochoi]

Thanks so much for the quick reply and fix. The fix works nicely in both my systems! 🎉 (Actually both systems I am using are 0.6)

Julia Version 0.6.0
Commit 903644385b (2017-06-19 13:05 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i5-4258U CPU @ 2.40GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

Julia Version 0.6.0
Commit 9036443 (2017-06-19 13:05 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2699 v3 @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)