OpenCL

  // create a compute context with GPU device  context = clCreateContextFromType(NULL, CL_DEVICE_TYPE_GPU, NULL, NULL, NULL);  // create a command queue  queue = clCreateCommandQueue(context, NULL, 0, NULL);  // allocate the buffer memory objects  memobjs = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float)*2*num_entries, srcA, NULL);  memobjs = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float)*2*num_entries, NULL, NULL);  // create the compute program  program = clCreateProgramWithSource(context, 1, &fft1D_1024_kernel_src, NULL, NULL);  // build the compute program executable  clBuildProgram(program, 0, NULL, NULL, NULL, NULL);  // create the compute kernel  kernel = clCreateKernel(program, "fft1D_1024", NULL);  // set the args values  clSetKernelArg(kernel, 0, sizeof(cl_mem),(void *)&memobjs);  clSetKernelArg(kernel, 1, sizeof(cl_mem),(void *)&memobjs);  clSetKernelArg(kernel, 2, sizeof(float)*(local_work_size+1)*16, NULL);  clSetKernelArg(kernel, 3, sizeof(float)*(local_work_size+1)*16, NULL);  // create N-D range object with work-item dimensions and execute kernel  global_work_size = num_entries;  local_work_size = 64;  clEnqueueNDRangeKernel(queue, kernel, 1, NULL, global_work_size, local_work_size, 0, NULL, NULL);

  // This kernel computes FFT of length 1024. The 1024 length FFT is decomposed into  // calls to a radix 16 function, another radix 16 function and then a radix 4 function  __kernel void fft1D_1024(__global float2 *in, __global float2 *out,                          __local float *sMemx, __local float *sMemy){    int tid = get_local_id(0);    int blockIdx = get_group_id(0) * 1024 + tid;    float2 data;    // starting index of data to/from global memory    in = in + blockIdx;  out = out + blockIdx;    globalLoads(data, in, 64); // coalesced global reads    fftRadix16Pass(data);      // in-place radix-16 pass    twiddleFactorMul(data, tid, 1024, 0);    // local shuffle using local memory    localShuffle(data, sMemx, sMemy, tid, (((tid & 15)* 65) +(tid >> 4)));    fftRadix16Pass(data);               // in-place radix-16 pass    twiddleFactorMul(data, tid, 64, 4); // twiddle factor multiplication    localShuffle(data, sMemx, sMemy, tid, (((tid >> 4)* 64) +(tid & 15)));    // four radix-4 function calls    fftRadix4Pass(data);      // radix-4 function number 1    fftRadix4Pass(data + 4);  // radix-4 function number 2    fftRadix4Pass(data + 8);  // radix-4 function number 3    fftRadix4Pass(data + 12); // radix-4 function number 4    // coalesced global writes    globalStores(data, out, 64);  }

Apple的网站上可以发现傅立叶变换的例子

使用 Python 3.x 搭配 PyOpenCL 与 NumPy

import ioimport randomimport numpy as npimport pyopencl as cldef dump_step(data, chunk_size):    """顯示排序過程"""    msg = io.StringIO('')    div = io.StringIO('')    for idx, item in enumerate(data):        if idx % chunk_size == 0:            if idx > 0:                msg.write(' ||')                div.write('   ')            div.write(' --')        else:            msg.write('   ')            div.write('------')        msg.write(' {:2d}'.format(item))    out = msg.getvalue()    if chunk_size == 1: print(' ' + '-' * (len(out) - 1))    print(out)    print(div.getvalue())    msg.close()    div.close()def cl_merge_sort_sbs(data_in):    """平行合併排序"""    # OpenCL kernel 函數程式碼    CL_CODE = '''    kernel void merge(int chunk_size, int size, global long* data, global long* buff) {        // 取得分組編號        const int gid = get_global_id(0);        // 根據分組編號計算責任範圍        const int offset = gid * chunk_size;        const int real_size = min(offset + chunk_size, size) - offset;        global long* data_part = data + offset;        global long* buff_part = buff + offset;        // 設定合併前的初始狀態        int r_beg = chunk_size >> 1;        int b_ptr = 0;        int l_ptr = 0;        int r_ptr = r_beg;        // 進行合併        while (b_ptr < real_size) {            if (r_ptr >= real_size) {                // 若右側沒有資料，取左側資料堆入緩衝區                buff_part = data_part;            } else if (l_ptr == r_beg) {                // 若左側沒有資料，取右側資料堆入緩衝區                buff_part = data_part;            } else {                // 若兩側都有資料，取較小資料堆入緩衝區                if (data_part < data_part) {                    buff_part = data_part;                } else {                    buff_part = data_part;                }            }            b_ptr++;        }    }    '''    # 配置計算資源，編譯 OpenCL 程式    ctx = cl.Context(dev_type=cl.device_type.GPU)    prg = cl.Program(ctx, CL_CODE).build()    queue = cl.CommandQueue(ctx)    mf = cl.mem_flags    # 資料轉換成 numpy 形式以利轉換為 OpenCL Buffer    data_np = np.int64(data_in)    buff_np = np.empty_like(data_np)    # 建立緩衝區，並且複製數值到緩衝區    data = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=data_np)    buff = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=buff_np)    # 設定合併前初始狀態    data_len = np.int32(len(data_np))    chunk_size = np.int32(1)    dump_step(data_np, chunk_size)    while chunk_size < data_len:        # 更新分組大小，每一回合變兩倍        chunk_size <<= 1        # 換算平行作業組數         group_size = ((data_len - 1) // chunk_size) + 1        # 進行分組合併作業        prg.merge(queue, (group_size,), (1,), chunk_size, data_len, data, buff)        # 將合併結果作為下一回合的原始資料        temp = data        data = buff        buff = temp        # 顯示此回合狀態        cl.enqueue_copy(queue, data_np, data)        dump_step(data_np, chunk_size)    queue.finish()    data.release()    buff.release()def main():    n = random.randint(5, 16)    data =     for i in range(n):        data.append(random.randint(1, 99))    cl_merge_sort_sbs(data)if __name__ == '__main__':    main()