a[rowIndex] = a_data + rowIndex * pitch / sizeof(float);
  然后在GPU设备端的内核函数中就能使用 
    a[i][j]的二维索引进行操作了,最后涉及到数据的拷贝则使用cudaMemcpy2D函数,函数中需要指定主机端和设备端的数据指针的每行内存大小,最后完整代码示例如下.
 
#include<iostream>
#include"book.h"
#include<cuda_runtime.h>
#include<device_launch_parameters.h>
#include"lock.h"
__global__ void testMat( float **a, int rows, int cols ){
    for(int tid_y = threadIdx.y + blockDim.y * blockIdx.y;tid_y < rows;
                                      tid_y += blockDim.y * gridDim.y){
        for(int tid_x = threadIdx.x + blockDim.x * blockIdx.x;tid_x < cols;
                                          tid_x += blockDim.x * gridDim.x){
            a[tid_y][tid_x] += 0.1;
int main( void ){
    int rows = 8;//定义矩阵行数M
    int cols = 4;//定义矩阵行数N
    //创建分配主机端和设备端的内存空间
    float **host_a,**dev_a;
    host_a = (float**)malloc(rows*sizeof(float*));
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a,rows*sizeof(float*) ) );
    float *host_a_data,*dev_a_data;size_t pitch;
    host_a_data = (float*)malloc(rows*cols*sizeof(float));
    HANDLE_ERROR( cudaMallocPitch( &dev_a_data, &pitch, cols*sizeof(float), rows ) );
    HANDLE_ERROR( cudaMemset2D( dev_a_data, pitch, 0, cols*sizeof(float), rows ) );
    for(int r = 0;r < rows;r++){
        host_a[r] = dev_a_data + r*pitch / sizeof(float);
    //主机端填充矩阵元素
    for(int i = 0;i < rows*cols;i++){
        host_a_data[i] = i;
    //主机端拷贝矩阵到设备端
    HANDLE_ERROR( cudaMemcpy( dev_a, host_a, rows*sizeof(float*), cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy2D( dev_a_data, pitch, host_a_data, cols*sizeof(float), 
                                cols*sizeof(float), rows, cudaMemcpyHostToDevice) );
    //设置线程网格和线程块的排布
    dim3 threads_rect(2,2);
    dim3 blocks_rect(2,2);
    //记录起始时间
    cudaEvent_t start,stop;float elapsedTime;
    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );
    //执行设备端的内核函数
    testMat<<<blocks_rect,threads_rect>>>(dev_a, rows, cols);
    //设备端拷贝矩阵到主机端
    HANDLE_ERROR( cudaMemcpy2D( host_a_data, cols*sizeof(float), dev_a_data, pitch, 
                                cols*sizeof(float), rows, cudaMemcpyDeviceToHost) );
    //获取结束时间，并显示计时结果
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, start, stop ) );
    printf( "Time to generate:  %3.2f ms\n", elapsedTime )




    
;
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
    //打印主机端矩阵数据中的计算结果
    for(int i = 0;i < rows*cols;i++){
        printf("%4.4f ",host_a_data[i]);
        if(i%cols == cols-1)
            printf("\n");
    //释放GPU设备端内存
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_a_data ) );
    //释放CPU主机端内存
    free( host_a );
    free( host_a_data );
二维矩阵按行求和
 
  对二维矩阵进行按行求和,也就是每一行单独求出一个和,最后有M个求和结果并以数组形式返回结果.首先需要设置线程网格中线程块的排布方式,这里有两种方式,第一种是M行1列的线程块,第二种是M行L列的线程块,L可取16或其它数字.第一种方式中,缺点是每行求和时的线程数量最大为1024,但是优点是更简单不需要进行多个线程块之间的原子累加操作.第二种方式中,优点是处理每行的线程数量可以远超1024,但是缺点是更复杂一点,需要进行原子操作,而原子操作的效率是很低的.具体如何使用,可以根据实际的数据,尝试不同的方式以确定最佳设置.最后完整代码示例如下.
 
#include<iostream>
#include"book.h"
#include<cuda_runtime.h>
#include<device_launch_parameters.h>
#include"lock.h"
__global__ void testMat_1( float **a, int rows, int cols ){
    /*对所有矩阵元素进行加0.1*/
    for(int tid_y = threadIdx.y + blockDim.y * blockIdx.y;tid_y < rows;
                                      tid_y += blockDim.y * gridDim.y){
        for(int tid_x = threadIdx.x + blockDim.x * blockIdx.x;tid_x < cols;
                                          tid_x += blockDim.x * gridDim.x){
            a[tid_y][tid_x] += 0.1;
__global__ void testMat_2( float **a,float *result, int rows, int cols ){
    /*对矩阵的每行求和,每个block线程块处理一行*/
    __shared__ float cache[1024];
    int cacheIndex = threadIdx.x;
    float temp = 0;
    int tid_y = blockIdx.x;
    for(int tid_x = threadIdx.x;tid_x < cols;tid_x += blockDim.x){
        //a[tid_y][tid_x] += 0.1;
        temp += a[tid_y][tid_x];
    cache[cacheIndex] = temp;
    __syncthreads();
    int i = blockDim.x / 2;
    while(i != 0){
        if(cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    if(cacheIndex == 0)
        result[tid_y] = cache[0];
__global__ void testMat_3( Lock lock,float **a,float *result, int rows, int cols ){
    /*对矩阵的每行求和,block线程块组成M行32列的线程网格*/
    __shared__ float cache[1024];
    int cacheIndex = threadIdx.x;
    float temp = 0;
    int tid_y = blockIdx.y;
    for(int tid_x = threadIdx.x + blockDim.x * blockIdx.x;tid_x < cols;
                                      tid_x += blockDim.x * gridDim.x){
    //for(int tid_x = threadIdx.x;tid_x < cols;tid_x += blockDim.x){
        //a[tid_y][tid_x] += 0.1;
        temp += a[tid_y][tid_x];
    cache[cacheIndex] = temp;
    __syncthreads();
    int i = blockDim.x / 2;
    while(i != 0){
        if(cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    if(cacheIndex == 0){
        lock.lock();
        result[tid_y] += cache[0];
        lock.unlock();
int main( void ){
    int rows = 1024;//定义矩阵行数M
    int cols = 131072;//定义矩阵列数N
    //创建分配主机端和设备端的内存空间
    float **host_a,**dev_a;
    host_a = (float**)malloc(rows*sizeof(float*));
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a,rows*sizeof(float




    
*) ) );
    float *host_a_data,*dev_a_data,*host_result,*dev_result;size_t pitch;
    host_a_data = (float*)malloc(rows*cols*sizeof(float));
    host_result = (float*)malloc(rows*sizeof(float));
    HANDLE_ERROR( cudaMallocPitch( &dev_a_data, &pitch, cols*sizeof(float), rows ) );
    HANDLE_ERROR( cudaMemset2D( dev_a_data, pitch, 0, cols*sizeof(float), rows ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_result,rows*sizeof(float) ) );
    for(int r = 0;r < rows;r++){
        host_a[r] = dev_a_data + r*pitch / sizeof(float);
        //host_a[r] = (float*)((char*)dev_a_data + r*pitch);
    //主机端填充矩阵元素
    for(int i = 0;i < rows*cols;i++){
        host_a_data[i] = 0.01;
    //主机端拷贝矩阵到设备端
    HANDLE_ERROR( cudaMemcpy( dev_a, host_a, rows*sizeof(float*), cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy2D( dev_a_data, pitch, host_a_data, cols*sizeof(float), 
                                cols*sizeof(float), rows, cudaMemcpyHostToDevice) );
    //设置线程网格和线程块的排布
    dim3 threads_rect(2,2);
    dim3 blocks_rect(2,2);
    //记录起始时间
    cudaEvent_t start,stop;float elapsedTime;
    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );
    //执行设备端的内核函数
    //testMat_1<<<blocks_rect,threads_rect>>>(dev_a, rows, cols);
    testMat_2<<<rows,1024>>>(dev_a,dev_result, rows, cols);
    //dim3 blocks_rect_2(1,rows);
    //Lock lock;
    //testMat_3<<<blocks_rect_2,1024>>>(lock, dev_a, dev_result, rows, cols);
    //设备端拷贝矩阵到主机端
    //HANDLE_ERROR( cudaMemcpy2D( host_a_data, cols*sizeof(float), dev_a_data, pitch, 
    //                            cols*sizeof(float), rows, cudaMemcpyDeviceToHost) );
    HANDLE_ERROR( cudaMemcpy( host_result, dev_result, rows*sizeof(float), cudaMemcpyDeviceToHost ) );
    //获取结束时间，并显示计时结果
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, start, stop ) );
    printf( "Time to generate:  %3.2f ms\n", elapsedTime );
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
    //打印主机端矩阵数据中的计算结果
    for(int i = 0;i < 20;i++){
        printf("%4.4f ",host_a_data[i]);
        if(i%cols == cols-1)
            printf("\n");
    printf("\n");
    for(int i = 0;i < 32;i++){
        printf("%f ",host_result[i]);
        if(i%4 == 3)
            printf("\n");
    //释放GPU设备端内存
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_a_data ) );
    //释放CPU主机端内存
    free( host_a );
    free( host_a_data );
二维矩阵乘法
 
  设矩阵都为n阶方阵,则设定 
    n×n的线程块网格,每个线程块1024个线程执行矢量点积运算,并执行求和归约,可以设定共享内存以提升内存读取速度.代码示例有.
 
__global__ void multiply(float **dist, float **x, float **result, int *sumColBoolIndexList, int Num, int m, int p){
    /*计算过程为q(:,index_Q)'*data(t).x./sum_Q(index_Q)'*/
    /*其中涉及到两个矩阵相乘A*B,设A为Num行m列,B为m行p列*/
    /*其中的Num行,等于非零列总数*/
    /*block线程块网格为x方向Num个,y方向p个*/
    /*每个线程块有1024个线程,负责A一行和B一列m个元素的矢量点积*/
    int tid_x = sumColBoolIndexList[blockIdx.x];
    int tid_y = blockIdx.y;
    __shared__ float cache[1024];
    int cacheIndex = threadIdx.x;
    float temp = 0;
    for(int tid_z = threadIdx.x;tid_z < m;tid_z += blockDim.x){
        /*负责长度为m的矢量点积*/
        temp += dist[tid_z][tid_x] * x[tid_y][tid_z];
    cache[cacheIndex] = temp;
    __syncthreads();
    int i = 0;
    while(i != 0){
        if(cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    if(cacheIndex == 0)
        result[tid_y][tid_x] = cache[0];
                                    通过矩阵按行求和与按列求和两个示例介绍CUDA并行算法设计的思路，希望对大家有所帮助。很多公司CUDA工程师面试也会考察这个题目。
1.        矩阵按行求和
矩阵A[M][N]，B[M];
B[i] =A[0]+A[1]+...+A[N-1];
串行代码：
  #define A(i)(j)  A[(i)*N+(j)]
void sum_row(T *A,T*B,int M,i...
                                      为更好地理解线程束执行的本质，将使用不同的执行配置分析下述的sumMatrixOnGPU2D核函数。使用nvprof配置指标，可以有助于理解为什么有些网格/块的维数组合比其他的组合更好。这些练习会提供网格和块的启发式算法，这是CUDA编程人员必备的技能。
二维矩阵求和的核函数如下所示：
__global__ void sumMatrixOnGPU2D(float *A, float *B, float *C, int NX, int NY)
   unsigned int ix = blockId
#include"iostream"
using namespace std;
/**********************************************************
本文通过一个矩阵加法的例子来说明如何使用网格和块来组织线程。
使用块和线程建立矩阵索引
通常情况下，一个矩阵用行优先的方法在全局内存中进行线性存储。如下图所示，这是一个8*6的矩阵。
在一个矩阵加法和核函数中，一个线程通常被分配一个数据元素来处理。首先要使用块和线程索引从全局内存中访问指定的数据。 接下来学习需要管理3种索引：
线程和块索引；
矩阵中给定点的坐标；
全局线性内存中的偏移量...
                                    CUDA矩阵元素求和在使用cuda 做图像高性能处理的时候，会很频繁的使用到矩阵的运算，于是打算做成一个相对完整的矩阵运算库，一旦使用调用即可，便再不用临时去写矩阵运算的规则和方法了。矩阵的加  减  乘  转置自然简单，在cuda中，矩阵的每个元素都对应着一个线程，于是矩阵加法转换为线程中的加法，减法变成了线程，乘法，转置也可以转化成类似的思路。当我写到求一个矩阵所有元素的和的时候，顺着思维惯性写
CUDA编程：GPU float 与 double 精度问题_longteng_guo的博客-CSDN博客
CUDA使用二维矩阵：能转化为一维矩阵尽量转化为一维矩阵，容易编程
在CUDA如何使用二维数组（**[M][N]）_Expressing Youself Using Code-CSDN博客
CUDA拷贝二维数组到GPU内存中_Dezeming的博客-CSDN博客
关于cuda_Malloc的第一个参数为什么是两个星星：
Cuda中的cuda_Malloc函数_Dez
                                    　　求矩阵每行的和？
　　可以把每行放入一个不同线程块，这样行与行之间进行粗粒度的并行。而对于每行，其对应的线程块中分配n个线程（对应行宽），使用共享存储器，让每个线程从显存中读取一个数至shared memory中，然后使用规约算法计算和。
代码如下：
#include "cuda_runtime.h" //CUDA运行时API
#include "device...
win10，PyCUDA: (2021, 1)，Python: 3.9.7
说明：调用时传入的grid、block参数，不能超过GPU的限制，所以需要先调用pycuda.driver库进行相关操作获取，网上有很多资料。
代码如下：
# -*- coding: utf-8 -*-
Created on Wed Feb 19 21:33:17