处理CUDA中的矩阵:理解基本概念

我正在构建一个CUDA内核来计算函数的数字N*Njacobian,使用有限差分; 在我提供的示例中,它是平方函数(向量的每个条目是平方的).主机编码在线性内存中分配,而我在内核中使用二维索引.

我的问题是,我还没有找到一种方法来对矩阵cudaMalloc编辑的对角线求和.我的尝试一直是使用该语句threadIdx.x == blockIdx.x作为对角线的条件,但它只评估true它们的两个0.

这是内核和编辑:我发布了整个代码作为答案,基于评论中的建议(main()基本相同,而内核不是)

template 
__global__ void jacobian_kernel (
                T * J,
                const T t0, 
                const T tn,
                const T h,
                const T * u0, 
                const T * un, 
                const T * un_old)
{
    T cgamma = 2 - sqrtf(2);
    const unsigned int t = threadIdx.x;
    const unsigned int b = blockIdx.x;
    const unsigned int tid = t + b * blockDim.x;
    /*__shared__*/ T temp_sx[BLOCK_SIZE][BLOCK_SIZE];
    /*__shared__*/ T temp_dx[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ T sm_temp_du[BLOCK_SIZE];

    T* temp_du = &sm_temp_du[0];

    if (tid < N )
    {
        temp_sx[b][t] = un[t]; 
        temp_dx[b][t] = un[t];

        if ( t == b )
        {
            if ( tn == t0 )
            {   
                temp_du[t] = u0[t]*0.001; 

                temp_sx[b][t] += temp_du[t]; //(*)
                temp_dx[b][t] -= temp_du[t];

                temp_sx[b][t] += ( abs( temp_sx[b][t] ) < 10e-6 ? 0.1 : 0 );
                temp_dx[b][t] += ( abs( temp_dx[b][t] ) < 10e-6 ? 0.1 : 0 );

                temp_sx[b][t] = ( temp_sx[b][t] == 0 ? 0.1 : temp_sx[b][t] );
                temp_dx[b][t] = ( temp_dx[b][t] == 0 ? 0.1 : temp_dx[b][t] );

            }

            else
            {
                temp_du[t] = MAX( un[t] - un_old[t], 10e-6 );
                temp_sx[b][t] += temp_du[t];
                temp_dx[b][t] -= temp_du[t];
            }
        }
        __syncthreads();

        //J = f(tn, un + du)
        d_func(tn, (temp_sx[b]), (temp_sx[b]), 1.f);
        d_func(tn, (temp_dx[b]), (temp_dx[b]), 1.f);

        __syncthreads();
        J[tid] = (temp_sx[b][t] - temp_dx[b][t]) * powf((2 * temp_du[t]), -1);

        //J[tid]*= - h*cgamma/2;
        //J[tid]+= ( t == b ? 1 : 0);
        //J[tid] = temp_J[tid];
    }
}   

计算雅可比的一般程序是

复制un到的每一行temp_sx和temp_dx

du作为一个0.01量级计算u0

综上所述du对角线的temp_sx,减去du从对角线temp_dx

计算temp_sx和的每个条目的平方函数temp_dx

减去它们并将每个条目除以 2*du

此过程可以概括为(f(un + du*e_i) - f(un - du*e_i))/2*du.

我的问题是要总结du的对角矩阵的中temp_sx和temp_dx像我在尝试(*).我怎样才能做到这一点？

编辑:现在调用1D块和线程; 事实上,.y在内核中根本没有使用轴.我用固定数量的共享内存调用内核

请注意,在int main()我调用内核时

#define REAL sizeof(float)
#define N 32
#define BLOCK_SIZE 16
#define NUM_BLOCKS ((N*N + BLOCK_SIZE - 1)/ BLOCK_SIZE)
...
dim3 dimGrid(NUM_BLOCKS,); 
dim3 dimBlock(BLOCK_SIZE); 
size_t shm_size = N*N*REAL;
jacobian_kernel <<< dimGrid, dimBlock, size_t shm_size >>> (...);

因此,我尝试处理块分割函数调用.在内核中对我使用的对角线求和if(threadIdx.x == blockIdx.x){...}.为什么这不正确？我问它是因为在调试和使代码打印语句时,它只评估true它们是否都是0.因此du[0]是唯一的数值,矩阵变为nan.请注意,这种方法适用于我构建的第一个代码,而我用内核调用了内核

jacobian_kernel <<< N, N >>> (...)

这样当threadIdx.x == blockIdx.x元素在对角线上时.这种方法不再适用,因为现在我需要处理更大的N(可能大于1024,这是每个块的最大线程数).

即使矩阵被分成块和线程,我应该在那里使用什么语句呢？

如果我应该分享一些其他信息,请告诉我.