正在求学CUDA 刚入门 写了一个矩阵求逆的次序,给矩阵优化时出错了 求大神指导,上面代码是应用分享内部存款和储

作者:操作系统    发布时间:2020-01-26 19:00     浏览次数 :

[返回]

#defineblock_size512#defineN1024___global__voidkernel1(floatA[][N],floatE[][N],intn){dim3block(1,N);dim3threads(N);inttx=threadIdx.x;intty=threadIdx.y;intbx=blockIdx.x;intby=blockIdx.y;for(intn=0;nN;n++){floattemp=A[n][n],c,d;for(inta=by*N,b=by*N;a=by*N+N-1;a+=32*32,b+=32*32){__shared__floatAs[32][32];__shared__floatBs[32][32];As[ty][tx]=A[n][a+tx];Bs[ty][tx]=E[n][b+tx];c=As[ty][tx]/temp;d=Bs[ty][tx]/temp;__syncthreads();}A[n][by*N+tx]=c;E[n][by*N+tx]=d;}}intsize=sizeof(float)*(N*N);float*h_A=(float*)malloc(size);float*h_E=(float*)malloc(size);float*h_F=(float*)malloc(size);float*d_A,*d_E;for(inti=0;iN;i++){for(intj=0;jN;j++){h_A[i*N+j]=A[i][j];h_E[i*N+j]=E[i][j];}}cudaMalloc((void**)d_A,size);cudaMalloc((void**)d_E,size);cudaMemcpy(d_A,h_A,size,cudaMemcpyHostToDevice);cudaMemcpy(d_E,h_E,size,cudaMemcpyHostToDevice);cudaEvent_tstart,stop;for(ints=0;s10;s++){cudaEventCreate(start);cudaEventCreate(stop);cudaEventRecord(start,0);kernel1N,N((float(*)[N])d_A,(float(*)[N])d_E,i);cudaMemcpy(h_F,d_E,size,cudaMemcpyDeviceToHost);cudaMemcpy(h_E,d_E,size,cudaMemcpyDeviceToHost);cudaEventRecord(stop,0);cudaEventSynchronize(stop);cudaEventElapsedTime(elapsedTime,start,stop);std::cout"GPU消耗时间:"elapsedTimestd::endl;printf("n");}for(inti=0;iN;i++){for(intj=0;jN;j++){printf("%.3f",h_F[i*N+j]);}printf("n");}free(h_A);free(h_E);cudaFree(d_A);cudaFree(d_E);cudaEventDestroy(start);cudaEventDestroy(stop);return0;}

图片 1

CUDA 是 NVIDIA 的 GPGPU 模型,它使用 C 语言为基础,可以直接以大多数人熟悉的 C 语言,写出在显示芯片上执行的程序,而不需要去学习特定的显示芯片的指令或是特殊的结构。”

各位大神,运行这个程序会在srand(time(NULL));处显示“标识符在其出现的地方是未被定义”的错误。随后我在头文件处加了#includetime.h运行后出现了如图所示的文件,以及运行结果,但是并不是想得到的结果,请问各位大神如何解决这个问题呢?谢谢啦#includeiostream#includecmath#includefstream#includestring#includesstream#includecstdlib#includetime.husingnamespacestd;constintna=-9999;//AllNA'sinthedataaresavedas-9999.Ourclaculationexcludesthesedatapoints.constintrange=151;//No.ofsnapshotsavailableconstintN=512;//Dimentionofeachsnapshot.E.g.amatrixof512x512cells.//-------------------------initializethematrix-----------------------////Thisfunctioncreatesonerandommatrixwithdensityxvoidcreate_random_matrix(int*A,floatx){inti,j;for(i=0;iN;i++){for(j=0;jN;j++){floatnumber=rand()/(float)RAND_MAX;if(numberx)A[i*N+j]=1;elseA[i*N+j]=0;}}}//-----------------------meancover-----------------------------------////Calculatesmeandensityofthematrixofsize"size"doubledensity(double*A,intsize){inti,j;doublecount=0,sum=0;for(i=0;isize;i++){for(j=0;jsize;j++){if(A[i*size+j]!=na){sum+=A[i*size+j];count+=1;}}}doubleaverage=sum/count;//cout"mean="averageendl;returnaverage;}//-----------------------coarse-graining-------------------------------////Thisfunctioncoarse-grainsthematrixwithcoarse-graininglength"cgLength"andreturnsthereducedmatrixofdimentionreducedbycgLengthvoidcoarse_grain(int*A,double*A_cg,intn,intcgLength){//nisthesizeofcoarse-grainedmatrixinti,j,k,l;for(i=0;in;i++){for(j=0;jn;j++){intinit_row=i*cgLength;intend_row=(i+1)*cgLength-1;intinit_col=j*cgLength;intend_col=(j+1)*cgLength-1;doublesum=0,count=0;for(k=init_row;kend_row+1;k++){for(l=init_col;lend_col+1;l++){if(A[k*N+l]==0||A[k*N+l]==1){sum+=A[k*N+l];count+=1;}}}doublereduced_mean=sum/count;A_cg[i*n+j]=reduced_mean;}}}//---------------------------------variance--------------------------------------//doublespVar(double*A,intsize,doublemean){doublesum=0,count=0;for(inti=0;isize;i++){for(intj=0;jsize;j++){if(A[i*size+j]!=na){sum=sum+A[i*size+j]*A[i*size+j];count+=1;}}}doublevar=sum/count-mean*mean;//cout"var="varendl;returnvar;}//--------------------------------skewness----------------------------------------------//floatcg_skew(double*A,intn,floatmean,floatvariance){//nisthedimensionofcoarsegrainingi.e.dimensionofthesubmatrixinti,j,k,l,count=0;floatskewness;floatsum_cube=0;for(i=0;in*n;i++){sum_cube+=A[i]*A[i]*A[i];}floatmean_cube=sum_cube/(float)(n*n);if(variance==0)skewness=0;elseskewness=(mean_cube-3.0*mean*variance-mean*mean*mean)/pow(variance,1.5);//coutmean"t"variance"t"skewnessendl;returnskewness;}//-------------------------------Correlationatlag1---------------------------------//doublespCor(double*A,intsize,doublemean,doublevariance){doublesum=0;inti,j;intdown,right;for(i=0;isize;i++){for(j=0;jsize;j++){down=(i+1)%size;right=(j+1)%size;sum=sum+(A[i*size+j]*A[down*size+j])+(A[i*size+j]*A[i*size+right]);}}doublecorLag1;if(variance==0)corLag1=0;elsecorLag1=(sum/(double)(2*size*size)-mean*mean)/variance;//cout"meansq="mean*mean"t""cor="corLag1endl;returncorLag1;}//--------------------------------------inttostringconverter----------------------------//stringIntToStr(intn){stringstreamresult;resultn;returnresult.str();}//----------------------------------mainfunction-----------------------------//intmain(){srand(time(NULL));floatinitial_density;doublemean,correlation,variance;int*mat=newint[N*N];int*full=newint[N*N*range];//~~~~~readingdatafromafile~~~~~~~////Thispartreadsthedatafromafileandsavesitinthearray"full"ifstreammyfile;myfile.open("saved_snapshots.dat");if(myfile.is_open()){for(intx=0;xrange;x++){for(inti=0;iN;i++){for(intj=0;jN;j++){myfilefull[x*N*N+i*N+j];}}}}elsecout"Unabletoopenfile";myfile.close();//~~~~~~~~~~~~~~~~~~////Loopoverdifferentcoarse-graininglengthsfor(intcgLength=33;cgLength65;cgLength++){//~~~declarereducedmatrixorcoarse-grainedmatrix~~~~~~~~~~//intrem;if(N%cgLength==0)rem=0;elserem=N%cgLength;intL=N-rem;intn=L/cgLength;double*mat_cg=newdouble[n*n];//~~~~~~~~~~~~~~~////Followinglooprunsoverallthesnapshotsfor(intx=0;xrange;x++){for(inti=0;iN;i++){for(intj=0;jN;j++){mat[i*N+j]=full[x*N*N+i*N+j];}}coarse_grain(mat,mat_cg,n,cgLength);mean=density(mat_cg,n);variance=spVar(mat_cg,n,mean);correlation=spCor(mat_cg,n,mean,variance);//~~~~~savethecorrelationandvariance~~~~~~~~~////Thispartcreatesdifferentfilesforsavingvarianceandcorrelationthesnapshotsatdifferentcoarse-graininglengthsstringfilename1="VarianceCg"+IntToStr(cgLength)+".dat";stringfilename2="CorrelationCg"+IntToStr(cgLength)+".dat";ofstreamoutdata1;ofstreamoutdata2;outdata1.open(filename1.c_str(),ios::app);outdata2.open(filename2.c_str(),ios::app);if(outdata1.is_open()){outdata1varianceendl;outdata2correlationendl;outdata1.close();outdata2.close();}elsecout"oudatanotopen";//~~~~~~~~~~~~~~~~~~~~~~//}delete[]mat_cg;}delete[]mat;delete[]full;return0;}

编者注:NVIDIA的GeFoce 8800GTX发布后,它的通用计算架构CUDA经过一年多的推广后,现在已经在有相当多的论文发表,在商业应用软件等方面也初步出现了视频编解码、金融、地质勘探、科学计算等领域的产品,是时候让我们对其作更深一步的了解。为了让大家更容易了解CUDA,我们征得Hotball的本人同意,发表他最近亲自撰写的本文。这篇文章的特点是深入浅出,也包含了hotball本人编写一些简单CUDA程序的亲身体验,对于希望了解CUDA的读者来说是非常不错的入门文章,PCINLIFE对本文的发表没有作任何的删减,主要是把一些台湾的词汇转换成大陆的词汇以及作了若干"编者注"的注释。

现代的显示芯片已经具有高度的可程序化能力,由于显示芯片通常具有相当高的内存带宽,以及大量的执行单元,因此开始有利用显示芯片来帮助进行一些计算工作的想法,即 GPGPU。CUDA 即是 NVIDIA 的 GPGPU 模型。

NVIDIA 的新一代显示芯片,包括 GeForce 8 系列及更新的显示芯片都支持 CUDA。NVIDIA 免费提供 CUDA 的开发工具(包括 Windows 版本和 Linux 版本)、程序范例、文件等等,可以在 CUDA Zone 下载。

GPGPU 的优缺点

使用显示芯片来进行运算工作,和使用 CPU 相比,主要有几个好处:

  1. 显示芯片通常具有更大的内存带宽。例如,NVIDIA 的 GeForce 8800GTX 具有超过 50GB/s 的内存带宽,而目前高阶 CPU 的内存带宽则在 10GB/s 左右。
  2. 显示芯片具有更大量的执行单元。例如 GeForce 8800GTX 具有 128 个 "stream processors",频率为 1.35GHz。CPU 频率通常较高,但是执行单元的数目则要少得多。
  3. 和高阶 CPU 相比,显卡的价格较为低廉。例如目前一张 GeForce 8800GT 包括 512MB 内存的价格,和一颗 2.4GHz 四核心 CPU 的价格相若。

当然,使用显示芯片也有它的一些缺点:

  1. 显示芯片的运算单元数量很多,因此对于不能高度并行化的工作,所能带来的帮助就不大。
  2. 显示芯片目前通常只支持 32 bits 浮点数,且多半不能完全支持 IEEE 754 规格, 有些运算的精确度可能较低。目前许多显示芯片并没有分开的整数运算单元,因此整数运算的效率较差。
  3. 显示芯片通常不具有分支预测等复杂的流程控制单元,因此对于具有高度分支的程序,效率会比较差。
  4. 目前 GPGPU 的程序模型仍不成熟,也还没有公认的标准。例如 NVIDIA 和 AMD/ATI 就有各自不同的程序模型。

整体来说,显示芯片的性质类似 stream processor,适合一次进行大量相同的工作。CPU 则比较有弹性,能同时进行变化较多的工作。

CUDA 架构

CUDA 是 NVIDIA 的 GPGPU 模型,它使用 C 语言为基础,可以直接以大多数人熟悉的 C 语言,写出在显示芯片上执行的程序,而不需要去学习特定的显示芯片的指令或是特殊的结构。

在 CUDA 的架构下,一个程序分为两个部份:host 端和 device 端。Host 端是指在 CPU 上执行的部份,而 device 端则是在显示芯片上执行的部份。Device 端的程序又称为 "kernel"。通常 host 端程序会将数据准备好后,复制到显卡的内存中,再由显示芯片执行 device 端程序,完成后再由 host 端程序将结果从显卡的内存中取回。

图片 2

由于 CPU 存取显卡内存时只能透过 PCI Express 接口,因此速度较慢(PCI Express x16 的理论带宽是双向各 4GB/s),因此不能太常进行这类动作,以免降低效率。

在 CUDA 架构下,显示芯片执行时的最小单位是 thread。数个 thread 可以组成一个 block。一个 block 中的 thread 能存取同一块共享的内存,而且可以快速进行同步的动作。

每一个 block 所能包含的 thread 数目是有限的。不过,执行相同程序的 block,可以组成 grid。不同 block 中的 thread 无法存取同一个共享的内存,因此无法直接互通或进行同步。因此,不同 block 中的 thread 能合作的程度是比较低的。不过,利用这个模式,可以让程序不用担心显示芯片实际上能同时执行的 thread 数目限制。例如,一个具有很少量执行单元的显示芯片,可能会把各个 block 中的 thread 顺序执行,而非同时执行。不同的 grid 则可以执行不同的程序(即 kernel)。

Grid、block 和 thread 的关系,如下图所示:

图片 3

每个 thread 都有自己的一份 register 和 local memory 的空间。同一个 block 中的每个 thread 则有共享的一份 share memory。此外,所有的 thread(包括不同 block 的 thread)都共享一份 global memory、constant memory、和 texture memory。不同的 grid 则有各自的 global memory、constant memory 和 texture memory。这些不同的内存的差别,会在之后讨论。

执行模式

由于显示芯片大量并行计算的特性,它处理一些问题的方式,和一般 CPU 是不同的。主要的特点包括:

  1. 内存存取 latency 的问题:CPU 通常使用 cache 来减少存取主内存的次数,以避免内存 latency 影响到执行效率。显示芯片则多半没有 cache(或很小),而利用并行化执行的方式来隐藏内存的 latency(即,当第一个 thread 需要等待内存读取结果时,则开始执行第二个 thread,依此类推)。
  2. 分支指令的问题:CPU 通常利用分支预测等方式来减少分支指令造成的 pipeline bubble。显示芯片则多半使用类似处理内存 latency 的方式。不过,通常显示芯片处理分支的效率会比较差。

因此,最适合利用 CUDA 处理的问题,是可以大量并行化的问题,才能有效隐藏内存的 latency,并有效利用显示芯片上的大量执行单元。使用 CUDA 时,同时有上千个 thread 在执行是很正常的。因此,如果不能大量并行化的问题,使用 CUDA 就没办法达到最好的效率了。

 

 

 

CUDA Toolkit的安装

目前 NVIDIA 提供的 CUDA Toolkit(可从这里下载)支持 Windows (32 bits 及 64 bits 版本)及许多不同的 linux 版本。

CUDA Toolkit 需要配合 C/C++ compiler。在 Windows 下,目前只支持 Visual Studio 7.x 及 Visual Studio 8(包括免费的 Visual Studio C++ 2005 Express)。Visual Studio 6 和 gcc 在 Windows 下是不支援的。在 Linux 下则只支援 gcc。

这里简单介绍一下在 Windows 下设定并使用 CUDA 的方式。

下载及安装

在 Windows 下,CUDA Toolkit 和 CUDA SDK 都是由安装程序的形式安装的。CUDA Toolkit 包括 CUDA 的基本工具,而 CUDA SDK 则包括许多范例程序以及链接库。基本上要写 CUDA 的程序,只需要安装 CUDA Toolkit 即可。不过 CUDA SDK 仍值得安装,因为里面的许多范例程序和链接库都相当有用。

CUDA Toolkit 安装完后,预设会安装在 C:/CUDA 目录里。其中包括几个目录:

  • bin -- 工具程序及动态链接库
  • doc -- 文件
  • include -- header 檔
  • lib -- 链接库档案
  • open64 -- 基于 Open64 的 CUDA compiler
  • src -- 一些原始码

安装程序也会设定一些环境变量,包括:

  • CUDA_BIN_PATH -- 工具程序的目录,默认为 C:/CUDA/bin
  • CUDA_INC_PATH -- header 文件的目录,默认为 C:/CUDA/inc
  • CUDA_LIB_PATH -- 链接库文件的目录,默认为 C:/CUDA/lib

在 Visual Studio 中使用 CUDA

CUDA 的主要工具是 nvcc,它会执行所需要的程序,将 CUDA 程序代码编译成执行档 (或 object 檔) 。在 Visual Studio 下,我们透过设定 custom build tool 的方式,让 Visual Studio 会自动执行 nvcc。

这里以 Visual Studio 2005 为例:

  1. 首先,建立一个 Win32 Console 模式的 project(在 Application Settings 中记得勾选 Empty project),并新增一个档案,例如 main.cu。
  2. 在 main.cu 上右键单击,并选择 Properties。点选 General,确定 Tool 的部份是选择 Custom Build Tool
  3. 选择 Custom Build Step,在 Command Line 使用以下设定:
    • Release 模式:"$(CUDA_BIN_PATH)/nvcc.exe" -ccbin "$(VCInstallDir)bin" -c -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)/$(InputName).obj $(InputFileName)
    • Debug 模式:"$(CUDA_BIN_PATH)/nvcc.exe" -ccbin "$(VCInstallDir)bin" -c -D_DEBUG -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)/$(InputName).obj $(InputFileName)
  4. 如果想要使用软件仿真的模式,可以新增两个额外的设定:
    • EmuRelease 模式:"$(CUDA_BIN_PATH)/nvcc.exe" -ccbin "$(VCInstallDir)bin" -deviceemu -c -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)/$(InputName).obj $(InputFileName)
    • EmuDebug 模式:"$(CUDA_BIN_PATH)/nvcc.exe" -ccbin "$(VCInstallDir)bin" -deviceemu -c -D_DEBUG -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)/$(InputName).obj $(InputFileName)
  5. 对所有的配置文件,在 Custom Build Step 的 Outputs 中加入 $(ConfigurationName)/$(InputName).obj。
  6. 选择 project,右键单击选择 Properties,再点选 Linker。对所有的配置文件修改以下设定:
    • General/Enable Incremental Linking:No
    • General/Additional Library Directories:$(CUDA_LIB_PATH)
    • Input/Additional Dependencies:cudart.lib

这样应该就可以直接在 Visual Studio 的 IDE 中,编辑 CUDA 程序后,直接 build 以及执行程序了。

 


CUDA和Visual C++2005 ide的设置比较复杂,OpenHero贡献了解决方案

CUDA VS2005 Wizard:

visual assist 支持cu文件:

语法高亮:

 

 

 

 

第一个CUDA程序

CUDA 目前有两种不同的 API:Runtime API 和 Driver API,两种 API 各有其适用的范围。由于 runtime API 较容易使用,一开始我们会以 runetime API 为主。

CUDA 的初始化

首先,先建立一个档案 first_cuda.cu。如果是使用 Visual Studio 的话,则请先按照这里的设定方式设定 project。

要使用 runtime API 的时候,需要 include cuda_runtime.h。所以,在程序的最前面,加上

#include <stdio.h>
#include <cuda_runtime.h>

接下来是一个 InitCUDA 函式,会呼叫 runtime API 中,有关初始化 CUDA 的功能:

bool InitCUDA()
{
    int count;

    cudaGetDeviceCount(&count);
    if(count == 0) {
        fprintf(stderr, "There is no device./n");
        return false;
    }

    int i;
    for(i = 0; i < count; i++) {
        cudaDeviceProp prop;
        if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
            if(prop.major >= 1) {
                break;
            }
        }
    }

    if(i == count) {
        fprintf(stderr, "There is no device supporting CUDA 1.x./n");
        return false;
    }

    cudaSetDevice(i);

    return true;
}

这个函式会先呼叫 cudaGetDeviceCount 函式,取得支持 CUDA 的装置的数目。如果系统上没有支持 CUDA 的装置,则它会传回 1,而 device 0 会是一个仿真的装置,但不支持 CUDA 1.0 以上的功能。所以,要确定系统上是否有支持 CUDA 的装置,需要对每个 device 呼叫 cudaGetDeviceProperties 函式,取得装置的各项数据,并判断装置支持的 CUDA 版本(prop.major 和 prop.minor 分别代表装置支持的版本号码,例如 1.0 则 prop.major 为 1 而 prop.minor 为 0)。

透过 cudaGetDeviceProperties 函式可以取得许多数据,除了装置支持的 CUDA 版本之外,还有装置的名称、内存的大小、最大的 thread 数目、执行单元的频率等等。详情可参考 NVIDIA 的 CUDA Programming Guide。

在找到支持 CUDA 1.0 以上的装置之后,就可以呼叫 cudaSetDevice 函式,把它设为目前要使用的装置。

最后是 main 函式。在 main 函式中我们直接呼叫刚才的 InitCUDA 函式,并显示适当的讯息:

int main()
{
    if(!InitCUDA()) {
        return 0;
    }

    printf("CUDA initialized./n");

    return 0;
}

这样就可以利用 nvcc 来 compile 这个程序了。使用 Visual Studio 的话,若按照先前的设定方式,可以直接 Build Project 并执行。

nvcc 是 CUDA 的 compile 工具,它会将 .cu 檔拆解出在 GPU 上执行的部份,及在 host 上执行的部份,并呼叫适当的程序进行 compile 动作。在 GPU 执行的部份会透过 NVIDIA 提供的 compiler 编译成中介码,而 host 执行的部份则会透过系统上的 C++ compiler 编译(在 Windows 上使用 Visual C++ 而在 Linux 上使用 gcc)。

编译后的程序,执行时如果系统上有支持 CUDA 的装置,应该会显示 CUDA initialized. 的讯息,否则会显示相关的错误讯息。

利用 CUDA 进行运算

到目前为止,我们的程序并没有做什么有用的工作。所以,现在我们加入一个简单的动作,就是把一大堆数字,计算出它的平方和。

首先,把程序最前面的 include 部份改成:

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

#define DATA_SIZE 1048576

int data[DATA_SIZE];

并加入一个新函式 GenerateNumbers:

void GenerateNumbers(int *number, int size)
{
    for(int i = 0; i < size; i++) {
        number[i] = rand() % 10;
    }
}

这个函式会产生一大堆 0 ~ 9 之间的随机数。

要利用 CUDA 进行计算之前,要先把数据复制到显卡内存中,才能让显示芯片使用。因此,需要取得一块适当大小的显卡内存,再把产生好的数据复制进去。在 main 函式中加入:

    GenerateNumbers(data, DATA_SIZE);

    int* gpudata, *result;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result, sizeof(int));
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);

上面这段程序会先呼叫 GenerateNumbers 产生随机数,并呼叫 cudaMalloc 取得一块显卡内存(result 则是用来存取计算结果,在稍后会用到),并透过 cudaMemcpy 将产生的随机数复制到显卡内存中。cudaMalloc 和 cudaMemcpy 的用法和一般的 malloc 及 memcpy 类似,不过 cudaMemcpy 则多出一个参数,指示复制内存的方向。在这里因为是从主内存复制到显卡内存,所以使用 cudaMemcpyHostToDevice。如果是从显卡内存到主内存,则使用 cudaMemcpyDeviceToHost。这在之后会用到。

接下来是要写在显示芯片上执行的程序。在 CUDA 中,在函式前面加上 __global__ 表示这个函式是要在显示芯片上执行的。因此,加入以下的函式:

__global__ static void sumOfSquares(int *num, int* result)
{
    int sum = 0;
    int i;
    for(i = 0; i < DATA_SIZE; i++) {
        sum += num[i] * num[i];
    }

    *result = sum;
}

在显示芯片上执行的程序有一些限制,例如它不能有传回值。其它的限制会在之后提到。

接下来是要让 CUDA 执行这个函式。在 CUDA 中,要执行一个函式,使用以下的语法:

    函式名称<<<block 数目, thread 数目, shared memory 大小>>>(参数...);

呼叫完后,还要把结果从显示芯片复制回主内存上。在 main 函式中加入以下的程序:

    sumOfSquares<<<1, 1, 0>>>(gpudata, result);

    int sum;
    cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);

    printf("sum: %d/n", sum);

因为这个程序只使用一个 thread,所以 block 数目、thread 数目都是 1。我们也没有使用到任何 shared memory,所以设为 0。编译后执行,应该可以看到执行的结果。

为了确定执行的结果正确,我们可以加上一段以 CPU 执行的程序代码,来验证结果:

    sum = 0;
    for(int i = 0; i < DATA_SIZE; i++) {
        sum += data[i] * data[i];
    }
    printf("sum (CPU): %d/n", sum);

编译后执行,确认两个结果相同。

计算运行时间

 

CUDA 提供了一个 clock 函式,可以取得目前的 timestamp,很适合用来判断一段程序执行所花费的时间(单位为 GPU 执行单元的频率)。这对程序的优化也相当有用。要在我们的程序中记录时间,把 sumOfSquares 函式改成:

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    int sum = 0;
    int i;
    clock_t start = clock();
    for(i = 0; i < DATA_SIZE; i++) {
        sum += num[i] * num[i];
    }

    *result = sum;
    *time = clock() - start;
}

把 main 函式中间部份改成:

    int* gpudata, *result;
    clock_t* time;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result, sizeof(int));
    cudaMalloc((void**) &time, sizeof(clock_t));
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);

    sumOfSquares<<<1, 1, 0>>>(gpudata, result, time);

    int sum;
    clock_t time_used;
    cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(&time_used, time, sizeof(clock_t),
        cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);

    printf("sum: %d time: %d/n", sum, time_used);

编译后执行,就可以看到执行所花费的时间了。

如果计算实际运行时间的话,可能会注意到它的执行效率并不好。这是因为我们的程序并没有利用到 CUDA 的主要的优势,即并行化执行。在下一段文章中,会讨论如何进行优化的动作。

 

在上一篇文章中,我们做了一个计算一大堆数字的平方和的程序。不过,我们也提到这个程序的执行效率并不理想。当然,实际上来说,如果只是要做计算平方和的动作,用 CPU 做会比用 GPU 快得多。这是因为平方和的计算并不需要太多运算能力,所以几乎都是被内存带宽所限制。因此,光是把数据复制到显卡内存上的这个动作,所需要的时间,可能已经和直接在 CPU 上进行计算差不多了。

不过,如果进行平方和的计算,只是一个更复杂的计算过程的一部份的话,那么当然在 GPU 上计算还是有它的好处的。而且,如果数据已经在显卡内存上(例如在 GPU 上透过某种算法产生),那么,使用 GPU 进行这样的运算,还是会比较快的。

刚才也提到了,由于这个计算的主要瓶颈是内存带宽,所以,理论上显卡的内存带宽是相当大的。这里我们就来看看,倒底我们的第一个程序,能利用到多少内存带宽。

 

程序的并行化

 

我们的第一个程序,并没有利用到任何并行化的功能。整个程序只有一个 thread。在 GeForce 8800GT 上面,在 GPU 上执行的部份(称为 "kernel")大约花费 640M 个频率。GeForce 8800GT 的执行单元的频率是 1.5GHz,因此这表示它花费了约 0.43 秒的时间。1M 个 32 bits 数字的数据量是 4MB,因此,这个程序实际上使用的内存带宽,只有 9.3MB/s 左右!这是非常糟糕的表现。

为什么会有这样差的表现呢?这是因为 GPU 的架构特性所造成的。在 CUDA 中,一般的数据复制到的显卡内存的部份,称为 global memory。这些内存是没有 cache 的,而且,存取 global memory 所需要的时间(即 latency)是非常长的,通常是数百个 cycles。由于我们的程序只有一个 thread,所以每次它读取 global memory 的内容,就要等到实际读取到数据、累加到 sum 之后,才能进行下一步。这就是为什么它的表现会这么的差。

由于 global memory 并没有 cache,所以要避开巨大的 latency 的方法,就是要利用大量的 threads。假设现在有大量的 threads 在同时执行,那么当一个 thread 读取内存,开始等待结果的时候,GPU 就可以立刻切换到下一个 thread,并读取下一个内存位置。因此,理想上当 thread 的数目够多的时候,就可以完全把 global memory 的巨大 latency 隐藏起来了。

要怎么把计算平方和的程序并行化呢?最简单的方法,似乎就是把数字分成若干组,把各组数字分别计算平方和后,最后再把每组的和加总起来就可以了。一开始,我们可以把最后加总的动作,由 CPU 来进行。

首先,在 first_cuda.cu 中,在 #define DATA_SIZE 的后面增加一个 #define,设定 thread 的数目:

#define DATA_SIZE    1048576
#define THREAD_NUM   256

接着,把 kernel 程序改成:

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    const int tid = threadIdx.x;
    const int size = DATA_SIZE / THREAD_NUM;
    int sum = 0;
    int i;
    clock_t start;
    if(tid == 0) start = clock();
    for(i = tid * size; i < (tid + 1) * size; i++) {
       sum += num[i] * num[i];
    }

    result[tid] = sum;
    if(tid == 0) *time = clock() - start;
}

程序里的 threadIdx 是 CUDA 的一个内建的变量,表示目前的 thread 是第几个 thread(由 0 开始计算)。以我们的例子来说,会有 256 个 threads,所以同时会有 256 个 sumOfSquares 函式在执行,但每一个的 threadIdx.x 则分别会是 0 ~ 255。利用这个变量,我们就可以让每一份函式执行时,对整个数据不同的部份计算平方和。另外,我们也让计算时间的动作,只在 thread 0(即 threadIdx.x = 0 的时候)进行。

同样的,由于会有 256 个计算结果,所以原来存放 result 的内存位置也要扩大。把 main 函式中的中间部份改成:

    int* gpudata, *result;
    clock_t* time;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
    cudaMalloc((void**) &time, sizeof(clock_t));
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);

    sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);

    int sum[THREAD_NUM];
    clock_t time_used;
    cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM, 
        cudaMemcpyDeviceToHost);
    cudaMemcpy(&time_used, time, sizeof(clock_t),
        cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);
    cudaFree(time);

可以注意到我们在呼叫 sumOfSquares 函式时,指定 THREAD_NUM 为 thread 的数目。最后,在 CPU 端把计算好的各组数据的平方和进行加总:

    int final_sum = 0;
    for(int i = 0; i < THREAD_NUM; i++) {
        final_sum += sum[i];
    }

    printf("sum: %d  time: %d/n", final_sum, time_used);

    final_sum = 0;
    for(int i = 0; i < DATA_SIZE; i++) {
        sum += data[i] * data[i];
    }
    printf("sum (CPU): %d/n", final_sum);

编译后执行,确认结果和原来相同。

这个版本的程序,在 GeForce 8800GT 上执行,只需要约 8.3M cycles,比前一版程序快了 77 倍!这就是透过大量 thread 来隐藏 latency 所带来的效果。

下一篇:没有了