2014-04-07 53 views
0

我有以下形式的示例循環。請注意,我的psi[i][j]依賴於psi[i+1][j], psi[i-1][j], psi[i][j+1] and psi[i][j-1],我只能計算內部矩陣的psi。現在,我嘗試在CUDA中編寫此代碼,但結果與順序不同。CUDA並行化依賴2D陣列

for(i=1;i<=leni-2;i++) 
for(j=1;j<=lenj-2;j++){ 
    psi[i][j]=(omega[i][j]*(dx*dx)*(dy*dy)+(psi[i+1][j]+psi[i-1][j])*(dy*dy)+(psi[i][j+1]+psi[i][j-1])*(dx*dx))/(2.0*(dx*dx)+2.0*(dy*dy)); 
} 

這是我的CUDA格式。

//KERNEL 
__global__ void ComputePsi(double *psi, double *omega, int imax, int jmax) 
{ 
int x = blockIdx.x; 
int y = blockIdx.y; 
int i = (jmax*x) + y; 
double beta = 1; 
double dx=(double)30/(imax-1); 
double dy=(double)1/(jmax-1); 

if((i)%jmax!=0 && (i+1)%jmax!=0 && i>=jmax && i<imax*jmax-jmax){ 
    psi[i]=(omega[i]*(dx*dx)*(dy*dy)+(psi[i+jmax]+psi[i-jmax])*(dy*dy)+(psi[i+1]+psi[i-1])*(dx*dx))/(2.0*(dx*dx)+2.0*(dy*dy)); 
} 
} 


//Code 
cudaMalloc((void **) &dev_psi, leni*lenj*sizeof(double)); 
cudaMalloc((void **) &dev_omega, leni*lenj*sizeof(double)); 
cudaMemcpy(dev_psi, psi, leni*lenj*sizeof(double),cudaMemcpyHostToDevice); 
cudaMemcpy(dev_omega, omega, leni*lenj*sizeof(double),cudaMemcpyHostToDevice); 
dim3 grids(leni,lenj); 
for(iterpsi=0;iterpsi<30;iterpsi++)   
    ComputePsi<<<grids,1>>>(dev_psi, dev_omega, leni, lenj); 

其中psi[leni][lenj] and omega[leni][lenj]和雙數組。

問題在於順序問題,CUDA代碼給出了不同的結果。代碼中是否需要修改?

+0

無關提示:我總是將輸入緩衝區聲明爲const指針,所以我不會搞砸。 – texasbruce

+0

相關提示:始終將輸入和輸出緩衝區分開,並且不要寫入全局輸入緩衝區 – texasbruce

回答

1

您正在全局內存中工作,並且您正在更改psi條目,而其他線程可能需要舊值。只需將新迭代的值存儲在單獨的變量中。但請記住,每次迭代後必須交換變量! 更復雜的方法是使用共享內存和空間域分配到獨立線程的解決方案。只需谷歌的CUDA教程解決熱/擴散方程,你會明白。

+0

感謝您的回覆,最初我在每次迭代後交換變量。但是由於我的迭代次數非常多,有人建議在cuda上交換需要很多時間,所以我正在尋找更快的替代方法。你的第二個解決方案看起來很有趣,我會看看它。謝謝。 – unkn0wn

+0

你不必真的交換。只需在你的內核中添加一個bool,並在每次迭代之後交換你的主機上的bool。那麼你的內核必須實現如 'if(swap){temp_psi = ... psi ...;} else {psi = ... temp_psi ...}'。 – OnWhenReady

+0

Helllo在那裏,我試圖尋找你第一個建議的交換方法,但可以找到任何解決方案值得我。您能否提出任何適合的解決方案或鏈接詳細說明使用交換變量來解決gauss-seidel問題。謝謝。 – unkn0wn

0
for(i=1;i<=leni-2;i++) 
    for(j=1;j<=lenj-2;j++){ 
    psi[i][j]= (omega[i][j]*(dx*dx)*(dy*dy) + 
       (psi[i+1][j]+psi[i-1][j]) * (dy*dy) + 
       (psi[i][j+1]+psi[i][j-1]) * (dx*dx) 
       )/(2.0*(dx*dx)+2.0*(dy*dy)); 
} 

我認爲,這個內核是不正確的順序之一:的psi[i][j]的值取決於操作的順序放在這裏 - 所以你將使用不更新psi[i+1][j]psi[i][j+1],但psi[i-1][j]psi[i][j-1]已更新在這個掃描

請確保使用CUDA的結果將會不同,其中操作順序不同。

爲了執行這樣的排序,如果可能的話,您需要插入如此多的同步,這可能對CUDA來說不值得。它真的是你需要做的嗎?

+0

使用這個內核的同步根本無濟於事。正如我所建議的OP應該使用交換變量。這種核函數(思想)是許多人(包括我)使用CUDA解決的微分方程的典型例子。並行解決PDE/ODE是CUDA確實存在的原因;-) – OnWhenReady

+0

是的,我認爲這是錯誤的順序 - 但是如果它真的是他打算做的......通常你需要編寫後向有限差分PDE內核爲'psi_new [i] [j] = f(psi)',然後用'psi_new'交換'psi'。並行和順序。讓我們看看OP如何看待它。 – Sigismondo

+0

這正是我說的話:-) – OnWhenReady