网友提问：我是一个cuda初学者，打算把现有的有限元程序加入cuda支持以提高运算速度。
有限元程序中，需要计算每个单元上的各个节点的内力分量，然后把这些节点力分量累加在一起。
比如，一个8节点长方体单元，可以计算出其8个节点的内力分量为
对每一个单元进行循环
LOOP KE=1 TO NELE
FORCE_ELE[1]=
FORCE_ELE[2]=
FORCE_ELE[3]=
…
FORCE_ELE[8]=
然后得到单元各个节点的编号
N1=NC[1][KE]
N2=NC[2][KE]
N3=NC[3][KE]
…
N8=NC[8][KE]
把这些力累加到节点上去
FORCE[N1]+=FORCE_ELE[1]
FORCE[N2]+=FORCE_ELE[2]
FORCE[N3]+=FORCE_ELE[3]
…
FORCE[N8]+=FORCE_ELE[8]
END LOOP
其中NC数组存储了网格的连接性， FORCE数组为gpu中全局数组。
因为单元网格是非结构化的，所以对于FORCE的读写就基本上是随机的（当然在执行效率上会非常差，这个以后再改进），所以在有多个BLOCK和THREAD的配置下执行时，应该会有同时读写同一块内存位置的时候。
当读FORCE的时候，CUDA能够保证其正确性，但是当写入时候，CUDA不能保证同时写入是否正确。
我也看了SKD中的histogram的例子，该例子通过用在每一个线程或wrap中分别存储结果后，然后再累加的方法解决了随机写会引起的冲突。
但是，在我的程序中，因为FORCE数组实在是太大，基本上应该在百万字节的大小，基本上不可能在每一个线程或wrap中分别全部存储。而且，实际上每个单元只有8个节点的数据需要更新，并没有必要分别每个线程或wrap中存储全部的FORCE数组。
这个问题，也许也有别人碰到过，请教有没有什么好地解决办法？
专家解答：
A:通过保留部分和再累加的方法避免写冲突的策略仍然是可行的。如果数据太多不能在shared memory中全部保存，一种简单的方法是launch多个kernel，在第一个kernel中进行局部数据更新（即更新FORCE_ELE)，然后把所有局部数据写到global memory里；然后用另一个kernel来根据局部数据更新全局的FORCE。从算法看起来FORCE对FORCE_ELE累加的顺序并没有要求，所以用两个kernel是可能解决这个问题的。
我根据上面有限元的步骤，按照类似的步骤写了一段模拟程序如下，请帮我看看该如何解决，多谢
程序：
#include <stdio.h>
#define BLOCK_NUM 3
#define THREAD_NUM 1
#define NUM_ELE 4
#define NUM_NOD 18
#define NC(j,N) *(NC+(N-1)*8+j-1)
#define FORCE(N) *(FORCE+N-1)
//global void init_force_kernel(float );
global void internal_kernel(int , float *);
int main(void)
{
int NC[8*NUM_ELE];
float FORCE[18];
NC(1,1)=1 ;
NC(2,1)=3 ;
NC(3,1)=9 ;
NC(4,1)=7 ;
NC(5,1)=2 ;
NC(6,1)=4 ;
NC(7,1)=10 ;
NC(8,1)=8 ;
NC(1,2)=3 ;
NC(2,2)=5 ;
NC(3,2)=11 ;
NC(4,2)=9 ;
NC(5,2)=4 ;
NC(6,2)=6 ;
NC(7,2)=12 ;
NC(8,2)=10 ;
NC(1,3)=7 ;
NC(2,3)=9 ;
NC(3,3)=15 ;
NC(4,3)=13 ;
NC(5,3)=8 ;
NC(6,3)=10 ;
NC(7,3)=16 ;
NC(8,3)=14 ;
NC(1,4)=9 ;
NC(2,4)=11 ;
NC(3,4)=17 ;
NC(4,4)=15 ;
NC(5,4)=10 ;
NC(6,4)=12 ;
NC(7,4)=18 ;
NC(8,4)=16 ;
int *dNC;
float *dFORCE;
cudaMalloc( (void**)&dNC,NUM_ELE8sizeof(int));
cudaMalloc( (void**)&dFORCE,NUM_NOD*sizeof(float));
cudaMemcpy(dNC,NC,NUM_ELE8sizeof(int),cudaMemcpyHostToDevice);
// initialize the force first
cudaMemset(dFORCE,0.0,NUM_NOD*sizeof(float));
internal_kernel < < <BLOCK_NUM,THREAD_NUM>>>(dNC,dFORCE);
cudaMemcpy(FORCE,dFORCE,NUM_NOD*sizeof(float),cudaMemcpyDeviceToHost);
for(int i=1;i <=NUM_NOD;i++)
printf(“NOD %d %f\n”,i,FORCE(i));
cudaFree(dNC);
cudaFree(dFORCE);
return 0;
}
global void internal_kernel(int* NC, float *FORCE)
{
const int tid=threadIdx.x;
const int bid= blockIdx.x;
const int THREADS=blockDim.x;
const int BLOCKS=gridDim.x;
float FORCE_ELE[8];
//compute internal force
for(int NELE=bidTHREADS+tid+1;NELE <=NUM_ELE;NELE+=BLOCKSTHREADS)
{
for(int j=1;j <=8;j++)
FORCE_ELE[j-1]=NELE10.0+j1.0;
for(int j=1;j <=8;j++)
{
int NOD;
NOD=NC(j,NELE);
FORCE(NOD)+=FORCE_ELE[j-1];
// printf(“NELE: %1d J:%1d NOD:%2d \t FORCE_ELE:%10.1f\t FORCE:%10.1f\n”,NELE,j,NOD,FORCE_ELE[j-1],FORCE(NOD));
}
}
return ;
}