-
Notifications
You must be signed in to change notification settings - Fork 68
/
reduction_untangling.jl
52 lines (50 loc) · 1.22 KB
/
reduction_untangling.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# issue 288
function not_a_reduction!(A, B)
@turbo for j ∈ 1:size(A.re, 2)
jre = B.re[j]
jim = B.im[j]
for i ∈ 1:size(A.re, 1)
ire = B.re[i]
iim = B.im[i]
cisim = iim * jre - ire * jim
cisre = ire * jre + iim * jim
ρre_i = A.re[i, j]
ρim_i = A.im[i, j]
re_out = ρre_i * cisre - ρim_i * cisim
im_out = ρre_i * cisim + ρim_i * cisre
A.re[i, j] = re_out
A.im[i, j] = im_out
end
end
return nothing
end
function not_a_reduction_noturbo!(A, B)
@turbo for j ∈ 1:size(A.re, 2)
jre = B.re[j]
jim = B.im[j]
for i ∈ 1:size(A.re, 1)
ire = B.re[i]
iim = B.im[i]
cisim = iim * jre - ire * jim
cisre = ire * jre + iim * jim
ρre_i = A.re[i, j]
ρim_i = A.im[i, j]
re_out = ρre_i * cisre - ρim_i * cisim
im_out = ρre_i * cisim + ρim_i * cisre
A.re[i, j] = re_out
A.im[i, j] = im_out
end
end
return nothing
end
@testset "Untangle reductions" begin
@show @__LINE__
N = 11
A1 = (re = rand(N, N), im = rand(N, N))
A2 = deepcopy(A1)
B = (re = rand(N), im = rand(N))
not_a_reduction!(A1, B)
not_a_reduction_noturbo!(A2, B)
@test A1.re ≈ A2.re
@test A1.im ≈ A2.im
end