斑竹您好：
我现在在做7575的矩阵的svd分解，采用的是单边雅克比的方法，您即使不了解也没关系，下面我会把必要的信息说明的；
其基本思想就是根据任意两列的各自的內积及相互的差积求出对应该两列的sin，cos值，对两列中相应元素做处理，使任意两列正交；
在这里d_a保存了7575的矩阵，d_v为单位的75*75的矩阵，根据d_cov中列求sin，cos，同时对d_cov和d_v相对应的列做相同的处理；
d_parity中保存了每次处理的列序号，
初始 （0,1）（2,3）（4,5）（6,7）。。。（72,73），74 flag=0，从0列开始取， 循环第一次时括号里的列队做处理，剩余的74列不处理，处理完后，对应的列号交换，
第二次 1（0,3）（2,5）（4,7）。。。。。。 （72,74） flag=1，从1列开始取， 第二次从第1个列号开始，每两列在一起处理，第一列不进行处理；
第三次（1，3）（0,5），（2,7）。。。。（ ，74） 72 flag=0，从0列开始取， 第三次在第二次的基础上交换处理的列号，剩余72列不处理；
。。。
这样经过75次后，每一列与其他列都进行了一次正交处理，希望我说清楚了，
我的程序再执行第前几次的时候感觉记过是对的，但执行次数多了后就不对了；
不知道是同步问题还是怎么了，希望斑竹能帮我看一下；
__global__ void svd_new(int *d_parity,float *d_a,float *d_v,int flag)
{
int y_id=blockIdx.y*blockDim.y+threadIdx.y;
int x_id=blockIdx.x*blockDim.x+threadIdx.x;
int threadIdx_y=threadIdx.y;
int threadIdx_x=threadIdx.x;
__shared__ int parity[75];
__shared__ float a[75*2];
__shared__ float aa[75*2];
__shared__ float a12[75];
__shared__ float v[75*2];
if (threadIdx_y<75&&threadIdx_x<1)
{
parity[threadIdx_y]=d_parity[threadIdx_y];
}
__syncthreads();
int col_num=x_id+flag;
a[threadIdx_x+threadIdx_y*2]=d_a[parity[col_num]+threadIdx_y*75];
aa[threadIdx_x+threadIdx_y*2]=a[threadIdx_x+threadIdx_y*2]*a[threadIdx_x+threadIdx_y*2];
v[threadIdx_x+threadIdx_y*2]=d_v[parity[col_num]+threadIdx_y*75];
__syncthreads();
if (threadIdx_y<75&&threadIdx_x<1)
{
a12[threadIdx_y]=a[threadIdx_y*2+threadIdx_x]*a[threadIdx_y*2+threadIdx_x+1];
}
__syncthreads();
int mid1=75/2;
int mid2=75%2;
int mid3;
while (mid1!=0&&threadIdx_y<mid1)
{
aa[threadIdx_y*2+threadIdx_x]=aa[threadIdx_y*2+threadIdx_x]+aa[(threadIdx_y+mid1+mid2)*2+threadIdx_x];
mid3=mid2+mid1;
mid1=mid3/2;
mid2=mid3%2;
__syncthreads();
}
mid1=75/2;
mid2=75%2;
while (mid1!=0&&threadIdx_y<mid1&&threadIdx_x<1)
{
a12[threadIdx_y]=a12[threadIdx_y]+a12[threadIdx_y+mid1+mid2];
mid3=mid2+mid1;
mid1=mid3/2;
mid2=mid3%2;
__syncthreads();
}
float dis1=aa[0];
float dis2=aa[1];
float dis12=a12[0];
__syncthreads();
// printf ("hello cuda blockId=%d threadId=%d dis1=%f dis2=%f dis12=%f\n",blockIdx.x, threadIdx.x, dis1,dis2,dis12);
float tao1,tan1,cos1,sin1;
float num1,num2;
int i=threadIdx_x;
tao1 = (dis1 - dis2) / (2 * dis12);
tan1 = mysign(tao1) / (fabs(tao1) + sqrt(1 + pow(tao1, 2)));
cos1 = 1 / sqrt(1 + pow(tan1, 2));
sin1 = cos1 * tan1;
int ii=int(pow(-1.0,i));
__syncthreads();
num1=a[threadIdx_y*2+i]*cos1+ii*a[threadIdx_y*2+i+ii]*sin1;
num2=v[threadIdx_y*2+i]*cos1+ii*v[threadIdx_y*2+i+ii]*sin1;
d_a[parity[col_num]+threadIdx_y*75]=num1;
d_v[parity[col_num]+threadIdx_y*75]=num2;
__syncthreads();
if (threadIdx_y<75-1&&threadIdx_x<1)
{
int col_id=threadIdx_y+flag;
d_parity[col_id]=parity[col_id+int(pow(-1.0,col_id+flag))];
}
__syncthreads();
}
有点长，希望斑竹能给点帮助，先谢谢了