用一个block里的多个线程从globle memory往一个shared memory中的数组里拷数据，然后处理。由于数据太多，所以打算用个for循环处理。代码大致结构如下
for (i=0;i<n;i++)
{
1 把部分数据复制到shared memory数组里，
2 __syncthreads()同步，
3 数据处理，
4 __syncthreads()同步。
}
但在第3步处理数据的时候发现，对每个线程而言，只有这个线程所在warp在第1步中往shared memory数组里所复制的值是可见的，这说明__syncthreads()没有起作用。如果把for循环删掉就能获得正确的结果。这个问题有没有啥官方的说法，官方的编程文档貌似没提到这个问题。
另外，一个必须串行进行的操作（有上万步），是在host里用循环语句多次调用小kernel函数快，还是在直接在一个大kernel函数中用循环语句快？
对于第一个问题，for循环中调用同步函数应该是没有问题的，你需要留意是否一个block中所有线程都会运行到此同步函数。如果没有，则会出现错误。你可以运行下面代码进行验证。
/***************************************************
*main.cu
*DATE: 2012/12/13
*
*Test if the __syncthreads can do sync in for loop
**************************************************/
#include <stdio.h>
__global__ void Kernel(int* a)
{
__shared__ int s_a[32*3];
int index = blockIdx.x * blockDim.x + threadIdx.x;
for(int i=0; i < 3; i++)
{
/* copy data to shared memory */
s_a[index] = a[index + i*32*3];
/* sync */
__syncthreads();
/* print value */
printf("loop = %d, index = %d, a = %d\n", i, index, s_a[(index + 32)%(32 * 3)]);
/* sync */
__syncthreads();
}
}
int main(){
/* Malloc and copy data */
int *h_a, *d_a;
h_a = (int *)malloc(3*3*32*sizeof(int));
cudaError eFlag = cudaMalloc(&d_a, 3*3*32*sizeof(int));
for(int i = 0; i < 3*3*32; i++)
{
h_a[i] = i;
}
cudaMemcpy(d_a,h_a,3*3*32*sizeof(int),cudaMemcpyHostToDevice);
/* Call Kernel */
Kernel<<<1,32*3>>>(d_a);
cudaDeviceReset();
return 0;
}
第二个问题，如果程序过大一般来说分成较小的kernel会速度快些，如果由于串行操作过长导致内部消耗了过多的寄存器资源造成Occupation 的降低可能导致程序实际并行执行的线程减少。另外较小的Kernel有助于你来分析到底最耗时的操作在哪个Kernel中有助于接下来针对性的优化
赞一下2# zehuanwang 详尽的回答。
关于第二个问题，如果LZ那个串行执行的任务每一步都可以并行，只是每一步之间必须严格顺序，那么不断的发射kernel估计效果还可以。如果是只能一个线程使劲跑那种，也许让CPU做效果更好。这个问题如需深入，还请LZ提供更为详细的信息，如果算法上能改为可并行的，效果会好很多。
祝您编码愉快！
同步的问题终于解决了，但又出了新问题。实验中发现CUDA的舍入误差会导致结果波动，而且波动幅度高达30%。kenerl中所有的数学函数都改用了双精度函数也没有改善。C++串行代码中如果全部使用float类型虽然也存在舍入误差，但误差远远小于CUDA的结果，而且误差基本恒定，不像CUDA中那样变来变去。不知这个问题有没有解决的方法。
CUDA的舍入误差没有那么高的，而且换成双精度也没改善更说明了这不是舍入精度的问题。
估计您的代码可能还存在BUG，请尽量寻找一下。
做的是一个Cholesky矩阵分解，计算法公式如下
for k = 1 : n
A(k, k) = sqrt(A(k, k))
for i = k + 1 : n
A(i, k) = A(i, k)/A(k, k)
end
for j = k + 1 : n
for i = j : n
A(i, j) = A(i, j)− A(i, k)A(j, k)
end
end
end
从这个公式可以看出，分解得到的三角矩阵是从左到右一列一列计算的，越靠近右侧的元素会经历越多的求和操作。
当一个矩阵的条件数较小时，所增加的填充元也较少，这意味这右侧元素所经历的求和操作越少，在这种情况下我代码的结果是正确并且稳定的。当矩阵条件数较大时，则右侧元素所经历的求和操作越来越多，此时我代码结果的波动从左到右依次增加。一步一步调试也能发现每一步相比串行计算都增加了一个很微小的随机误差，看起来很像是累加误差，但奇怪的是我用了kahan求和之后结果并没有改善。所以现在有点点怀疑是不是CUDA里的sqrt（）函数不准确，待会用一个线程完全串行的跑一次试试。
根据programming guide （CUDA 4.2）的Table C-1，sqrtf()函数的maximum ULP error在计算能力>=2并且编译参数为 -prec-sqrt=true时，为0；其他情况为3。
Table C-2，sqrt()函数 maximum ULP error 为0（IEEE754 round-to-nearest-even）。
其他一些sqrt变种也都是IEEE 754标准相容的。
所以说这个精度应该是可以的。
不过您的CPU实现里面的浮点运算，有可能是跑的80X87的模式，中间过程是80bit精度的，那个精度应该更好一些。
总之，还需要您进一步研究才行。
祝您编码愉快~
弄了好几天，终于找到原因了，但问题还是没解决。分解之所以出现偏差是由于前面所得到的矩阵中的微小偏差被分解放大所带来的。前面的矩阵是有限元系统矩阵，有限元系统矩阵是由一个个小的单元系统矩阵累加所得到的，而累加过程中会舍弃掉一些末尾的有效数字，由于CUDA中个线程的执行顺序是随机的，单元系统矩阵累加的顺序也是随机的，因此舍弃掉的尾数值也是随机的，最后会造成末尾有效数字上会出现波动。这个问题理论上应该可以通过kahan求和改善，但由于累加过程中必须用到原子操作，kahan求和没法应用。
有限元算法的该累加过程依赖于累加顺序么？如何说明CPU累加的顺序就是正确的？
理论上来说累加的结果和累加的顺序无关，但计算机中都存在一个舍入误差的的问题（在累加过程中，最低位的有效数字被舍去，节省出来的有效数字位数被用于表示累加结果中的新出现的高位值），考虑到这个因素，累加的结果是和累加的顺序是有关的。在CPU中，由于累加的顺序是确定的，因此累加的结果也是确定的。CPU计算中后来我加上了kahan求和来改善累加的精度，最后的结果和matlab double精度的计算结果基本一致。现在在尝试通过对单元的分组来避免原子操作，然后运用kahan求和。
好的，多谢您的反馈。
某来源建议了一种本地实现kahan修正项的计算方法，但未经验证，我大致叙述如下，供您参考。
根据programming guide的说法，atomicAdd()函数的返回值是该累加地址累加前的值。我们试图利用这一特性来生成每个线程自己的kahan修正项。
float old=atomicAdd (&total_sum,a); //old 拿到total_sum原子加以前的值
float local_sum=old+a; //在本地再做一次加法，假定本地加法和原子加法精度特性相同
float c=a-(local_sum-old); //c即是本线程kahan算法的修正项
这样每个线程都保留了kahan算法的修正项（这里举例是进行一次kahan相加，多次应该也可以，需要略作改动），只是最后还要将这些修正项加在一起才行，考虑到这些修正项都比较小，不易发生大数吃小数的问题。如果是CPU循环执行，每次可以把上一步的积累的修正项更新一下，似乎要更好一点。也可以考虑将这些修正项选若干个规约求和以后，再用上述方法加给total_sum，并得到一个新的修正项，等等。
此方法未经验证，仅供您参考。
本文参考了http://cudazone.nvidia.cn/forum/forum.php?mod=viewthread&tid=1275&page=1
以及和樟树在QQ交谈中得到的建议。
祝您编码愉快~
结果不稳定的问题最终通过我前面提到的方法解决了，即将单元分成没有组内不共点的几个组，然后一个组一个组的累加，这样就不需要原子操作，可以用原版的kahan算法了。
另外关于循环与同步问题还作了个小实验，结果貌似是这样的：
单重循环内部的__syncthreads()是有效的，即
for (i=0;i<n;i++)
{
操作1
__syncthreads()；
操作2
__syncthreads()；
}
这样是可以的。
但嵌套循环中内层循环内部的__syncthreads()是无效的，即
for (i=0;i<n;i++)
{
for (j=0;i<m;j++)
{
操作1
__syncthreads()；
操作2
__syncthreads()；
}
}
这样的同步操作无效。
不知我得出的结论是否正确。
感谢您的反馈，我也不清楚多层嵌套循环内部线程同步有没有什么限制，期待原厂支持进一步揭示此问题。
祝您编码愉快~
我把2#的代码简单改成二重循环试验了一下，结果是符合预期的，没有发现异常。所以12#给出的循环内的同步应该是可用的，如有问题可能在其他方面。
#include <stdio.h>
__global__ void sync_Kernel(int *a)
{
__shared__ int s_a[32*3];
int tid = threadIdx.x+blockIdx.x*blockDim.x;
for (int j=0;j<4;++j)
{
for (int i=0;i<3;++i)
{
//copy data to shared memory
s_a[tid]=a[tid+i*32*3+j*32*3*3];
__syncthreads();
//print value
//printf("j=%d, i=%d, tid=%d, a[%d]=%d \n",j,i,tid,(tid+i*32*3+j*32*3*3),s_a[tid]);
if (tid%32==0)
printf("j=%d, i=%d, tid=%d, a=%d \n",j,i,tid,s_a[(tid+32)%(32*3)]);
__syncthreads();
}
if(tid==0)
printf("\n");
}
return;
}
int main()
{
//malloc
int *h_a=NULL,*d_a=NULL;
h_a=(int*)malloc(32*3*3*4*sizeof(int));
cudaMalloc(&d_a,32*3*3*4*sizeof(int));
//memset
cudaMemset(d_a,0,32*3*3*4*sizeof(int));
memset(h_a,0,32*3*3*4*sizeof(int));
for (int i=0;i<32*3*3*4;++i)
{
h_a[i]=i;
}
cudaMemcpy(d_a,h_a,32*3*3*4*sizeof(int),cudaMemcpyHostToDevice);
//invoke kernel
sync_Kernel<<<1,32*3>>>(d_a);
cudaDeviceReset();
free(h_a);
cudaFree(d_a);
h_a=NULL;
d_a=NULL;
return 0;
}
祝您编码愉快