Comment puis-je ignorer le quatrième élément d'un float4 lorsque j'utilise cublas sgemv?

Question

Une partie du code sur lequel je travaille nécessite d'effectuer une multiplication vectorielle aussi rapide que possible, par exemple en utilisant une bibliothèque tierce optimisée comme cublas (bien que le même principe s'applique à n'importe quel cpu blas).

Le problème est qu'il y a une sorte de foulée entre les éléments dans le vecteur comme ceci:

La matrice est stockée sous forme de tableau 1D 3Nx3N de flotteurs.

Le vecteur est stocké en tant que tableau N 1D de float4s, mais seuls les trois premiers éléments de chaque float4 doivent être utilisés, le quatrième doit être ignoré.

N est de l'ordre de plusieurs millions.

Si le vecteur étaient stockés comme float3 au lieu de float4 je pouvais jeter le pointeur à flotter, comme dans cet exemple de travail:

//Compile with nvcc test.cu -O3 -lcublas -o test 

/* 
Multiply a 3Nx3N float matrix, M, by a vector, X, of N float3 elements 

The result, Y, is a 3N float vector 
----------------------- 

What if X is a vector of N float4? 

How can I tell cublas to skip the forth element? 

*/ 

#include<iostream> 
#include<thrust/device_vector.h> 
#include<cuda_runtime.h> 
#include<cublas_v2.h> 

using namespace std; 

int main(){ 

    int N = 3; 

    thrust::device_vector<float3> X(N); 

    thrust::device_vector<float> Y(3*N); 

    for(int i=0; i<N; i++) 
    X[i] = make_float3(1,1,1); //make_float4(1,1,1,0); //in the case of float4 i.e., The result should be the same 

    thrust::device_vector<float> M(3*N*3*N, 1); 


    cublasHandle_t handle; 

    cublasCreate(&handle); 

    float beta = 0.0f; 
    float alpha = 1.0f; 
    cublasSgemv(handle, CUBLAS_OP_T, 
      3*N, 3*N, 
      &alpha, 
      thrust::raw_pointer_cast(&M[0]), 3*N, 
      (float*) thrust::raw_pointer_cast(&X[0]), 1, 
      &beta, 
      thrust::raw_pointer_cast(&Y[0]), 1); 

    cout<<"Performed Y = M·X\n\tX = "; 
    for(int i=0; i<N; i++){ 
    float3 Xi = X[i]; 
    cout<<Xi.x<<" "<<Xi.y<<" "<<Xi.z<<" "; 
    } 
    cout<<"\n\tY = "; 
    for(int i=0; i<3*N; i++){ 
    cout<<Y[i]<<" "; 
    } 
    cout<<endl; 

    return 0; 
} 

Mais, comment puis-je effectuer cette opération si le vecteur X est stocké comme float4 s? Etant donné que float4 * peut être interprété comme un float * avec 4 fois plus d'éléments, la question pourrait être plus générale (bien que je ne m'intéresse qu'au cas float4); S'il y a une foulée entre chaque 3 éléments "utiles". Je veux dire à cublas que le tableau n'est pas coalescent en mémoire. Mais quelque chose comme: Il y a 3 éléments au début, les trois suivants sont des éléments de "stride" après cela, etc. Similaire à ce que vous pouvez faire dans OpenGL avec des objets de tableau de sommets.

EDIT:

Les réponses ont suggéré que la méthode la plus viable est de copier tout le réseau strided dans un réseau temporel, transformé, float3 que cublas comprenne.

Les deux options au moment de le faire sont:

1. Use cudaMemcpy2D 
2. Use a thrust transformation 
3. Use a custom copy kernel 

J'ai écrit ce code pour tester les trois cas:

//Compile with Compile with: nvcc test.cu -O3 -lcublas -o test 
#include<iostream> 
#include<thrust/device_vector.h> 
#include<cuda.h> 
#include<cuda_runtime.h> 
#include<cublas_v2.h> 

using namespace std; 


struct Timer{ 
    cudaEvent_t start, stop; 
    float time; 

    void tic(){ 
    cudaEventCreate(&start); 
    cudaEventCreate(&stop); 
    cudaEventRecord(start, 0); 
    } 
    float toc(){ 
    cudaEventRecord(stop, 0); 
    cudaEventSynchronize(stop); 
    cudaEventElapsedTime(&time, start, stop); 

    cudaEventDestroy(start); 
    cudaEventDestroy(stop); 

    return time; 
    } 

}; 



struct copy_functor{ 
    copy_functor(){} 
    __device__ float3 operator() (const float4& X4){ 
    return make_float3(X4.x, X4.y, X4.z); 
    } 
}; 


__global__ void copy_kernel(const float4* __restrict__ X4, float3* __restrict__ X3, int N){ 
    int id = blockIdx.x*blockDim.x + threadIdx.x; 
    if(id < N){ 
    float4 x4 = X4[id]; 
    X3[id] = make_float3(x4.x, x4.y, x4.z); 
    } 
} 

int main(){ 

    int N = 1000000; 
    int Ntest = 1000; 

    Timer t; 

    thrust::device_vector<float3> X3(N, make_float3(0,0,0)); 
    thrust::device_vector<float4> X4(N, make_float4(1,1,1,10)); 


    /*************************CUDAMEMCPY2D*******************/ 
    t.tic(); 

    for(int i= 0; i<Ntest; i++){ 
    cudaMemcpy2DAsync(thrust::raw_pointer_cast(&X3[0]), 
       3*sizeof(float), 
       thrust::raw_pointer_cast(&X4[0]), 
       4*sizeof(float), 
       3*sizeof(float), 
       N, 
       cudaMemcpyDeviceToDevice); 
    cudaDeviceSynchronize(); 
    } 
    printf ("Time for cudaMemcpy2DAsync: %f ms\n", t.toc()/(float)Ntest); 


    /************************THRUST***********************/ 
    t.tic(); 

    for(int i= 0; i<Ntest; i++){ 
    transform(X4.begin(), X4.end(), X3.begin(), copy_functor()); 
    cudaDeviceSynchronize(); 
    } 

    printf ("Time for thrust transformation: %f ms\n", t.toc()/(float)Ntest); 

    /*********************COPY KERNEL*****************************/ 

    t.tic(); 
    for(int i= 0; i<Ntest; i++){ 
    copy_kernel<<< N/128 + 1, 128 >>>(thrust::raw_pointer_cast(&X4[0]), 
         thrust::raw_pointer_cast(&X3[0]), N); 
    cudaDeviceSynchronize(); 
    } 
    printf ("Time for copy kernel: %f ms\n", t.toc()/(float)Ntest); 


return 0; 
} 

Notez que je joue la moyenne de 1000 exemplaires.

La sortie de ce code dans une GTX 980 est la suivante:

Time for cudaMemcpy2DAsync: 1.465522 ms 
Time for thrust transformation: 0.178745 ms 
Time for copy kernel: 0.168507 ms 

cudaMemcpy2D est un ordre de grandeur plus lent que le reste.

poussée et le noyau copie sont très similaires et la plus rapide

Ce comportement semble rester avec un certain nombre d'éléments.

EDIT2:

D'autres réponses suggèrent que GEMM pourrait être utilisé pour communiquer la foulée. Sans la nécessité d'un tableau temporel.

Interprétation du vecteur matriciel mul. comme Matrix Matrix mul. serait fait comme ceci:

cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_T, 
       3*N, 1 /*m*/, 3*N, 
       &alpha, 
       thrust::raw_pointer_cast(&M[0]), 3*N, 
       (float*) thrust::raw_pointer_cast(&X3[0]), 1 /*ldb*/, 
       &beta, 
       thrust::raw_pointer_cast(&Y[0]), 3*N); 

Cependant, à ce stade, je ne sais pas comment passer X4 au lieu de X3. La solution semble être dans les paramètres m et ldb.

Answer 1

Vous pouvez traiter votre 1-D vecteur float4 comme NX3 matrice de flotteur 2-D avec une foulée de la ligne de 4, et utiliser cudaMemcpy2DAsync pour changer la foulée de 4 à 3 avec

cudaMemcpy2DAsync(dst, 
        3*sizeof(float), 
        src, 
        4*sizeof(float), 
        3*sizeof(float), 
        N, 
        cudaMemcpyDeviceToDevice); 

Puis le dst peut être traité comme un vecteur flottant 1N 1N et transmis directement à gemv().

Compte tenu de l'échelle de votre N, l'heure de la copie n'est pas perceptible par rapport à gemv().

---

EDIT

résultat de référence de @Apo montre qu'il est préférable d'utiliser un noyau de copie au lieu de cudaMemcpy2DAsync. J'étais sur-attendu sur cudaMemcpy2DAsync et je pensais que ce serait bien optimisé et avoir les meilleures performances pour tous les cas.