PyOpenCL - IntelとNVidiaで異なる結果

Question

私は、ガウスパルス伝播に基づくシミュレーションを実行しようとしています。私はi5 4590 & GTX 970（最新のドライバ）と私の2015年代初頭のMacbook Airを使って、私のWindowsデスクトップとの間でクロス開発しています。

私のメインコードを実行すると、私のデスクトップ上ではまともな結果が得られませんでした（計算が分岐しました）が、私のMacでは結果は大丈夫でした。

さらに問題を調べるために、シンプルガウス伝播を実行しようとしました。 MacBookの結果は多かれ少なかれOKですが、デスクトップでは完全に混乱します。

私は両方のマシンで同じコードを実行していますが、どちらもPython（2.7.10）とそれぞれのモジュールの同じディストリビューションを持っています。ここで

は私のコード

import scipy as sp 
import pyopencl as cl 
import matplotlib.pyplot as plot 

ctx = cl.create_some_context() 
queue = cl.CommandQueue(ctx) 
MF = cl.mem_flags 

dx = 0.01 
X = sp.arange(0.0, 100, dx) 
N = len(X) 

A_h = (sp.exp(-(X-50)**2/10.)*sp.exp(-1j*1000*X)).astype(sp.complex64) 
A_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=A_h) 

plot.figure() 
plot.plot(X, abs(A_h) ** 2/(abs(A_h) ** 2).max()) 

Source = """ 
    #define complex_ctr(x, y) (float2)(x, y) 
    #define complex_mul(a, b) complex_ctr(mad(-(a).y, (b).y, (a).x * (b).x), mad((a).y, (b).x, (a).x * (b).y)) 
    #define complex_unit (float2)(0, 1) 

    __kernel void propagate(__global float2 *A){ 
     const int gid_x = get_global_id(0); 
     float EPS = 0.1f; 
     A[gid_x] = A[gid_x] + EPS*complex_mul((A[gid_x+1] + A[gid_x-1]), complex_unit); 
    } 
""" 
prg = cl.Program(ctx, Source).build() 
for i in range(3000): 
    print i 
    event = prg.propagate(queue, (N,), None, A_d) 
    event.wait() 
cl.enqueue_copy(queue, A_h, A_d) 

plot.plot(X, abs(A_h) ** 2/(abs(A_h) ** 2).max()) 

plot.show() 

そして、ここでの結果は

デスクトップ結果である：

マック結果：

グリーンラインCORR繁殖後のガウスに対応し、ブルーラインは初期ガウスである

この問題はNVidia側で何が起こる可能性がありますか？私は私がこれを防ぐために重要なステップが欠落し、それがややによる運にMac上で実行されていることと思う

EDIT

これは、ユーザーの提案に基づいて取り組んでいる私の最終的なコードである

import scipy as sp 
import pyopencl as cl 
import matplotlib.pyplot as plot 

ctx = cl.create_some_context() 
queue = cl.CommandQueue(ctx) 
MF = cl.mem_flags 

dx = sp.float32(0.001) 
X = sp.arange(0.0, 100, dx).astype(sp.float32) 
N = len(X) 

A_h = (sp.exp(-(X-50)**2/10.)*sp.exp(1j*1000*X)).astype(sp.complex64) 
A_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=A_h) 
B_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=A_h) 

plot.figure() 
plot.plot(X, abs(A_h) ** 2/(abs(A_h) ** 2).max()) 

Source = """ 
    #define complex_ctr(x, y) (float2)(x, y) 
    #define complex_mul(a, b) complex_ctr((a).x * (b).x - (a).y * (b).y, (a).x * (b).y + (a).y * (b).x) 
    #define complex_unit (float2)(0, 1) 

    __kernel void propagate(__global float2 *A){ 
     const int gid_x = get_global_id(0); 
     float EPS = 0.1f; 
     float2 a1, a2; 
     a1 = A[gid_x-1]; 
     a2 = A[gid_x+1]; 
     barrier(CLK_GLOBAL_MEM_FENCE); 
     A[gid_x] += EPS*complex_mul((a1 + a2), complex_unit); 
    } 
""" 

prg = cl.Program(ctx, Source).build() 
for i in range(12000): 
    print i 
    evolve = prg.propagate(queue, (N,), None, A_d) 
    evolve.wait() 
cl.enqueue_copy(queue, A_h, A_d) 

plot.plot(X, abs(A_h) ** 2) 

plot.show() 

Answer 1

編集：ああ、ちょうど@talonmiesコメントを読んで、これと同じ解決策です。

このコードは、OpenCLの中に安全ではない、それはデータ競合の問題があります。

A[gid_x] = A[gid_x] + EPS*complex_mul((A[gid_x+1] + A[gid_x-1]), complex_unit); 

すべての作業項目xはx+1とx-1を使用しています。作業項目のスケジュールに応じて、結果は異なります。

使用2バッファの代わりに、Aから読み、簡単に、Bに書き込む：

import scipy as sp 
import pyopencl as cl 
import matplotlib.pyplot as plot 

ctx = cl.create_some_context() 
queue = cl.CommandQueue(ctx) 
MF = cl.mem_flags 

dx = 0.01 
X = sp.arange(0.0, 100, dx) 
N = len(X) 

A_h = (sp.exp(-(X-50)**2/10.)*sp.exp(-1j*1000*X)).astype(sp.complex64) 
A_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=A_h) 
B_d = cl.Buffer(ctx, MF.READ_WRITE) 

plot.figure() 
plot.plot(X, abs(A_h) ** 2/(abs(A_h) ** 2).max()) 

Source = """ 
    #define complex_ctr(x, y) (float2)(x, y) 
    #define complex_mul(a, b) complex_ctr(mad(-(a).y, (b).y, (a).x * (b).x), mad((a).y, (b).x, (a).x * (b).y)) 
    #define complex_unit (float2)(0, 1) 

    __kernel void propagate(__global float2 *A, __global float2 *B){ 
     const int gid_x = get_global_id(0); 
     float EPS = 0.1f; 
     B[gid_x] = A[gid_x] + EPS*complex_mul((A[gid_x+1] + A[gid_x-1]), complex_unit); 
    } 
""" 
prg = cl.Program(ctx, Source).build() 
for i in range(3000): 
    print i 
    event = prg.propagate(queue, (N,), None, A_d, B_d) 
    A_d, B_d = B_d, A_d #Swap buffers, so A always has results 
    #event.wait() #You don't need this, this is slowing terribly the execution, enqueue_copy already waits 
cl.enqueue_copy(queue, A_h, A_d) 

plot.plot(X, abs(A_h) ** 2/(abs(A_h) ** 2).max()) 

plot.show()