Une manière plus rapide de faire l'addition de matrice multidimensionnelle?

Question

J'ai une matrice A de taille (m * l * 4) et la taille de m est d'environ 100 000 et l = 100. la taille de la liste est toujours égale à n et n < = m. Je voulais faire une addition matricielle d'une liste donnée d'index. J'ai écrit cette fonction et dois appeler cette fonction un grand nombre de fois.

void MatrixAddition(int l, int n, vector<int>& list, int ***A,int ***C,int cluster) 
{ 
    for(int i=0;i<l;i++) 
    { 
     for(int j=0;j<4;j++) 
       C[cluster][i][j]=0; 
    } 

for (int i = 0; i < l; i++) 
{ 
     for(int j=0;j<n;++j) 
    { 
     for(int k=0;k<4;k++) 
      C[cluster][i][k]+=A[list[j]][i][k]; 
    } 
} 

} 

J'utilise gprof pour calculer combien de temps prend chaque morceau de fonction dans le code entier et je trouve mon 60% du temps pris par fonction addition matricielle. Existe-t-il un autre moyen d'écrire cette fonction afin que mon temps d'exécution diminue.

secondes de temps secondes appels ms/ms appel/appel Nom
52,00 7,85 7,85 20 392,60 405,49 addition matricielle (int, int, std :: vector> &, int ***, int ***, int)

Answer 1

Échangez la boucle par i et la boucle par j dans la deuxième partie. Cela rendra la fonction plus facile à mettre en cache.

for(int j=0;j<n;++j) 
{ 
    for (int i = 0; i < l; i++) 
    { 
     for(int k=0;k<4;k++) 
      C[cluster][i][k]+=A[list[j]][i][k]; 
    } 
} 

En outre, j'espère que vous n'avez pas oublié le drapeau -O3.

Answer 2

(mise à jour: une version antérieure avait une mauvaise indexation Cette version s'auto-vectorise assez facilement).

Utilisez un tableau multidimensionnel C (plutôt qu'un tableau de pointeurs vers des pointeurs) ou un tableau 1D plat indexé avec i*cols + j, afin que la mémoire soit contiguë. Cela fera une énorme différence dans l'efficacité de la prélecture de matériel pour faire bon usage de la bande passante mémoire. Les charges avec une adresse provenant d'une autre charge aspirent vraiment à la performance, ou inversement, avoir des adresses prévisibles bien connues à l'avance aide beaucoup car les charges peuvent démarrer bien avant qu'elles ne soient nécessaires (grâce à une exécution dans le désordre).

En outre, la réponse de @ user31264 est correcte, vous devez échanger les boucles afin que la boucle j soit la plus extérieure. C'est bien mais loin d'être suffisant en soi.

Ceci permettra également au compilateur de s'auto-vectoriser. En fait, j'ai eu du mal à faire gcc auto-vectoriser. (Mais c'est probablement parce que je suis l'indexation mal, parce que je ne regardais le code la première fois. Donc, le compilateur ne savait pas que nous bouclez mémoire contiguë.)

---

J'ai joué avec elle sur le Godbolt compiler explorer.

J'ai finalement obtenu une bonne belle sortie du compilateur de cette version, qui prend A et C sous forme de tableaux 1D plat et fait l'indexation elle-même:

void MatrixAddition_contiguous(int rows, int n, const vector<int>& list, 
           const int *__restrict__ A, int *__restrict__ C, int cluster) 
    // still auto-vectorizes without __restrict__, but checks for overlap and 
    // runs a scalar loop in that case 
{ 
    const int cols = 4; // or global constexpr or something 
    int *__restrict__ Ccluster = C + ((long)cluster) * rows * cols; 

    for(int i=0;i<rows;i++) 
    //#pragma omp simd 
    for(int k=0;k<4;k++) 
     Ccluster[cols*i + k]=0; 

    for(int j=0;j < cols;++j) { // loop over clusters in A in the outer-most loop 
    const int *__restrict__ Alistj = A + ((long)list[j]) * rows * cols; 
    // #pragma omp simd // Doesn't work: only auto-vectorizes with -O3 
    // probably only -O3 lets gcc see through the k=0..3 loop and treat it like one big loop 
    for (int i = 0; i < rows; i++) { 
     long row_offset = cols*i; 
     //#pragma omp simd // forces vectorization with 16B vectors, so it hurts AVX2 
     for(int k=0;k<4;k++) 
     Ccluster[row_offset + k] += Alistj[row_offset + k]; 
    } 
    } 
} 

levage manuellement le list[j] certainement aidé le compilateur se rendre compte que les magasins en C ne peut pas affecter les indices qui seront chargés à partir de list[j]. Il n'était probablement pas nécessaire de soulever manuellement les autres objets.

Le levage A[list[j]], plutôt que seulement list[j], est un artefact d'un previous approach where I had the indexing wrong. Tant que nous levons la charge de list[j] autant que possible, le compilateur peut faire un bon travail même s'il ne sait pas que list ne chevauche pas C.

La boucle intérieure, avec gcc 5.3 ciblant x86-64 -O3 -Wall -march=haswell -fopenmp (et -fverbose-asm) est:

.L26: 
    vmovdqu ymm0, YMMWORD PTR [r8+rax]  # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B] 
    vpaddd ymm0, ymm0, YMMWORD PTR [rdx+rax] # vect__71.75, MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: vectp.73_90, index: ivtmp.91_26, offset: 0B] 
    add  r12d, 1 # ivtmp.88, 
    vmovdqu YMMWORD PTR [r8+rax], ymm0  # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], vect__71.75 
    add  rax, 32 # ivtmp.91, 
    cmp  r12d, r9d # ivtmp.88, bnd.66 
    jb  .L26  #, 

Il est donc fait huit ajoute à la fois, avec AVX2 vpaddd, avec des charges non alignés et un magasin non alignés en arrière -vectorizing, il devrait faire du bon code avec ARM NEON, ou PPC Altivec, ou tout ce qui peut faire l'addition 32bit emballés.

Je ne pouvais pas obtenir gcc pour me dire quoi que ce soit avec -ftree-vectorizer-verbose=2, mais -Rpass-analysis=loop-vectorize de clang était légèrement utile.