CODEWITHSUNDEEP
cudamatrix-multiplication
Baback Kh
Baback Kh

Reputation: 1

Cuda multiplication of an m-by-n matrix in an n-by-1 vector

The following kernel multiplies two n-by-n matrices:

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

I use the above kernel in MATLAB like this:

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

but the result is not zero! How can I change this kernel so that I can multiply an m-by-n matrix in an n-by-1 vector?

Upvotes: 0

Views: 1247

Answers (1)

talonmies
talonmies

Reputation: 72350

Because it transpires you are using Matlab (and arrays are stored in column major order), you will have to make a very minor modification to the kernel code:

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

The kernel requires one thread per entry in the output matrix, so you will need to run 3 x 1 threads in you (3 x 4) * (4 x 1) example. When run you should see this:

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

which is correct for column major order storage.

Upvotes: 1

Related Questions