__global__ void KernelV(int *dev_XMatrix,double *dev_V,double *dev_VK)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bid*blockDim.x+tid;
if(xIndex<XROW)
{
if(xIndex==0)
{
for(int j=0;j<XCOL;j++)
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[j+xIndex*XCOL]*dev_VK[j];
__syncthreads();
}
}
else
{
for(int j=0;j<XCOL;j++)
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[j+xIndex*XCOL]*dev_VK[j];
__syncthreads();
dev_V[xIndex*XROW]=dev_V[xIndex*XROW]+dev_XMatrix[j+xIndex*XCOL]*dev_XMatrix[j]*dev_VK[j];
}
__syncthreads();
}
}
}
__global__ void KernelV2(int *dev_ZZ,double *dev_V,double *dev_VKK)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bid*blockDim.x+tid;
if(xIndex<XROW)
{
if(xIndex==0)
{
for(int j=0;j<(XCOL*(XCOL-1)/2);j++)
{
dev_V[xIndex]=dev_V[xIndex]+dev_ZZ[j*XROW]*dev_ZZ[xIndex+j*XROW]*dev_VKK[j];
__syncthreads();
}
}
else
{
for(int j=0;j<(XCOL*(XCOL-1)/2);j++)
{
dev_V[xIndex]=dev_V[xIndex]+dev_ZZ[j*XROW]*dev_ZZ[xIndex+j*XROW]*dev_VKK[j];
__syncthreads();
dev_V[xIndex*XROW]=dev_V[xIndex*XROW]+dev_ZZ[xIndex+j*XROW]*dev_ZZ[j*XROW]*dev_VKK[j];
}
__syncthreads();
}
}
}
__global__ void KernelVOne(double *dev_V,double dev_VK,int *dev_XMatrix,int j,double s0)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bid*blockDim.x+tid;
if(xIndex<XROW)
{
if(xIndex==0)
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[xIndex*XCOL+j]*(dev_VK-s0);
__syncthreads();
}
else
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[xIndex*XCOL+j]*(dev_VK-s0);
__syncthreads();
dev_V[xIndex*XROW]=dev_V[xIndex*XROW]+dev_XMatrix[xIndex*XCOL+j]*dev_XMatrix[j]*(dev_VK-s0);
__syncthreads();
}
}
}
int main(void)
{
culaStatus status;
status = culaInitialize();
cudaError_t cudaStatus;
clock_t t1,t2,t3,t4,t5,t6,t7,t8,t9; //设定计时器变量
readGen_sim();
readPhe_sim();
double *V=(double*)malloc(sizeof(double)*XROW*XROW);
for(int i=0;i<XROW*XROW;i++)
{
V=0;
}
double VK[XCOL];
for(int i=0;i<XCOL;i++)
{
VK=1e-30;
}
t1=clock(); //设定计时器变量t1开始计时
int *dev_XMatrix;
double *dev_V,*dev_VK;
cudaMalloc((void**)&dev_XMatrix,XSIZE);
cudaMalloc((void**)&dev_V,XROW*XROW*sizeof(double));
cudaMalloc((void**)&dev_VK,XCOL*sizeof(double));
cudaMemcpy(dev_XMatrix,XMatrix,XROW*XCOL*sizeof(int),cudaMemcpyHostToDevice);
cudaMemcpy(dev_V,V,XROW*XROW*sizeof(double),cudaMemcpyHostToDevice);
cudaMemcpy(dev_VK,VK,XCOL*sizeof(double),cudaMemcpyHostToDevice);
KernelV<<<50,300>>>(dev_XMatrix,dev_V,dev_VK);
cudaMemcpy(V,dev_V,XROW*XROW*sizeof(double),cudaMemcpyDeviceToHost);
cudaFree(dev_VK);
cudaFree(dev_XMatrix);
double *VKK=(double*)malloc(sizeof(double)*XCOL*XCOL);
for(int i=0;i<XCOL*XCOL;i++)
{
VKK=1e-30;
}
int *ZZ=(int*)malloc(sizeof(int)*XROW*(XCOL*(XCOL-1)/2));
int uu=0;
for(int i=0;i<(XCOL-1);i++)
{
for(int j=i+1;j<XCOL;j++)
{
for(int k=0;k<XROW;k++)
{
uu=((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+k;
ZZ[uu]=XMatrix[i+k*XCOL]*XMatrix[j+k*XCOL];
}
}
}
double *dev_VKK;
int *dev_ZZ;
cudaMalloc((void**)&dev_VKK,XCOL*XCOL*sizeof(double));
cudaMalloc((void**)&dev_ZZ,sizeof(int)*XROW*(XCOL*(XCOL-1)/2));
cudaMemcpy(dev_V,V,XROW*XROW*sizeof(double),cudaMemcpyHostToDevice);
cudaMemcpy(dev_VKK,VKK,XCOL*XCOL*sizeof(double),cudaMemcpyHostToDevice);
cudaMemcpy(dev_ZZ,ZZ,sizeof(int)*XROW*(XCOL*(XCOL-1)/2),cudaMemcpyHostToDevice);
KernelV2<<<50,300>>>(dev_ZZ,dev_V,dev_VKK);
cudaMemcpy(V,dev_V,XROW*XROW*sizeof(double),cudaMemcpyDeviceToHost);
cudaFree(dev_ZZ);
cudaFree(dev_VKK);
cudaFree(dev_V);
for(int i=0;i<XROW;i++)
{
V[i*XROW+i]=V[i*XROW+i]+v0;
}
double *VI=(double*)malloc(sizeof(double)*XROW*XROW);
for(long int i=0;i<XROW*XROW;i++)
{
VI=V;
}
t2=clock(); //设定计时器变量t2开始计时
inverse(VI,XROW);
t3=clock(); //设定计时器变量t3开始计时
int iter=0;
double err=1000;
double R[XROW];
FILE *p1;
p1=fopen("iter_simulat.txt","w");
fprintf(p1,"iter error b resid_var\n");
cudaMemcpy(dev_XMatrix,XMatrix,XSIZE,cudaMemcpyHostToDevice);
while ((iter<iter_max) && (err>err_max)) //start the big loop
{
t4=clock(); //设定计时器变量t4开始计时
iter=iter+1;
double VK0[XCOL];
for(int i=0;i<XCOL;i++)
{
VK0=VK;
}
double VISUM=0.0;
double VcolSUM[XROW];
double VcolYPhe=0.0;
for(long int i=0;i<XROW*XROW;i++)
{
VISUM=VISUM+VI;
}
for(int i=0;i<XROW;i++)
{
VcolSUM=0.0;
}
for(int j=0;j<XROW;j++)
{
for(int i=0;i<XROW;i++)
{
VcolSUM[j]=VcolSUM[j]+VI[j+i*XROW];
}
}
for(int i=0;i<XROW;i++)
{
VcolYPhe=VcolYPhe+VcolSUM*YPhe;
}
b=(1/VISUM)*VcolYPhe; //update the fixed effects
//printf("%.7lf\n",b);
double s0=v0;
for(int i=0;i<XROW;i++)
{
R=YPhe-b;
}
double RTVI[XROW];
double RTVIRSUM=0.0;
for(int i=0;i<XROW;i++)
{
RTVI=0.0;
}
for(int i=0;i<XROW;i++)
{
for(int j=0;j<XROW;j++)
{
RTVI=RTVI+R[j]*VI[i+XROW*j];
}
}
for(int i=0;i<XROW;i++)
{
RTVIRSUM=RTVIRSUM+RTVI*R;
}
v0=(s0*RTVIRSUM)/XROW;
for(int i=0;i<XROW;i++)
{
V[i*XROW+i]=V[i*XROW+i]+(v0-s0);
}
for(int i=0;i<XROW*XROW;i++)
{
VI=V;
}
inverse(VI,XROW); //update the residual variance
t5=clock(); //设定计时器变量t5开始计时
cudaMemcpy(dev_V,V,XROW*XROW*sizeof(double),cudaMemcpyHostToDevice);
for(int j=0;j<XCOL;j++) //start update main effect QTL variance
{
for(int ione=0;ione<XROW*XROW;ione++)
{
VItemp1[ione]=VI[ione];
}
s0=VK[j];
double ZVZ=0.0;
double YVZ=0.0;
double ZTVI[XROW];
for(int i=0;i<XROW;i++)
{
ZTVI=0.0;
}
for(int i=0;i<XROW;i++)
{
for(int k=0;k<XROW;k++)
{
ZTVI=ZTVI+XMatrix[j+k*XCOL]*VI[i+k*XROW];
}
}
for(int i=0;i<XROW;i++)
{
ZVZ=ZVZ+ZTVI*XMatrix[j+i*XCOL];
}
for(int i=0;i<XROW;i++)
{
RTVI=0.0;
}
for(int i=0;i<XROW;i++)
{
for(int k=0;k<XROW;k++)
{
RTVI=RTVI+R[k]*VI[i+k*XROW];
}
}
for(int i=0;i<XROW;i++)
{
YVZ=YVZ+RTVI*XMatrix[j+i*XCOL];
}
double dum=0.0; //using explicit algorithm
dum=Maxmium((YVZ*YVZ-ZVZ)/(ZVZ*ZVZ)+s0,1e-30);
VK[j]=dum;
/*for(int i=0;i<XROW;i++)
{
for(int k=0;k<XROW;k++)
{
V[i*XROW+k]=V[i*XROW+k]+XMatrix[i*XCOL+j]*XMatrix[k*XCOL+j]*(VK[j]-s0);
}
}*/
KernelVOne<<<50,300>>>(dev_V,VK[j],dev_XMatrix,j,s0);
cudaStatus = cudaDeviceSynchronize();
if (cudaStatus != cudaSuccess)
{
fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching KernelV11!\n", cudaStatus);
}
cudaMemcpy(V,dev_V,XROW*XROW*sizeof(double),cudaMemcpyDeviceToHost);
for(int i=0;i<XROW*XROW;i++)
{
VI=V;
}
if((VK[j]-s0)!=0)
{
inverse(VI,XROW);
for(int ione=0;ione<XROW*XROW;ione++)
{
VItemp1[ione]=VI[ione];
}
}
else
{
for(int ione=0;ione<XROW*XROW;ione++)
{
VI[ione]=VItemp1[ione];
}
}
} //end update main effect QTL variance
t6=clock(); //设定计时器变量t6开始计时
for(int i=0;i<XCOL-1;i++) //start update epistatic effect QTL variance
{
for(int j=i+1;j<XCOL;j++)
{
for(int itwo=0;itwo<XROW*XROW;itwo++)
{
VItemp2[itwo]=VI[itwo];
}
s0=VKK[j+i*XCOL];
double ZVZ=0.0;
double YVZ=0.0;
double ZTVI[XROW];
for(int n=0;n<XROW;n++)
{
ZTVI[n]=0.0;
}
for(int k=0;k<XROW;k++)
{
for(int n=0;n<XROW;n++)
{
ZTVI[k]=ZTVI[k]+ZZ[((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+n]*VI[k+n*XROW];
}
}
for(int n=0;n<XROW;n++)
{
ZVZ=ZVZ+ZTVI[n]*ZZ[((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+n];
}
for(int n=0;n<XROW;n++)
{
RTVI[n]=0.0;
}
for(int k=0;k<XROW;k++)
{
for(int n=0;n<XROW;n++)
{
RTVI[k]=RTVI[k]+R[n]*VI[k+n*XROW];
}
}
for(int n=0;n<XROW;n++)
{
YVZ=YVZ+RTVI[n]*ZZ[((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+n];
}
double dum=0.0; //using explicit algorithm
dum=Maxmium((YVZ*YVZ-ZVZ)/(ZVZ*ZVZ)+s0,1e-30);
VKK[j+i*XCOL]=dum;
for(int n=0;n<XROW;n++)
{
for(int k=0;k<XROW;k++)
{
V[n*XROW+k]=V[n*XROW+k]+ZZ[((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+n]*ZZ[((2*XCOL*i-i*i-i)*XROW/2)+(j-i-1)*XROW+k]*(VKK[j+i*XCOL]-s0);
}
}
for(int n=0;n<XROW*XROW;n++)
{
VI[n]=V[n];
}
if((VKK[j+i*XCOL]-s0)!=0)
{
inverse(VI,XROW);
for(int itwo=0;itwo<XROW*XROW;itwo++)
{
VItemp2[itwo]=VI[itwo];
}
}
else
{
for(int itwo=0;itwo<XROW*XROW;itwo++)
{
VI[itwo]=VItemp2[itwo];
}
}
}
} //end update epistatic effect QTL variance
t7=clock(); //设定计时器变量t7开始计时
double v01=v0;
double VK1[XCOL];
for(int n=0;n<XCOL;n++)
{
VK1[n]=VK[n];
}
err=0.0;
for(int n=0;n<XCOL;n++)
{
err=err+(VK0[n]-VK1[n])*(VK0[n]-VK1[n]);
}
err=err/XCOL;
// printf("%2d %.10lf %.7lf %.7lf\n",iter,err,b,v01);
fprintf(p1,"%2d %e %e %e\n",iter,err,b,v01);
}
下面是我写的一个小程序，在逐步实现并行的过程中遇到一个问题，就是我写了3个KERNEL，第3个KERNEL不能线程全部同步，如果仅仅有前两个KERNEL，运算结果和串行是一摸一样的，但是加了第3个KERNEL（KernelVOne）就出错了，实在是找不出原因。我能想到的错误：1.线程没有同步，线程没有全部结束；2.数组越界；3.参数没有传到设备端；4.有关设备端内存分配、传递、KERNEL的位置放的不对。第2、3个错误应该不会发生，因为我觉得自己检查的挺认真，现在我认为最可能的就是第1、4个错误。其实，每一个KERNEL我都单独拿出来试过，跟SAS运算结果一样，所以我认为KERNEL本身应该没有错误，其实就是把for循环展开了而已，我的KERNEL代码上面屏蔽的代码就是它的串行代码。为了便于查看结果的正确与否，我把真实数据缩小了很多，后面的代码也有好多地方可以展开，但是第3个KERNEL我就卡住了，没有办法再去写第4个第5个。。。，这个月底老板就要汇报了，着实的着急，这个问题已经弄了第3天了，还请各位耐心帮帮我，感激不尽！！！
以下是我的头文件，我知道发这么长的代码很烦人，不过我确实是想不到其他办法了，还请版主理解和见谅！
#include “cuda_runtime.h”
#include “device_launch_parameters.h”
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cula_lapack.h>
#include <cublas_v2.h>
#include “time.h”
#define XSIZE (sizeof(int)*5000)
#define XROW 100
#define XCOL 50
#define iter_max 50
#define err_max 1e-8
1：虽然代码模式要比之前清晰很多，但是400多行的代码，基本无注释，完全缩进，还大把的循环和循环嵌套，实在无法目测出哪里有问题。
2：您的kernelVOne 目测未发现明显问题，在xIndex！=0的情况下，有两次赋值，其中第二次是否会越界？是否是您算法所需要的行为，这请您自己保证。按照您的说法，并无越界。
3：关于您错误的原因，您提到4点。第一点是线程没有同步或线程没有全部结束。请问您是如何得出此结论的？如果需要block内的线程同步，请根据需要添加__syncthreads().如果要grid内线程全部同步，只能结束kernel。以及如果您是使用回拷的数据判断，这并不存在线程没有结束的问题。
4：您提出第二点，数组越界，并自己否认了这一点。这需要您自己根据算法实现来保证，我无法评论。
5：您提出第三点，参数没有传到设备端。请问依据何在？
6：您提出第四点，有关设备端内存分配、传递、KERNEL的位置放的不对。内存分配和copy数值需要您自己保证使用参数的正确性。如果参数是正确的，机器不会出错的。而对于“kernel的位置放得不对”，我无法理解。您是说代码位置中kernel启动时机不对么？这是需要您自己保证的。还是说在GPU上kernel执行的位置不对？这一点，你是无法干预的，而且一般也是不会出问题的。
7：您没有提供任何详细的出错信息和调试信息，所以无法回答您哪里有问题。
8：您提出每个kernel单独测试都是好的，那么根据逻辑推断，自然最大可能就是调用kernel的逻辑有问题，您可以根据这一点继续检查。当然并不能排除单个kernel有之前考虑不周的地方。您还是需要检查kernel和调用kernel的逻辑，有没有背叛您本来的意思。
9：您只说第三个kernel“卡住”，并没有提供详细的调试信息或者报错信息。这样回答的人也会被“卡住”。
最后，希望您提供更详尽的出错和调试信息，提供分段和缩进清晰，有注释的代码。这样才更方便于检查问题。您赶时间的心情可以理解，但是查错还得按照客观规律办事。
祝您编码愉快。
是这样的。我描述的可能不是很清楚，visual studio平台下面并没有提示错误信息，能够通过编译并且运行，但是输出到.txt文件里的结果是错误的。
对于您提到的3：，我是昨晚去计算机学院（我是农学院的）请教了一个做GPU加速方向的同学，他说他之前写代码的时候，曾经出现过前两个KERNEL正确执行，而第三个KERNEL出现问题了，他说KERNEL是没有错误的，原因是线程没有结束，然后我就想我的是不是也是因为线程没有结束？他告诉我在KERNEL后面加上 cudaDeviceSynchronize();
对于您提到的5：，我是通过在__global__函数里面printf()传递的值，我发现当我前两个KERNEL存在的话，比如说j输出的都是错误的，当我把前两个KERNEL屏蔽了，j的值输出是正确的。
对于您提到的9：，我说的“卡住”并不是不能执行，我说的意思是我把串行代码用第三个KERNEL代替，就不能输出正确结果了。如果三个KERNEL全部存在的话，屏幕刷屏，不断输出cudaDeviceSynchronize returned error code 30 after launching KernelV11!
非常感谢您的耐心回答，我再看看KERNEL逻辑关系吧。
1：编译阶段并不能发现所有的错误，编译通过的程序也可能在运行时出现各种错误。而如果是您本身编写程序逻辑上的问题，那么计算机按照代码执行完毕之后，结果依然可以是错的。
2：在kernel后面加上cudaDeviceSynchronize()可以保证此前的所有device操作都结束，可以用来在kernel后面等待kernel执行结束，然后查看返回值看有什么问题没有。但你在kernel后面有cudamemcpy，这其实保证了kernel要执行完才能回拷数据。（当然也不排除kernel直接挂掉，cuda里面一次kernel执行挂了，后面的所有device操作都会挂掉）
3：您曾经专门发帖询问传递参数的问题，我有2楼回复，未见您回复反馈，不知对此处可有帮助？
4：您在最后才给出了最有价值的DEBUG信息。您代码本来就有用于debug的判断，但您没有在第一时间透露给我们。反复刷屏输出这样的错误信息，正表明了您反复启动的kernel已经挂掉了，也许在他之前就挂掉了（如果之前的每次kernel启动后面都有cudadeviceSynchronize，并且返回值都没问题的话，那么就是您循环这里的kernel挂了）。您可以根据返回的错误值用一个 cudageterrorstring神马的函数得到具体的错误类型（当然也可能是error unknow），然后看看。
既然这里kernel已经挂了，那么cudamemcpy也不会执行，自然不可能在host端得到正确的结果。
综上，建议您检查第三个kernel和调用该kernel的代码，看看有什么问题，先保证所有kernel的执行都不会报错。能从device端正常把结果拷回host端。
之后host端结果如果还有问题，再考虑您实现的算法逻辑上，哪里出现了偏差。
大致如上，祝您编码愉快~
对于您2楼的回复我检查过了，应该没有问题。其实我第三个KERNEL分离出来就是要执行如下的功能，这是一个简化版，应该比较清晰。
__global__ void KernelV11(double *dev_V,double dev_VK,int *dev_XMatrix,int j,double s0)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bid*blockDim.x+tid;
if(xIndex<XROW)
{
if(xIndex==0)
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[xIndex*XCOL+j]*(dev_VK-s0);
__syncthreads();
}
else
{
dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[xIndex*XCOL+j]*(dev_VK-s0);
__syncthreads();
dev_V[xIndex*XROW]=dev_V[xIndex*XROW]+dev_XMatrix[j]*dev_XMatrix[xIndex*XCOL+j]*(dev_VK-s0);
__syncthreads();
}
}
__syncthreads();
}
double Maxmium(double a, double b) //比较大小
{
if (a>b)
return a;
else
return b;
}
int main(void)
{
clock_t t1,t2;
double *V=(double*)malloc(XROW*XROW*sizeof(double));
int *XMatrix=(int*)malloc(XROW*XCOL*sizeof(int));
double *VK=(double*)malloc(XCOL*sizeof(double));
for(int i=0;i<XCOL;i++)
{
VK[i]=0.5;
}
for(int i=0;i<XROW*XCOL;i++)
{
XMatrix[i]=1+rand()%Mrand;
}
for(int i=0;i<XROW*XROW;i++)
{
V[i]=0.0;
}
double s0=0.1;
double *dev_V;
int *dev_XMatrix;
cudaMalloc((void**)&dev_V,XROW*XROW*sizeof(double));
cudaMalloc((void**)&dev_XMatrix,XROW*XCOL*sizeof(int));
//cudaMalloc((void**)&dev_VK,XCOL*sizeof(double));
t1=clock();
int k=0;
while(k<5)
{
for(int i=0;i<XROW;i++)
{
V[i*XROW+i]=V[i*XROW+i]+s0;
}
cudaMemcpy(dev_V,V,XROW*XROW*sizeof(double),cudaMemcpyHostToDevice);
//cudaMemcpy(dev_VK,VK,XCOL*sizeof(double),cudaMemcpyHostToDevice);
cudaMemcpy(dev_XMatrix,XMatrix,XROW*XCOL*sizeof(int),cudaMemcpyHostToDevice);
for(int j=0;j<XCOL;j++)
{
double dum=0.0;
dum=Maxmium(s0,0.01*j);
VK[j]=dum;
KernelV11<<<50,300>>>(dev_V,VK[j],dev_XMatrix,j,s0);
}
cudaMemcpy(V,dev_V,XROW*XROW*sizeof(double),cudaMemcpyDeviceToHost);
k++;
}
t2=clock();
for (int i=0;i<10;i++)
{
for(int j=0;j<6;j++)
{
printf("%f ",V[6*i+j]);
}
printf("\n");
}
printf("GPU所用的时间为：%d毫秒\n",(t2-t1));
cudaFree(dev_V);
//cudaFree(dev_VK);
cudaFree(dev_XMatrix);
free(V);
free(XMatrix);
free(VK);
}
这段代码跑出的结果跟SAS软件是一样的，说明KERNEL函数应该没有问题，KernelV11<<<50,300>>> 在那段400多行的代码中的运行跟这段代码类似，您可以给XROW,XCOL赋一下具体值，运行一下看看是不是有什么是我自己没发现的错误，因为我不会用工具调试，通常我喜欢printf（），辛苦您了！
重新目测了7#的代码，然眼力有限，实无法代替调试，同时也无法理解您的算法模型，无法给出进一步的建议了。
依然建议您先解决每步都报错的问题，使得程序无错运行，之后如有错误结果，再考虑逻辑问题。只能建议到这一步了，欢迎其他人补充。
祝您编码愉快~
我建议你可以先在内核执行后面全部加上:
printf(“%s\n”, cudaGetErrorString(cudaGetLashError()));
来确定内核执行是否出错。
好的，我想告诉大家一个不幸的消息，我写的内核函数有问题，只不过刚开始因为赋了一些特殊值，结果没有影响~非常感谢！
好像是cudaGetLastError(),我在每个KERNEL后面加上这句话之后，屏幕上全部是unkown error! 没有提供有用的信息
其实我一直怀疑您的kernel里面第二个赋值不代表您的真实意图，并且会写越界，但是您一直否认这一点。
可否检查下kernel里面这个地方呢？
祝您编码愉快~
好的，我再检查一下
我已经将代码修改如下所示，前两个KERNEL如果不屏蔽，第三个KERNEL就会出错，但是将前两个KERNEL屏蔽了，下面这个KERNEL代替串行部分就不会报错了！GOOGLE了好久，没发现我想要的答案，现在这个问题没有想明白，我现在只能将前两个KERNEL屏蔽了，因为它们可能不如第三个KERNEL节省的时间多，第三个在循环里面，所用的时间应该是单个KERNEL*循环次数，如何能让三个KERNEL同时存在呢？难道真的是越界问题？我再检查一下吧
__global__ void KernelVOne(double *dev_V,double dev_VK,int *dev_XMatrix,int j,double s0)
{
int k=threadIdx.x;
int i=blockIdx.x;
if((i<XROW)&&(k<XROW))
{
dev_V[i*XROW+k]=dev_V[i*XROW+k]+dev_XMatrix[i*XCOL+j]*dev_XMatrix[k*XCOL+j]*(dev_VK-s0);
}
}
没看明白您想表达什么？
1：CUDA并没有限制说只能两个kernel，第三个就一定会报错，所以必然还是您的程序有问题。
2：根据CUDA的规定，只要有一个kernel/API函数出错，后面就会一直出错。所以也可能是第三个kernel之前的某处有问题，建议您一步步地检查返回值。
3：什么叫“如何能让三个kernel同时存在”？是指都用上么？如果是的话，您把代码全部改对就可以。
4：以及希望您能在确实理解算法和CUDA实现的基础上，写好代码，找出BUG，解决问题。而不是头疼医头脚疼医脚，遇到问题赶紧问个大概，然后赶紧代进去，然后粗略地测试一下，然后报告成功或失败。这里面每一步都没有在坚实的基础上，不但会使得外人无法判断，也将使得您的问题变得越来越没有头绪。以上面第三个kernel为例，我给您的实现和您之前自己的实现，绝对是不等价的。但是您却表示，两者单独测试都能得到正确的结果，这一点我深表怀疑。
最后，麻烦您写汉字的时候也稍微排下版吧，无缩进+无标号+无空格看上去十分费劲啊。
祝您编码愉快。
好的。
就是向您第3：，说的那样，都用上；
对于您说的第4：，我解释一下，因为月底就汇报了，自己本身也是菜鸟，所以感觉没有那么多时间去看一些CUDA方面的资料，等我汇报完了，我认真看一下官方提供的一些资料，为了方便我是把赵开勇等人编的《高性能运算之CUDA》看了一遍。另外，是我没有表述明白，您给我的KERNEL和我原来的KERNEL结果是不一样的，您的结果是正确的，而我的结果是错误的。
最后，对于一些很唐突的问题，我深表歉意，很多时候是我力不从心，自己解决不了，汇报又在即，所以我就直接到论坛里发问了。如果我自己倒腾恐怕汇报的时候很难有什么进展，然后我还是再问一下吧，吃晚饭的时候我一直在想如果我的XROW>1024，那该如何应对呢？现在真是急了，还要做PPT，因为Fermi构架下，block中的线程不能大于1024，我的k=threadIdx.x,那样的话该采用什么办法来利用线程呢？比如说XROW=60000.
忙完这段时间，我会静下心来好好阅读一些CUDA方面的文献，体会和理解一下如何运用。
谢谢！
17#可能言重了，抱歉。
对于XROW==60000的情况可以这样。
假定每个block是512线程，一维。
对grid进行调整，使用2维的grid，将grid的x dimension和threadIdx一起来表示k；用grid的y dimension表示i。
具体为：
dim3 blocknum (60000/512+1,60000); //blockIdx.x is 0~60000/512,blockIdx.y is 0~59999
dim3 threadmun (512); //each block has 512 threads,threadIdx.x is 0~511
__global__ void BStatistics_3 (...)
{
int k=blockIdx.x*512+threadIdx.x;
int i=blockIdx.y;
if((i<XROW)&&(k<XROW))
{
dev_V[i*XROW+k]=dev_V[i*XROW+k]+dev_XMatrix[i*XCOL+j]*dev_XMatrix[k*XCOL+j]*(dev_VK-s0);
}
}
int main()
{
......
BStatistics_3 <<<blocknum,threadmun>>> (...)
......
}
根据programming guide，FERMI上grid支持3个维度，每个维度最大为65535，即blockIdx.x，.y,.z都是可以0~65535-1。
神马？还不够？
您可以进一步拆分成两个kernel来实现，或者买kepler，kepler在x维度上最大支持到（2^31）-1；y，z维度仍然是65535。您可以酌情修改算法实现。
祝您编码愉快~
:):):):):)
:3_48::3_48::3_48::3_48::3_48::3_48: