када пакетная холеская факторизация

Question

Я вроде понимаю, как бороться с 2D CUDA. Но пакетный холецкий имеет 4D к концу алгоритма. Я приложил cholesky и мой код cuda, если кто-нибудь мог дать мне подсказку.

int i, k, m, n;
    // Batched Cholesky factorization.
    for (i = 0; i < batch; i++) {

            float *pA = &dA[i*N*N];

            // Single Cholesky factorization.
            for (k = 0; k < N; k++) {

                    // Panel factorization.
                    pA[k*N+k] = sqrtf(pA[k*N+k]);

                    for (m = k+1; m < N; m++)
                            pA[k*N+m] /= pA[k*N+k];

                    // Update of the trailing submatrix.
                    for (n = k+1; n < N; n++)
                            for (m = n; m < N; m++)
                                    pA[n*N+m] -= (pA[k*N+n]*pA[k*N+m]);
            }


    }

Cuda:

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int k = blockIdx.y * blockDim.y + threadIdx.y;
    int m = blockIdx.z * blockDim.z + threadIdx.z;
    int n = blockIdx.z * blockDim.z + threadIdx.z;

    if( k >= N || m >= N || n >= N || i >= batch ) return;
    float *pA = &dA[i*N*N];
    pA[k*N+k] = sqrtf(pA[k*N+k]);
    pA[k*N+m] /= pA[k*N+k];
    pA[n*N+m] -= (pA[k*N+n]*pA[k*N+m]);

стартер:

    dim3 dimBlock( (batch+31)/32, (n+31)/32, (n+31)/32 );
    dim3 dimGrid( 32, 32, 32);

    spotrf_batched_kernel<<< dimBlock, dimGrid, 0, stream>>>(n, batch, dA);

Answer 1

Я собираюсь оставить это здесь без особых комментариев. Код относительно понятен. Эта реализация полностью соответствует вашей последовательной версии со следующими функциями:

1. Каждый блок выполняет ровно одну факторизацию в пакете. Запустите столько блоков, сколько есть пакетных матриц для разложения.
2. Поскольку факторизация выполняется в пределах блока, возможна синхронизация между параллельными операциями, поэтому соблюдается порядок операций факторизации
3. Единственный параллелизм, который раскрывает алгоритм, находится внутри строковых операций факторизации и операций обновления
4. Размеры блоков должны соответствовать количеству строк в матричном размере пакета в круглых кратных размеру основы (на сегодняшний день 32 на всех устройствах с поддержкой CUDA)

Код, приведенный ниже, был чрезвычайно легко протестирован и не гарантирует его работоспособность или правильность. Используйте на свой страх и риск:

#include <iostream>
#include <algorithm>

__global__ 
void batchkernel(float** batches, int nbatches, int N, int LDA)
{
    if (blockIdx.x < nbatches) {
        float* pA = batches[blockIdx.x];
        for (int k = 0; k < N; k++) {

            // Panel factorization.
            if (threadIdx.x == 0)  {
                pA[k*LDA+k] = sqrtf(pA[k*LDA+k]);
            }
            __syncthreads();

            for (int m = threadIdx.x; ((m < N) && (threadIdx.x > k)); m+=blockDim.x) {
                pA[k*LDA+m] /= pA[k*LDA+k];
            }
            __syncthreads();

            // Update of the trailing submatrix.
            for (int n = k+1; (n < N); n++) {
                for (int m = threadIdx.x; ((m < N) && (threadIdx.x >= n)); m+=blockDim.x) {
                    pA[n*LDA+m] -= pA[k*LDA+n] * pA[k*LDA+m];
                }
            }
            __syncthreads();
        }
    }
}

void refCholeskey(float* pA, int N)
{
    int k, m, n;

    // Single Cholesky factorization.
    for (k = 0; k < N; k++) {
        // Panel factorization.
        pA[k*N+k] = sqrtf(pA[k*N+k]);

        for (m = k+1; m < N; m++)
            pA[k*N+m] /= pA[k*N+k];

        // Update of the trailing submatrix.
        for (n = k+1; n < N; n++)
            for (m = n; m < N; m++)
                pA[n*N+m] -= (pA[k*N+n]*pA[k*N+m]);
    }
}


int main()
{
    // B = np.random.random((10,10))
    // SPDmatrix = (0.5*(B+B.T)) + B.shape[0]*np.eye(B.shape[0])
    const int N = 10;
    const int LDA = 10;
    float SPDmatrix[LDA*N] = {
    10.22856331,   0.17380577,   0.61779525,   0.66592082,   0.46915566,
     0.09946502,   0.69386511,   0.35224291,   0.53155506,   0.51441469,
     0.17380577,  10.67971161,   0.34481401,   0.64766522,   0.22372943,
     0.55896022,   0.59083588,   0.48872497,   0.54049871,   0.74764959,
     0.61779525,   0.34481401,    10.229388,   0.40904432,    0.5015491,
     0.52152334,   0.19684814,   0.28262256,   0.04384535,   0.61919751,
     0.66592082,   0.64766522,   0.40904432,  10.78410647,   0.12708693,
      0.3241063,    0.6984497,   0.65074097,   0.08027563,   0.56332844,
     0.46915566,   0.22372943,    0.5015491,   0.12708693,  10.52234091,
     0.76346103,   0.80932473,    0.8234331,   0.52737611,   0.65777357,
     0.09946502,   0.55896022,   0.52152334,    0.3241063,   0.76346103,
    10.54906761,   0.32865411,   0.32467483,   0.80720007,   0.36287463,
     0.69386511,   0.59083588,   0.19684814,    0.6984497,   0.80932473,
     0.32865411,  10.29729551,   0.34707933,   0.69379356,   0.87612982,
     0.35224291,   0.48872497,   0.28262256,   0.65074097,    0.8234331,
     0.32467483,   0.34707933,  10.42929929,   0.78849458,     0.159371,  
     0.53155506,   0.54049871,   0.04384535,   0.08027563,   0.52737611,
     0.80720007,   0.69379356,   0.78849458,  10.49604818,   0.43871288,
     0.51441469,   0.74764959,   0.61919751,   0.56332844,   0.65777357,
     0.36287463,   0.87612982,     0.159371,   0.43871288,  10.94535485 };

    const int nbatches = 8;
    float** batches;
    cudaMallocManaged((void **)&batches, nbatches * sizeof(float*));

    for(int i=0; i<nbatches; i++) {
        cudaMallocManaged((void **)&batches[i], N * LDA * sizeof(float));
        cudaMemcpy(batches[i], SPDmatrix, N * LDA * sizeof(float), cudaMemcpyDefault);
    }

    int blocksz = 32;
    int nblocks = nbatches;

    batchkernel<<<nblocks, blocksz>>>(batches, nbatches, N, LDA);
    refCholeskey(SPDmatrix, N);
    cudaDeviceSynchronize();

    float maxabsrelerror = 0.0f;
    for(int i = 0; i < N*N; i++) {
        float absrelerror = std::fabs(SPDmatrix[i] - batches[0][i]) / std::fabs(SPDmatrix[i]);
        maxabsrelerror = std::max(absrelerror, maxabsrelerror);
    }
    std::cout << "Maximum absolute relative error = " << maxabsrelerror << std::endl;

    cudaDeviceReset();

    return 0;
}