由於CPU的摩爾定律已經不再適用，目前加速程式的最佳選擇就是通過GPU並行。經過幾天的摸索後，完成了這個利用GPU加速積分演算法的小實驗。

數值積分中最常用的方法之一就是辛普森積分法，首先我們寫出一段三階辛普森積分的小程式：

double Simpson_integ (int n_steps, double a, double b, double (*func)( double t))
{
    int n;
    double sum=0;
    double tk=0., dt;
    double f1, f2, f3;
   
    dt=(b-a)/n_steps;
    for(n=0;n<n_steps;n++)
    {
        tk=a+n*dt;
        f1=(*func)(tk );
        f2=(*func)(tk+dt/2);
        f3=(*func)(tk+dt);
        sum+=dt*(f1+4*f2+f3)/6;
    }
  
    return sum;
}
 

改寫成CUDA程式：

__device__
double integrant_func(double t)
{
    return 1.0/(exp(30.1*t*t/(1.23*1.23*300.0))-1);
}

__global__
void integrant(int n, double a, double dt, double *sum)
{
    double tk;
    double f1, f2, f3;

    int index = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    for(int i = index; i < n; i += stride)
    {
        tk=a+i*dt;
        f1=integrant_func(tk );
        f2=integrant_func(tk+dt/2);
        f3=integrant_func(tk+dt);
        sum[i]=dt*(f1+4*f2+f3)/6;
    }

}

double Simpson_integ (int n_steps, double a, double b, void (*func)(int n, double a, double dt, double *sum))
{
    double *sum;
    double dt;
    
    cudaMallocManaged(&sum, n_steps * sizeof(double));

    dt=(b-a)/n_steps;
    
    int blockSize = 256;
    int numBlock = (n_steps + blockSize - 1) / blockSize;

    integrant<<<numBlock, blockSize>>>(n_steps, a, dt, sum);
        
    cudaDeviceSynchronize();
    
    double Sum=0.;
    for(int i = 0; i < n_steps; i++)
        Sum += sum[i];

    cudaFree(sum);
    return Sum;
}
 

其中device函式就是我們要計算的被積函式，而global宣告的就是被積函式對應的kernel函式。

從主函式中呼叫：

include <iostream>
#include <sys/time.h>

#include "Integ.h"

using namespace std;

int main()
{
	int steps=1e6;
    double sum;
	struct timeval start,end;
	long timeuse;
	
	for(int i=0;i<10;i++)
	{
    	gettimeofday(&start, NULL);
        sum = Simpson_integ(steps, 1e-4, 3.1415926, integrant);
        gettimeofday(&end, NULL);
    	timeuse = 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec-start.tv_usec;
    	cout<<"GPU 計算結果："<<sum<<"	用時:"<<timeuse/1e6<<"s"<<endl;

    }
    return 0;
}
 

可能會問為什麼要執行10遍。其實是因為第一次執行時和原來利用CPU計算時做個了個比較，發現GPU計算速度在步數取時居然還沒有CPU算的快，大驚之下感覺又運行了幾次，發現一個很奇怪的現象，GPU做第一次計算時的速度是比較慢的，但第二次開始就十分快了，難道GPU和汽車一樣還有一個起步加速階段？

前面是用了一個GPU的程式，因為手邊可利用的GPU有4個，所以重新寫了一下：

double Simpson_integ (int n_steps, double a, double b, void (*func)(int n, double a, double dt, double *sum))
{
    int N_GPU = 4;

    TGPUplan plan[N_GPU];

    int np = (n_steps + N_GPU - 1) / N_GPU;

    double dt;

    dt=(b-a)/n_steps;

    for(int i = 0; i < N_GPU; i++)
    {
        cudaSetDevice(i);
        cudaStreamCreate(&plan[i].stream);

        cudaMalloc((void **)&plan[i].d_sum, np * sizeof(double));
        plan[i].h_sum = (double *)malloc(np * sizeof(double));
    }

    int blockSize = 256;
    int numBlock = (np + blockSize - 1) / blockSize;

    for(int i = 0; i < N_GPU; i++)
    {
        cudaSetDevice(i);

        integrant<<<numBlock, blockSize, 0, plan[i].stream>>>(np, a + i * dt * np, dt, plan[i].d_sum);
        
        cudaMemcpyAsync(plan[i].h_sum, plan[i].d_sum, np * sizeof(double), cudaMemcpyDeviceToHost, plan[i].stream);
    }
    
    
    
    double Sum=0.;
    for(int i = 0; i < N_GPU; i++)
    {
        cudaSetDevice(i);

        cudaStreamSynchronize(plan[i].stream);

        for(int j = 0; j < np; j++)
            Sum += plan[i].h_sum[j];

        cudaFree(plan[i].d_sum);
        free(plan[i].h_sum);

        cudaStreamDestroy(plan[i].stream); //Destroy the stream
    }
}

注意其中在執行kernel函式時初始值不再是，而是。

就執行結果來說，使用4個GPU後第一次加速很明顯，但隨後的幾次用時和1個GPU加速後用時相差不大，但都只有CPU用時的十分之一左右。總體來說還算是理想的。可能更負責的計算才能體現出GPU加速的優越性吧。