gpt4 book ai didi

matrix - 处理 CUDA 中的矩阵 : understanding basic concepts

转载 作者:行者123 更新时间:2023-12-04 21:07:37 25 4
gpt4 key购买 nike

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

我的问题是 我还没有找到对角线求和的方法 矩阵 cudaMalloc 'ed。我一直在尝试使用声明 threadIdx.x == blockIdx.x作为对角线的条件,但它的计算结果为 true只为他们俩在 0 .

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

template <typename T>
__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
  • 总和 dutemp_sx 的对角线, 减去 du来自 temp_dx 的对角线
  • 计算 temp_sx 的每个条目的平方函数和 temp_dx
  • 减去它们并将每个条目除以 2*du

  • 这个过程可以总结为 (f(un + du*e_i) - f(un - du*e_i))/2*du .

    我的问题是求和 du到矩阵的对角线 temp_sxtemp_dx就像我在 (*) 中尝试过的一样.我怎样才能做到这一点?

    编辑:现在调用一维块和线程;事实上, .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,这是每个块的最大线程数)。

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

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

    最佳答案

    根据对答案的评论中的建议,这是我设法解决问题的方法。该示例是可编译的,前提是您输入 helper_cuda.hhelper_string.h在同一目录中或添加 -I CUDA 示例的指令包含路径,与 CUDA 工具包一起安装。相关的变化只在内核中; main() 有一个小的变化虽然,因为我调用了双倍的资源来执行内核,但是 .y线程块网格的轴甚至根本没有使用,所以它没有产生任何错误。

    #include <stdio.h> 
    #include <stdlib.h>
    #include <iostream>
    #include <math.h>
    #include <assert.h>

    #include <cuda.h>
    #include <cuda_runtime.h>
    #include "helper_cuda.h"
    #include "helper_string.h"
    #include <fstream>

    #ifndef MAX
    #define MAX(a,b) ((a > b) ? a : b)
    #endif
    #define REAL sizeof(float)
    #define N 128
    #define BLOCK_SIZE 128
    #define NUM_BLOCKS ((N*N + BLOCK_SIZE - 1)/ BLOCK_SIZE)

    template <typename T>
    inline void printmatrix( T mat, int rows, int cols);
    template <typename T>
    __global__ void jacobian_kernel ( const T * A, T * J, const T t0, const T tn, const T h, const T * u0, const T * un, const T * un_old);
    template<typename T>
    __device__ void d_func(const T t, const T u[], T res[], const T h = 1);
    template<typename T>

    int main ()
    {
    float t0 = 0.; //float tn = 0.;
    float h = 0.1;

    float* u0 = (float*)malloc(REAL*N); for(int i = 0; i < N; ++i){u0[i] = i+1;}
    float* un = (float*)malloc(REAL*N); memcpy(un, u0, REAL*N);
    float* un_old = (float*)malloc(REAL*N); memcpy(un_old, u0, REAL*N);
    float* J = (float*)malloc(REAL*N*N);
    float* A = (float*)malloc(REAL*N*N); host_heat_matrix(A);

    float *d_u0;
    float *d_un;
    float *d_un_old;
    float *d_J;
    float *d_A;

    checkCudaErrors(cudaMalloc((void**)&d_u0, REAL*N)); //printf("1: %p\n", d_u0);
    checkCudaErrors(cudaMalloc((void**)&d_un, REAL*N)); //printf("2: %p\n", d_un);
    checkCudaErrors(cudaMalloc((void**)&d_un_old, REAL*N)); //printf("3: %p\n", d_un_old);
    checkCudaErrors(cudaMalloc((void**)&d_J, REAL*N*N)); //printf("4: %p\n", d_J);
    checkCudaErrors(cudaMalloc((void**)&d_A, REAL*N*N)); //printf("4: %p\n", d_J);
    checkCudaErrors(cudaMemcpy(d_u0, u0, REAL*N, cudaMemcpyHostToDevice)); assert(d_u0 != NULL);
    checkCudaErrors(cudaMemcpy(d_un, un, REAL*N, cudaMemcpyHostToDevice)); assert(d_un != NULL);
    checkCudaErrors(cudaMemcpy(d_un_old, un_old, REAL*N, cudaMemcpyHostToDevice)); assert(d_un_old != NULL);
    checkCudaErrors(cudaMemcpy(d_J, J, REAL*N*N, cudaMemcpyHostToDevice)); assert(d_J != NULL);
    checkCudaErrors(cudaMemcpy(d_A, A, REAL*N*N, cudaMemcpyHostToDevice)); assert(d_A != NULL);

    dim3 dimGrid(NUM_BLOCKS); std::cout << "NUM_BLOCKS \t" << dimGrid.x << "\n";
    dim3 dimBlock(BLOCK_SIZE); std::cout << "BLOCK_SIZE \t" << dimBlock.x << "\n";
    size_t shm_size = N*REAL; //std::cout << shm_size << "\n";

    //HERE IS A RELEVANT CHANGE OF THE MAIN, SINCE I WAS CALLING
    //THE KERNEL WITH A 2D GRID BUT WITHOUT USING THE .y AXIS,
    //WHILE NOW THE GRID IS 1D
    jacobian_kernel <<< dimGrid, dimBlock, shm_size >>> (d_A, d_J, t0, t0, h, d_u0, d_un, d_un_old);

    checkCudaErrors(cudaMemcpy(J, d_J, REAL*N*N, cudaMemcpyDeviceToHost)); //printf("4: %p\n", d_J);

    printmatrix( J, N, N);

    checkCudaErrors(cudaDeviceReset());
    free(u0);
    free(un);
    free(un_old);
    free(J);

    }

    template <typename T>
    __global__ void jacobian_kernel (
    const T * A,
    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;
    T* temp_du = &sm_temp_du;

    //HERE IS A RELEVANT CHANGE (*)
    if ( t < BLOCK_SIZE && b < NUM_BLOCKS )
    {
    temp_sx[b][t] = un[t]; //printf("temp_sx[%d] = %f\n", t,(temp_sx[b][t]));
    temp_dx[b][t] = un[t];
    //printf("t = %d, b = %d, t + b * blockDim.x = %d \n",t, b, tid);

    //HERE IS A NOTE (**)
    if ( t == b )
    {
    //printf("t = %d, b = %d \n",t, b);
    if ( tn == t0 )
    {
    *temp_du = u0[t]*0.001;

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

    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 = MAX( un[t] - un_old[t], 10e-6 );
    temp_sx[b][t] += *temp_du;
    temp_dx[b][t] -= *temp_du;
    }
    ;
    }
    //printf("du[%d] %f\n", tid, (*temp_du));
    __syncthreads();
    //printf("temp_sx[%d][%d] = %f\n", b, t, temp_sx[b][t]);
    //printf("temp_dx[%d][%d] = %f\n", b, t, temp_dx[b][t]);

    //d_func(tn, (temp_sx[b]), (temp_sx[b]), 1.f);
    //d_func(tn, (temp_dx[b]), (temp_dx[b]), 1.f);
    matvec_dev( tn, A, (temp_sx[b]), (temp_sx[b]), N, N, 1.f );
    matvec_dev( tn, A, (temp_dx[b]), (temp_dx[b]), N, N, 1.f );
    __syncthreads();
    //printf("temp_sx_later[%d][%d] = %f\n", b, t, (temp_sx[b][t]));
    //printf("temp_sx_later[%d][%d] - temp_dx_later[%d][%d] = %f\n", b,t,b,t, (temp_sx[b][t] - temp_dx[b][t]) / 2 * *temp_du);
    //if (t == b ) printf( "2du[%d]^-1 = %f\n",t, powf((2 * *temp_du), -1));
    J[tid] = (temp_sx[b][t] - temp_dx[b][t]) / (2 * *temp_du);
    }
    }

    template<typename T>
    __device__ void d_func(const T t, const T u[], T res[], const T h )
    {
    __shared__ float temp_u;
    temp_u = u[threadIdx.x];
    res[threadIdx.x] = h*powf( (temp_u), 2);
    }

    template <typename T>
    inline void printmatrix( T mat, int rows, int cols)
    {
    std::ofstream matrix_out;
    matrix_out.open( "heat_matrix.txt", std::ofstream::out);
    for( int i = 0; i < rows; i++)
    {
    for( int j = 0; j <cols; j++)
    {
    double next = mat[i + N*j];
    matrix_out << ( (next >= 0) ? " " : "") << next << " ";
    }
    matrix_out << "\n";
    }
    }

    相关变更在 (*) .在我使用之前 if (tid < N)这有两个缺点:
  • 首先,它是错误的,因为它应该是tid < N*N ,因为我的数据是二维的,而 tid是一个跟踪所有数据的全局索引。
  • 即使我写了 tid < N*N ,因为我将函数调用分成块,t < BLOCK_SIZE && b < NUM_BLOCKS我似乎更清楚索引在代码中的排列方式。

  • 此外,声明 t == b(**) 实际上是对对角元素进行操作的正确方法 的矩阵。它被评估的事实 true仅在 0是因为我上面的错误。

    感谢您的建议!

    关于matrix - 处理 CUDA 中的矩阵 : understanding basic concepts,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/41764331/

    25 4 0
    Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
    广告合作:1813099741@qq.com 6ren.com