Умножение Куды матрицы m-на-n в векторе n-на-1

Question

Следующее ядро ​​умножает две n-на-n матрицы:

    __global__ void matrixMultiplication(const double *A, const double *B, double *C, int N)
{
    int i = blockDim.y * blockIdx.y + threadIdx.y;
    int j = blockDim.x * blockIdx.x + threadIdx.x;
    double value = 0;
    for(int k = 0; k < N; k++){
    value += A[k * N + j] * B[i * N + k];
    }
    C[i * N + j] = value;
    }

Я использую вышеуказанное ядро ​​в MATLAB следующим образом:

k = parallel.gpu.CUDAKernel('matrixMultiplication.ptx', 'matrixMultiplication.cu');
A = rand(3,4);
b = rand(4,1);
C = zeros(3,1);
k.ThreadBlockSize = [3,4,1];
k.GridSize = [1, 1];
D = A*b;
C = feval(k,A,b,C,4);
D-C

но результат не ноль! Как я могу изменить это ядро, чтобы умножить матрицу m-на-n в векторе n-на-1?

Answer 1

Поскольку выясняется, что вы используете Matlab (а массивы хранятся в основном порядке столбцов), вам придется внести очень незначительные изменения в код ядра:

#include <thrust/iterator/counting_iterator.h>
#include <thrust/copy.h>
#include <thrust/device_vector.h>
#include <iostream>

__global__ void matrixMultiplication(const double *A, const double *B, double *C, int M, int N)
{
    int i = blockDim.y * blockIdx.y + threadIdx.y;
    int j = blockDim.x * blockIdx.x + threadIdx.x;
    double value = 0;
    for(int k = 0; k < N; k++){
        value += A[k * M + j] * B[i * M + k];
    }
    C[i * N + j] = value;
}

int main()
{
    const int M = 3, N = 4, K = 1;
    thrust::device_vector<double> A(M*N), B(N*K), C(M*K);

    thrust::counting_iterator<double> counter(1.0);
    thrust::copy(counter, counter + (M*N), A.begin());
    thrust::copy(counter, counter + (N*K), B.begin());

    dim3 grid(1,1), block(M,K);

    matrixMultiplication<<<grid, block>>>( thrust::raw_pointer_cast(A.data()),
                                           thrust::raw_pointer_cast(B.data()),
                                           thrust::raw_pointer_cast(C.data()),
                                           M, N );

    cudaDeviceSynchronize();

    for(int i=0; i<M*K; i++)
        std::cout << C[i] << std::endl;

    return 0;
}

Ядру требуется по одному потоку на каждую запись в выходной матрице, поэтому вам нужно запустить 3 x 1 потока в вашем примере (3 x 4) * (4 x 1). При запуске вы должны увидеть это:

$ nvcc -arch=sm_52 -std=c++11 -o spoonfull spoonfull.cu 
$ ./spoonfull 
70
80
90

, что верно для хранения основного заказа столбца.