記一次CUDA編程任務

　　

　　這個月6號開始，着手解決一個具有實際意義的計算任務。任務數據有9879896條，每條包含30個整數，任務是計算每兩條數據之間的斯皮爾相關系數及其P值。原始數據只有500+MB，因此我並不認為這是個多么大的計算任務。隨后稍加計算，我還是很驚呆的，要計算(9879896×9879895)÷2≈4.88億億組數據，但此時這還只是個數字概念，我也沒意識到時間復雜度和空間復雜度的問題。

 

1. 計算規模初體驗

數據格式：9879896行，30列，每列之間以空格符隔開，例如：

0 2 0 2 0 0 0 0 0 0 0 40 0 0 35 0 0 53 0 44 0 0 0 0 0 0 0 0 0 0
0 0 1 148 0 0 0 0 0 0 0 0 0 0 1133 0 1 0 0 1820 0 0 0 2 0 0 0 1 0 0
0 0 0 33 1 0 0 0 0 0 0 0 0 0 231 0 0 0 0 402 0 0 0 0 0 0 0 0 0 0
0 0 6 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 6 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
... ...
... ...

空間復雜度：單純計算下結果大概有多大吧，每組計算結果包含相關系數和P值，若都以float（占4字節）精度存儲，需要占用內存：4.88億億×8B≈400TB，當然，我們不具備這么大內存，因此無論以何種方式計算，都需要一批批地重復將數據載入內存、計算、存入硬盤這個過程，直到運算完成。那么，存入硬盤的結果會有400TB嗎？不然，P值小於或等於0.05的結果才會需要輸出，因此實際上會遠遠小於這個值，具體會小多少，先運行一批數據后才能做出估算。

時間復雜度：計算的組數規模是(n×(n-1))÷2，那么就看程序能跑多快了。我想先看看MATLAB多線程、Python多線程、Spark分布式計算能跑多快，是否能在最快時間內解決問題。

 

2. MATLAB多線程

　　MATLAB寫起來最簡單，計算相關系數和P值都不用操心，一行自帶的函數調用就完成。打開MATLAB左下角的並行池，MATLAB將會自動尋找到機子上有的物理核心，並分配與物理核心數相同的worker。比如我的電腦是4核8線程，它只能開4個worker，不識別虛擬核心。

　　代碼如下：

t1 = clock;
disp('>> loading ...');
A = importdata('D:/MASTER2016/5.CUDA/data-ID-top30-kv3.txt');
b = A'; %由於MATLAB只計算列與列之間的相關系數，因此需要轉置操作
disp(etime(clock,t1));

num = size(b, 2);
disp('>> calculating ...');
fid = fopen('D:/MASTER2016/5.CUDA/result-matlab.txt', 'wt');

for i = 1 : num
    for j = i+1 : num
        [m, n] = corr(b(:, i), b(:, j), 'type', 'Spearman', 'tail', 'both');
        if isnan(n) || n>0.05
            continue;
        end
		
        fprintf(fid, 'X%d\tX%d\t%d\t%d\n', i, j, m, n);
    end
end
fclose(fid);
disp('>> OK!');

　　這里我並沒有考慮內存空間不夠的問題，因為我只是想說明MATLAB的計算速度。開了多顆核心的情況下，MATLAB並沒能完全壓榨出所有的CPU性能，計算速度緩慢無比，更要命的是，它會越算越慢。據我估算，即使空間復雜度足夠，MATLAB也要用超過20年的時間才能算完，這還是不考慮越算越慢的情況。

　　好了，此方案僅是打醬油。

 

3. Python多線程

　　Python語言由於本身的體質問題，Cython下不能調用多核，只能用多線程。理論上是這樣，但還是有很多擴展包能夠充分壓榨出多核CPU性能，例如multiprocessing是其中的佼佼者。multiprocessing用起來也非常簡單，考慮到CPU的多核運算下，每顆核心的算力還是很可觀的，所有不能把每個計算組都拆成並行線程，那樣內存的讀寫開銷反而會使CPU一直在等待狀態，不能一直滿負載工作。鑒於此，我設計9879895組線程，每組代表某個特定行與剩下的各個數據行形成的數據組。這樣每組線程下的運算量還是比較大的，能使CPU盡可能全在滿負載狀態。

　　代碼如下：

# coding=utf-8

import math
import multiprocessing
import time

import scipy.stats as stats


def calculate2(i, X, all_glb, data_array_glb):
    all = all_glb.value
    result = []
    for j in range(i + 1, all):
        x = X
        y = data_array_glb[j]
        if math.fsum(x) == 0 or math.fsum(y) == 0:
            continue
        corr, p = stats.spearmanr(x, y)
        if p > 0.05:
            continue
        result.append([i + 1, j + 1, corr, p])
    return result


if __name__ == "__main__":

    multiprocessing.freeze_support()

    input_file = 'D:/MASTER2016/5.CUDA/data-ID-top30-kv3.txt'
    output_file = 'D:/MASTER2016/5.CUDA/result-python.txt'

    print '>> loading ...'
    start = time.clock()
    data = open(input_file)
    data_array = []
    for line in data:
        data_array.append(map(int, line.strip().split(' ')))
    data.close()
    print time.clock()-start, 's'

    print '>> calculating ...'
    results = []
    pool_size = 8
    pool = multiprocessing.Pool(processes=pool_size)
    all = len(data_array)
    manager = multiprocessing.Manager()
    all_share = manager.Value('i', int(all))
    data_array_share = manager.list(data_array)
    for i in range(all):
        data_X = data_array[i]
        results.append(pool.apply_async(calculate2, args=(i, data_X, all_share, data_array_share)))
    pool.close()
    pool.join()
    print time.clock() - start, 's'
    data_array = None

    print '>> saving ...'
    data2 = open(output_file, 'w')
    for res in results:
        temp_list = res.get()
        for temp in temp_list:
            data2.write('X'+str(temp[0])+'\t'+'X'+str(temp[1])+'\t'+str(temp[2])+'\t'+str(temp[3])+'\n')
    print time.clock()-start, 's'
    data2.close()

　　

 

      這里，我依然沒有考慮空間復雜度問題，因為要先看看計算能力是否能滿足任務要求。Python的這個多線程下，確實能充分榨干CPU性能，風扇呼呼響，要命的是也存在越算越慢的問題。但是，即使CPU一直這么滿負載運算，我粗略估算了下，也得要個14年+才能算完，也不算越算越慢的情況。

　　所以，此方案是打醬油2號。

 

4. Spark方案

　　Spark方案我並沒有寫完，因為寫着寫着就感覺到。。。肯定還是不行，CPU的算力也就那樣了。就算調12台機器一起跑，也不適合用CPU下的線程模型解決問題了。

這種高並行的計算，要想取得最快計算速度，非GPU莫屬。

 

5. CUDA方案

　　CUDA方案下，首先必須清晰地設計好線程模型，即：我需要用到幾塊GPU？我需要在每塊GPU上設計多少個block？每個block設計多少個線程？每個線程分配多少運算量？這四個問題基本決定了CUDA程序的性能和復雜度。

　　CUDA是一種異構並行解決方案，即CPU用於控制，GPU用於主運算的方案。一個GPU有一個grid，每個grid里有大量block，每個block里有大量thread。在運算時，每個thread都是完全獨立並行地運算，每個線程里的運算靠內核函數控制，這也是CUDA編程的核心，目前只能用CUDA C編寫。因此JCUDA和PyCUDA做的只是內存分配這些CPU端控制的事情，還不能代替GPU端的CUDA C代碼。

 

　　如上圖，左邊列是Host端，即CPU上執行的控制端，用於分配GPU內存空間，拷貝內存數據到GPU顯存等等操作。右邊列是Device端，即GPU上的並行模型，由grid，block，thread三者構成。不同型號GPU的最大block數和每個block中的最大thread不同，但是可以查詢。在安裝好CUDA Toolkit后，windows用戶可以進入C:\ProgramData\NVIDIA Corporation\CUDA Samples\v8.0\1_Utilities\deviceQuery目錄，打開相應版本的項目，執行運行查詢。

　　比如我的機器：

　　基於此，我設計的線程模型是：比如數據是ROWS行，COLS列，那么有((ROWS-1)×ROWS)÷2組計算，每一行都要與從這行開始后面的每一行進行計算。開辟(ROWS-1)個block，編號0~(ROWS-1)對應着數據的行號。所以，對於第一行，行號是0，要與1~(ROWS-1)的每一行進行計算，一共有(ROWS-1)組，這些計算任務分配給第一塊block的1024個線程上計算。依此類推。這樣做並不是最佳的任務分配方案，因為不是公平分配，編號越靠后的block分配的任務越少。但是，這樣做的好處是便於利用共享內存，加速每一個block內的計算。

　　比如第一行，將數據第一行存入共享內存，那么它在與其他行分別計算的時候，直接從每個block內的共享內存讀取數據，遠遠比從顯存上的全局內存讀取速度快得多。需要注意的是，每塊block內的共享內存的大小也有硬件限制，上面截圖中可以看到，GTX 950M的共享內存是49152B。

　　Talk is cheap. Show me the code:

#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

// 定義總數據矩陣的行數和列數
#define ROWS 15000
#define COLS 30

// 定義每一塊內的線程個數，GT720最多是1024（必須大於總矩陣的列數：30）
#define NUM_THREADS 1024


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;
}

__device__ float meanForRankCUDA(int num)
{
    float sum = 0;
    for (int i = 0; i <= num; i++) {
        sum += i;
    }
    return sum / (num + 1);
}


__device__ float meanForArrayCUDA(float array[], int len)
{
    float sum = 0;
    for (int i = 0; i < len; i++) {
        sum += array[i];
    }
    return sum / len;
}


__device__ float spearmanKernel(int Xarray[], int Yarray[])
{
    //1，對原先的數據進行排序，相同的值取平均值
    float Xrank[30];
    float Yrank[30];
    int col = 30;

    for (int i = 0; i < col; i++) {
        int bigger = 1;
        int equaer = -1;
        for (int j = 0; j < col; j++) {
            if (Xarray[i] < Xarray[j]) {
                bigger = bigger + 1;
            }
            else if (Xarray[i] == Xarray[j]) {
                equaer = equaer + 1;
            }
        }
        Xrank[i] = bigger + meanForRankCUDA(equaer);
    }
    for (int i = 0; i < col; i++) {
        int bigger = 1;
        int equaer = -1;
        for (int j = 0; j < col; j++) {
            if (Yarray[i] < Yarray[j]) {
                bigger = bigger + 1;
            }
            else if (Yarray[i] == Yarray[j]) {
                equaer = equaer + 1;
            }
        }
        Yrank[i] = bigger + meanForRankCUDA(equaer);
    }

    //2，計算斯皮爾曼相關性系數
    float numerator = 0;
    float denominatorLeft = 0;
    float denominatorRight = 0;
    float meanXrank = meanForArrayCUDA(Xrank, col);
    float meanYrank = meanForArrayCUDA(Yrank, col);
    for (int i = 0; i < col; i++) {
        numerator += (Xrank[i] - meanXrank) * (Yrank[i] - meanYrank);
        denominatorLeft += powf(Xrank[i] - meanXrank, 2);
        denominatorRight += powf(Yrank[i] - meanYrank, 2);
    }
    float corr = 0;
    if ((denominatorLeft != 0) && (denominatorRight != 0)) {
        corr = numerator / sqrtf(denominatorLeft * denominatorRight);
    }
    return corr;
}


__global__ static void spearCUDAShared(const int* a, size_t lda, float* c, size_t ldc, float* d, size_t ldd)
{
    extern __shared__ int data[];
    const int tid = threadIdx.x;
    const int row = blockIdx.x;
    int i, j;
    // 同步第1行~倒數第二行到共享內存，行數由block個數（總數據矩陣的行數-1）控制，每個block共享一行數據
    if (tid < 30) {
        data[tid] = a[row * lda + tid];
    }
    __syncthreads();

    int cal_per_block = gridDim.x - row; // 每個塊分擔的計算量
    int cal_per_thread = cal_per_block / blockDim.x + 1; // 每個線程分擔的計算量
    // 分配各線程計算任務，通過for循環控制在一個線程需要計算的組數
    for (i = row + cal_per_thread * tid; i < (row + cal_per_thread * (tid + 1)) && i < gridDim.x; i++) {
        int j_row[30]; // 存放總數據矩陣的第j行
        for (j = 0; j < 30; j++) {
            j_row[j] = a[(i + 1)*lda + j];
        }
        float corr = spearmanKernel(data, j_row);
        c[row * ldc + (i + 1)] = corr;
        float t_test = 0;
        if (corr != 0) t_test = corr*(sqrtf((30 - 2) / (1 - powf(corr, 2))));
        d[row * ldd + (i + 1)] = t_test;
        //printf("block號：%d, 線程號：%d, 計算組：%d-%d, id號：%d, block個數：%d, 每塊線程個數：%d, 該塊總計算量：%d, 該塊中每個線程計算量：%d, corr: %lf, %d, %d, %d - %d, %d, %d\n", row, tid, row, i + 1, (row*blockDim.x + tid), gridDim.x, blockDim.x, cal_per_block, cal_per_thread, corr, data[0], data[1], data[29], j_row[0], j_row[1], j_row[29]);
    }
}


clock_t matmultCUDA(const int* a, float* c, float* d)
{
    int *ac;
    float *cc, *dc;
    clock_t start, end;
    start = clock();

    size_t pitch_a, pitch_c, pitch_d;
    // 開辟a、c、d在GPU中的內存
    cudaMallocPitch((void**)&ac, &pitch_a, sizeof(int)* COLS, ROWS);
    cudaMallocPitch((void**)&cc, &pitch_c, sizeof(float)* ROWS, ROWS);
    cudaMallocPitch((void**)&dc, &pitch_d, sizeof(float)* ROWS, ROWS);
    // 復制a從CPU內存到GPU內存
    cudaMemcpy2D(ac, pitch_a, a, sizeof(int)* COLS, sizeof(int)* COLS, ROWS, cudaMemcpyHostToDevice);

    spearCUDAShared << <ROWS - 1, NUM_THREADS, sizeof(int)* COLS >> > (ac, pitch_a / sizeof(int), cc, pitch_c / sizeof(float), dc, pitch_d / sizeof(float));

    cudaMemcpy2D(c, sizeof(float)* ROWS, cc, pitch_c, sizeof(float)* ROWS, ROWS, cudaMemcpyDeviceToHost);
    cudaMemcpy2D(d, sizeof(float)* ROWS, dc, pitch_d, sizeof(float)* ROWS, ROWS, cudaMemcpyDeviceToHost);
    cudaFree(ac);
    cudaFree(cc);

    end = clock();
    return end - start;
}


void print_int_matrix(int* a, int row, int col) {
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            printf("%d\t", a[i * col + j]);
        }
        printf("\n");
    }
}


void print_float_matrix(float* c, int row, int col) {
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            printf("%f\t", c[i * col + j]);
        }
        printf("\n");
    }
}

void read_ints(int* a) {
    FILE* file = fopen("D:\\MASTER2016\\5.CUDA\\data-ID-top30-kv.txt", "r");
    int i = 0;
    int count = 0;

    fscanf(file, "%d", &i);
    while (!feof(file))
    {
        a[count] = i;
        count++;
        if (count == ROWS*COLS) break;
        fscanf(file, "%d", &i);
    }
    fclose(file);
}


int main()
{
    int *a; // CPU內存中的總數據矩陣，ROWS行，COLS列
    float *c; // CPU內存中的相關系數結果矩陣，ROWS行，ROWS列
    float *d; // CPU內存中的T值結果矩陣，ROWS行，ROWS列
    a = (int*)malloc(sizeof(int)* COLS * ROWS);
    c = (float*)malloc(sizeof(float)* ROWS * ROWS);
    d = (float*)malloc(sizeof(float)* ROWS * ROWS);

    clock_t start = clock();
    printf(">> loading ... rows: %d, cols: %d", ROWS, COLS);
    read_ints(a);
    clock_t end = clock() - start;
    printf("\nTime used: %.2f s\n", (double)(end) / CLOCKS_PER_SEC);

    //print_int_matrix(a, ROWS, COLS);
    //printf("\n");

    printf(">> calculating ... ");
    printf("\n---------------------------------------");
    printf("\ntotal groups: %lld", (long long)ROWS*(ROWS - 1) / 2);
    printf("\ntotal threads: %d (blocks) * 1024 = %d", (ROWS - 1), (ROWS - 1) * 1024);
    printf("\ntotal space complexity: %lld MB", (long long)((ROWS / 1024) * (ROWS / 1024) * 8));
    printf("\n---------------------------------------");
    if (!InitCUDA()) return 0;
    clock_t time = matmultCUDA(a, c, d);
    double sec = (double)(time + end) / CLOCKS_PER_SEC;
    printf("\nTime used: %.2f s\n", sec);

    printf(">> saving ... ");
    FILE *f = fopen("D:\\MASTER2016\\5.CUDA\\result-c-2.txt", "w");
    for (int i = 0; i < ROWS; i++) {
        for (int j = i + 1; j < ROWS; j++) {
            float t_test = d[i * ROWS + j];
            if (t_test >= 2.042) {
                fprintf(f, "X%d\tX%d\t%f\t%lf\n", i + 1, j + 1, c[i * ROWS + j], t_test);
            }
        }
    }
    fclose(f);
    end = clock() - start;
    printf("OK\nTime used: %.2f s\n", (double)(end) / CLOCKS_PER_SEC);

    //printf(">> 相關系數結果矩陣: \n");
    //print_float_matrix(c, ROWS, ROWS);
    //printf(">> T值結果矩陣: \n");
    //print_float_matrix(d, ROWS, ROWS);

    getchar();
    return 0;
}

 

　　需要指出的是，上面程序保存為filename.cu文件，執行nvcc -o filename filename.cu編譯，執行filename即可運行。其中ROWS是從總數據文件中讀取的行數，用於控制數據規模調試程序，如果ROWS大於或等於總數據行數，那么就是讀取整個文件了。

　　由於空間復雜度太高，也就是最開始提到的，那么下面做些調整，加個控制參數，每次只計算一定的行數，使顯存滿載但不超出即可。相應地，內核函數中的索引號，保存文件的函數都需要做些微調，代碼如下：

#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "cuda_runtime.h"  
#include "device_launch_parameters.h"

// 定義總數據矩陣的行數和列數
#define ROWS 1000
#define COLS 30

// 控制一次計算占用顯存的大小：CONTROL_ROWS*ROWS*8（字節）< 顯存
#define CONTROL_ROWS 45

// 定義每一塊內的線程個數，GT720最多是1024
#define NUM_THREADS 1024


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;
}

__device__ float meanForRankCUDA(int num)
{
    float sum = 0;
    for (int i = 0; i <= num; i++) {
        sum += i;
    }
    return sum / (num + 1);
}


__device__ float meanForArrayCUDA(float array[], int len)
{
    float sum = 0;
    for (int i = 0; i < len; i++) {
        sum += array[i];
    }
    return sum / len;
}


__device__ float spearmanKernel(int Xarray[], int Yarray[])
{
    //1，對原先的數據進行排序，相同的值取平均值
    float Xrank[30];
    float Yrank[30];
    int col = 30;

    for (int i = 0; i < col; i++) {
        int bigger = 1;
        int equaer = -1;
        for (int j = 0; j < col; j++) {
            if (Xarray[i] < Xarray[j]) {
                bigger = bigger + 1;
            }
            else if (Xarray[i] == Xarray[j]) {
                equaer = equaer + 1;
            }
        }
        Xrank[i] = bigger + meanForRankCUDA(equaer);
    }
    for (int i = 0; i < col; i++) {
        int bigger = 1;
        int equaer = -1;
        for (int j = 0; j < col; j++) {
            if (Yarray[i] < Yarray[j]) {
                bigger = bigger + 1;
            }
            else if (Yarray[i] == Yarray[j]) {
                equaer = equaer + 1;
            }
        }
        Yrank[i] = bigger + meanForRankCUDA(equaer);
    }

    //2，計算斯皮爾曼相關性系數
    float numerator = 0;
    float denominatorLeft = 0;
    float denominatorRight = 0;
    float meanXrank = meanForArrayCUDA(Xrank, col);
    float meanYrank = meanForArrayCUDA(Yrank, col);
    for (int i = 0; i < col; i++) {
        numerator += (Xrank[i] - meanXrank) * (Yrank[i] - meanYrank);
        denominatorLeft += powf(Xrank[i] - meanXrank, 2);
        denominatorRight += powf(Yrank[i] - meanYrank, 2);
    }
    float corr = 0;
    if ((denominatorLeft != 0) && (denominatorRight != 0)) {
        corr = numerator / sqrtf(denominatorLeft * denominatorRight);
    }
    return corr;
}


__global__ static void spearCUDAShared(const int* a, size_t lda, float* c, size_t ldc, float* d, size_t ldd, int cols, int start)
{
    extern __shared__ int data[];
    const int tid = threadIdx.x;
    const int row = blockIdx.x;

    int i, j;
    // 同步第1行~倒數第二行到共享內存，行數由block個數控制，每個block共享一行數據
    if (tid < 30) {
        data[tid] = a[(start + row) * lda + tid];
    }
    __syncthreads();

    int cal_per_block = cols - (start + row); // 每個塊分擔的計算量
    int cal_per_thread = cal_per_block / blockDim.x + 1; // 每個線程分擔的計算量
    // 分配各線程計算任務，通過for循環控制在一個線程需要計算的組數
    for (i = row + cal_per_thread * tid; i < (row + cal_per_thread * (tid + 1)) && i < cols; i++) {
        int j_row[30]; // 存放總數據矩陣的第j行
        for (j = 0; j < 30; j++) {
            j_row[j] = a[(start + i + 1)*lda + j];
        }
        float corr = spearmanKernel(data, j_row);
        c[row * ldc + (start + i + 1)] = corr;
        float t_test = 0;
        if (corr != 0) t_test = corr*(sqrtf((30 - 2) / (1 - powf(corr, 2))));
        d[row * ldd + (start + i + 1)] = t_test;
        //printf("block號：%d, 線程號：%d, 計算組：%d-%d, id號：%d, block個數：%d, 每塊線程個數：%d, 該塊總計算量：%d, 該塊中每個線程計算量：%d, corr: %lf, %d, %d, %d - %d, %d, %d\n", row, tid, row, i + 1, (row*blockDim.x + tid), gridDim.x, blockDim.x, cal_per_block, cal_per_thread, corr, data[0], data[1], data[29], j_row[0], j_row[1], j_row[29]);
    }
}


clock_t matmultCUDA(const int* a, float* c, float* d, int start_index, int control_rows)
{
    int *ac;
    float *cc, *dc;
    clock_t start, end;
    start = clock();

    size_t pitch_a, pitch_c, pitch_d;
    // 開辟a、c、d在GPU中的內存
    cudaMallocPitch((void**)&ac, &pitch_a, sizeof(int)* COLS, ROWS);
    cudaMallocPitch((void**)&cc, &pitch_c, sizeof(float)* ROWS, control_rows);
    cudaMallocPitch((void**)&dc, &pitch_d, sizeof(float)* ROWS, control_rows);
    // 復制a從CPU內存到GPU內存
    cudaMemcpy2D(ac, pitch_a, a, sizeof(int)* COLS, sizeof(int)* COLS, ROWS, cudaMemcpyHostToDevice);

    spearCUDAShared << <control_rows, NUM_THREADS, sizeof(int)* COLS >> > (ac, pitch_a / sizeof(int), cc, pitch_c / sizeof(float), dc, pitch_d / sizeof(float), ROWS - 1, start_index);
    
    cudaMemcpy2D(c, sizeof(float)* ROWS, cc, pitch_c, sizeof(float)* ROWS, control_rows, cudaMemcpyDeviceToHost);
    cudaMemcpy2D(d, sizeof(float)* ROWS, dc, pitch_d, sizeof(float)* ROWS, control_rows, cudaMemcpyDeviceToHost);
    cudaFree(ac);
    cudaFree(cc);
    cudaFree(dc);

    end = clock();
    return end - start;
}


void print_int_matrix(int* a, int row, int col) {
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            printf("%d\t", a[i * col + j]);
        }
        printf("\n");
    }
}


void print_float_matrix(float* c, int row, int col) {
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            printf("%f\t", c[i * col + j]);
        }
        printf("\n");
    }
}

void read_ints(int* a, char *input_file) {
    FILE* file = fopen(input_file, "r");
    int i = 0;
    int count = 0;

    fscanf(file, "%d", &i);
    while (!feof(file))
    {
        a[count] = i;
        count++;
        if (count == ROWS*COLS) break;
        fscanf(file, "%d", &i);
    }
    fclose(file);
}

void clear_ints(char * out_file) {
    FILE *f = fopen(out_file, "w");
    fclose(f);
}

void cal_and_save(int i, int *a, char *out_file, int control_rows) {
    float *c; // CPU內存中的相關系數結果矩陣，ROWS行，ROWS列
    float *d; // CPU內存中的T值結果矩陣，ROWS行，ROWS列
    c = (float*)malloc(sizeof(float)* control_rows * ROWS);
    d = (float*)malloc(sizeof(float)* control_rows * ROWS);

    clock_t time = matmultCUDA(a, c, d, i, control_rows);

    FILE *f = fopen(out_file, "a");
    for (int m = 0; m < control_rows; m++) {
        for (int n = i + m + 1; n < ROWS; n++) {
            float t_test = d[m * ROWS + n];
            if (t_test >= 2.042) {
                fprintf(f, "X%d\tX%d\t%f\t%lf\n", i + m + 1, n + 1, c[m * ROWS + n], t_test);
            }
        }
    }
    fclose(f);

    //printf(">> 相關系數結果矩陣: \n");
    //print_float_matrix(c, CONTROL_ROWS, ROWS);
    //printf(">> T值結果矩陣: \n");
    //print_float_matrix(d, CONTROL_ROWS, ROWS);

    free(c);
    free(d);
}

int main()
{
    int *a; // CPU內存中的總數據矩陣，ROWS行，COLS列
    a = (int*)malloc(sizeof(int)* COLS * ROWS);

    char *input_file = "D:\\MASTER2016\\5.CUDA\\data-ID-top30-kv.txt";
    char *out_file = "D:\\MASTER2016\\5.CUDA\\result-c.txt";

    clock_t start = clock();
    printf(">> loading ... rows: %d, cols: %d", ROWS, COLS);
    read_ints(a, input_file);
    clear_ints(out_file);
    clock_t end = clock() - start;
    printf("\nTime used: %.2f s\n", (double)(end) / CLOCKS_PER_SEC);

    //print_int_matrix(a, ROWS, COLS);
    //printf("\n");

    printf(">> calculating ... ");
    printf("\n---------------------------------------");
    printf("\ntotal groups: %lld", (long long)ROWS*(ROWS - 1) / 2);
    printf("\ntotal threads: %d (blocks) * 1024 = %d", (ROWS - 1), (ROWS - 1) * 1024);
    printf("\ntotal space complexity: %lld MB", (long long)((CONTROL_ROWS / 1024) * (ROWS / 1024) * 8));
    printf("\n---------------------------------------");

    if (!InitCUDA()) return 0;

    int i;
    for (i = 0; i < ROWS - 1; i += CONTROL_ROWS) {
        printf("\n>> calculating and saving ... id: %d ... ", i);
        cal_and_save(i, a, out_file, CONTROL_ROWS);
        end = clock() - start;
        printf("Time used: %.2f s", (double)(end) / CLOCKS_PER_SEC);
    }

    // 不能整除的非整數部分需要計算
    //i -= CONTROL_ROWS;
    //int control_rows = ROWS - 1 - i;
    //printf("\n%d", control_rows);
    //if (control_rows > 0) {
    //    printf("\n>> calculating and saving ... id: %d ... ", i);
    //    cal_and_save(i, a, out_file, control_rows);
    //    end = clock() - start;
    //    printf("Time used: %.2f s", (double)(end) / CLOCKS_PER_SEC);
    //}

    printf("\nFinished.\n");

    getchar();
    return 0;
}

 

　　到現在，由於空間復雜度過高而顯存不夠的問題通過增加時間復雜度的方法基本解決了。當然在顯存足夠的情況下，還是一次性算完是最快的，實測CUDA提速100+倍，數據量越大提速越明顯。原因一是你必須被逼着按照CUDA的並行模型來寫程序，二是GPU的架構設計確實更適合超大並行程序的加速。游戲畫面渲染就是這樣，你可以想成一個block控制一塊屏幕的渲染，每塊block的每個線程控制幾個像素格的渲染，而這些圖像渲染完全可以是獨立並行的，GPU的設計初衷，即是增加核心數，不玩命升頻率，增加顯存帶寬，使成百上千的核心數的並行計算能力得到充分釋放。

　　但是目前的程序當然也不是完美的，我沒有考慮如何隱藏內存與顯存之間數據的傳輸延遲，沒有考慮多塊GPU如何聯動運算。后面我會思考這些。

　　並行計算是計算的未來。異構並行計算，也將是所有架構師必須增加的學習庫。