Уменьшение времени работы ядра cuda: динамическое распределение памяти в матрицах ядра

Question

Я хочу выполнить подгонку OLS для очень большого числа матриц меньшего размера, параллельно выполняя матричные операции на графическом процессоре. Я написал код, который, кажется, работает, но он медленнее, чем ожидалось. В настоящее время требуется меньше времени, чтобы запустить его в одном потоке на CPU, несмотря на параллельные вычисления на GPU. Nvidia Visual Profiler, похоже, указывает на то, что выделение памяти занимает много времени. Я подозреваю, что виновником является динамическое распределение памяти в матрицах разных размеров внутри ядра. Мне нужен совет и помощь по ускорению работы ядра.

Я попытался использовать new и delete для каждой матрицы, созданной в цикле.

Вот ядро:

__global__
void comb_ols(double *y, double *X, double *R2 ,const unsigned int M, const unsigned int N, int* sub_col, int *sub_size, int* cumulative_size, const unsigned int numberOfCalculations){

    int size;   
    int start_index;

    int index = blockIdx.x*blockDim.x+threadIdx.x;
    int stride = blockDim.x*gridDim.x;  
    for(int i = index; i < numberOfCalculations; i+=stride){    

        size = sub_size[i];
        start_index = cumulative_size[i];             

        double *sub_matrix = new double[M*(1+size)];


            for(int j = 0; j < size; j++){
            for(int k  = 0; k<M; k++){
                sub_matrix[k] = 1;
                sub_matrix[k + M * (1 +  j)] = X[k + M * (sub_col[start_index+j]+1)];                                           
                                            }       
            }
        }

        R2[i] = getR2(y,sub_matrix,M,size+1);


        delete [] sub_matrix;
    }
}

В функции устройства getR2 имеем следующее:

__device__
double getR2(double *y, double *X ,const unsigned int M, const unsigned int N) {

    // Initilize values
    double R2, numerator;
    double* A = new double[N*N];
    double* IA = new double[N*N];
    double* yX = new double[N];  
    // Generate all components
    XtX(X, A, M, N);
    LUPDecompose(A, N);
    LUPInvert(A, N, IA);
    yTX(y, X, yX, M, N);
    // Calc R2
    numerator = olsR2numerator(yX, IA, N);
    R2 = numerator / yTy(y, M);
    //R2 = yTy(y,M);

    delete[] A;
    delete[] IA;
    delete[] yX;

    return R2;
}

Фактический вызов ядра выглядит так:

com_ols<<<numBlocks, blockSize >>>(Y,X,R2,M,N,sub_columns, sub_size, cumulative_size, numberOfCalculations);

В настоящее время время работы ядра составляет примерно 1,4 секунды, тогда как в однопоточном процессоре оно составляет 0,7 секунды. Я ожидаю, что время выполнения ядра будет намного быстрее, так как оно только зацикливает много итераций матричных операций, которые должны быть подходящими для gpu. Есть что-то неэффективное в том, как распределяется память матриц разных размеров. Что вы, ребята, говорите о динамическом хранении матриц разных размеров внутри ядра? Как это сделать наиболее эффективным способом?

Любые другие отзывы о данном коде приветствуются.

Answer 1

Мне кажется, здесь применимы три очень простых эмпирических правила:

1. Динамическое выделение памяти всегда дорого, независимо от платформы, на которой вы программируете.
2. Код исполнителя никогда не использует динамическое выделение памяти, если только это не является абсолютно необходимым.
3. Если динамическое выделение памяти является абсолютно необходимым, предварительно выделите память и используйте ее как можно больше

Если вы посмотрите на свой код, он нарушает все эти три понятия.

Вы четко знаете (или могли бы просто рассчитать), какое максимальное значение sub_size перед запуском ядра. Используйте это априорное знание в своих интересах - предварительно выделите память кучи для вычислений, которая достаточно велика, чтобы обработать самую большую проблему в наборе данных, и повторно использовать ее в течение срока службы потока. Ваше ядро ​​может очень легко выглядеть примерно так:

__global__
void comb_ols(double *y, double *X, double *R2 ,const unsigned int M, 
             const unsigned int N, int* sub_col, int *sub_size, int* cumulative_size, 
             const unsigned int numberOfCalculations, const int max_size){

    int size;   
    int start_index;

    int index = blockIdx.x*blockDim.x+threadIdx.x;
    int stride = blockDim.x*gridDim.x;

    double *sub_matrix = new double[M*(1+max_size)];
    R2scratch temp(1+max_size);

    for(int i = index; i < numberOfCalculations; i+=stride){    

        size = sub_size[i];
        start_index = cumulative_size[i];             
        for(int j = 0; j < size; j++){
            for(int k  = 0; k<M; k++){
                sub_matrix[k] = 1;
                sub_matrix[k + M * (1 +  j)] = X[k + M * (sub_col[start_index+j]+1)];                                           
                                            }       
            }
        }
        R2[i] = getR2(y,sub_matrix,M,size+1,temp);
    }
    delete [] sub_matrix;
}

и устройство работает примерно так:

struct R2scratch
{
    double* A;
    double* IA;
    double* yX;  

    __device__
    R2scratch(int N) {
        A = new double[N*N];
        IA = new double[N*N];
        yX = new double[N];  
    };

    __device__
    ~R2scratch() {
        delete[] A;
        delete[] IA;
        delete[] yX;
    };
};

__device__
double getR2(double *y, double *X ,const unsigned int M, const unsigned int N, 
             R2scratch &scratch) {

    // Initilize values
    double R2, numerator;
    double* A = scratch.A;
    double* IA = scratch.IA;
    double* yX = scratch.yX;

    // Generate all components
    XtX(X, A, M, N);
    LUPDecompose(A, N);
    LUPInvert(A, N, IA);
    yTX(y, X, yX, M, N);
    // Calc R2
    numerator = olsR2numerator(yX, IA, N);
    R2 = numerator / yTy(y, M);
    //R2 = yTy(y,M);

    return R2;
}

[Код явно написан в браузере, никогда не компилируется и не тестируется, используйте на свой страх и риск].

Делая это, вы амортизируете стоимость однократного выделения памяти по многим вычислениям, что должно быть намного более эффективным, чем ваш текущий подход.