Salut a tous j'ai le bout de code suivant :

template <typename T>
Matrix<T> Matrix<T>::dot512(const Matrix<T> &a)const {
        if(this->m_width != a.m_length)
                throw std::invalid_argument("Bad size for multiplication\n");

        Matrix<T> ret(this->m_length, a.m_width);
        for(size_t i = 0u; i < a.m_width; i++) {
                std::vector<T> colA(this->m_width);
                for(size_t j = 0u; j < this->m_width; j++)
                        colA[j] = a.m_tab[j][i];
                        
                for(size_t j = 0u; j < this->m_length; j++) {
                        int tmp = 0;

                        {
                                __m256 vecA = _mm256_set_ps(m_tab[j][k], m_tab[j][k+1], m_tab[j][k+2], m_tab[j][k+3], m_tab[j][k+4], m_tab[j][k+5], m_tab[j][k+6], m_tab[j][k+7]);
                                __m256 vecB = _mm256_set_ps(colA[k], colA[k+1], colA[k+2], colA[k+3],colA[k+4], colA[k+5], colA[k+6], colA[k+7]);
                                __m256 res  = _mm256_mul_ps(vecA, vecB);
                                tmp += sum8(res);
                        }
                        /*
                        for(size_t k = 0; k < this->m_width; k++) {
                                tmp += this->m_tab[j][k] * colA[k];
                        }*/
                        ret.m_tab[j][i] = tmp;
                }
        }
        return ret;
}

Donc ce code permet de faire une multiplication de matrices assez classique.

Le probleme c'est que si je commente la boucle

for(size_t k = 0; k < this->m_width; k+=8)

 et que je decommente l'autre boucle, mon code s'execute environ 4 fois plus rapidement, avec l'option -O3

Quelqu'un aurait une idée??

J'utilise cette multiplication dans un strassen que j'ai parlleliser. Peut-etre que le les threads doivent attendre pour avoir acces aux registres mmx??

Merci de votre aide.

Bonjour,

Attention tu as coupé la ligne 14 au lieu de la copier.

Il y a des conversions de type dans tes boucles (tmp est int, les SIMD sont float et les éléments des tableaux des T que tu n'a pas précisé). Sachant qu'une conversion float<=>int est plus couteuse qu'une multiplication, la plupart du CPU est dans les conversions et dans la fonction sum8() que je ne connais pas.
Suivant le type T que tu utilises, il y auras des gains très différents. L'optimiseur a certainement mieux géré les conversions et cela ne m'étonne pas qu'il fasse mieux.

Salut,

J'avoue, je ne connais pas __m256 et ce genre de types. Car tout ce qui est spécifique à certains processeurs, je préfère laisser le compilo l'optimiser.

As tu essayé de comparer ta version à un code sans __m256, pour voir si le compilo arrive à faire aussi bien, voir mieux ?

Car parfois, en forçant l'optimisation, on débranche des optimisations automatique qui sont bien plus performantes.

Après il faut connaître quelques règles : les instructions SIMD sont invoquées dans les for quand il n'y a pas de if dedans par exemple. Et c'est le cas dans ton code. 

Il y a un (voire deux) problèmes d’inefficacités au niveau de la gestion de la mémoire dynamique avant d'aller chercher la réécriture avec une utilisation explicite de la vectorisation.

Si vraiment tu veux des perfs, je ne vois aucune raison de ne pas employer une bibliothèque dédiée (eigen/blaze/armadillo) et de plugger le noyau de calculs sur la MKL que l'on peut obtenir légalement et gratuitement aujourd'hui.

Salut,

Merici a vous deux pour vos reponses. Voici le code de la fonction sum8(je l'ai reprise sur stackOverflow).

La fonction permet de faire la somme des termes d'un registre simd.

float sum8(__m256 x) {
    // hiQuad = ( x7, x6, x5, x4 )
    const __m128 hiQuad = _mm256_extractf128_ps(x, 1);
    // loQuad = ( x3, x2, x1, x0 )
    const __m128 loQuad = _mm256_castps256_ps128(x);
    // sumQuad = ( x3 + x7, x2 + x6, x1 + x5, x0 + x4 )
    const __m128 sumQuad = _mm_add_ps(loQuad, hiQuad);
    // loDual = ( -, -, x1 + x5, x0 + x4 )
    const __m128 loDual = sumQuad;
    // hiDual = ( -, -, x3 + x7, x2 + x6 )
    const __m128 hiDual = _mm_movehl_ps(sumQuad, sumQuad);
    // sumDual = ( -, -, x1 + x3 + x5 + x7, x0 + x2 + x4 + x6 )
    const __m128 sumDual = _mm_add_ps(loDual, hiDual);
    // lo = ( -, -, -, x0 + x2 + x4 + x6 )
    const __m128 lo = sumDual;
    // hi = ( -, -, -, x1 + x3 + x5 + x7 )
    const __m128 hi = _mm_shuffle_ps(sumDual, sumDual, 0x1);
    // sum = ( -, -, -, x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 )
    const __m128 sum = _mm_add_ss(lo, hi);
    return _mm_cvtss_f32(sum);
}

Pour l'instant je ne m'occupe que des types de 32 bits.

Aussi, j'ai réecri ma boucle pour ne pas a avoir a faire de conversion d'int vers float :

                        for(size_t k = 0; k < this->m_width; k+=8) {
                                __m256i f = _mm256_set_epi32(this->m_tab[i][j], this->m_tab[i][j+1], this->m_tab[i][j+2], this->m_tab[i][j+3], this->m_tab[i][j+4], this->m_tab[i][j+5], this->m_tab[i][j+6], this->m_tab[i][j+7]);
                                __m256i s = _mm256_set_epi32(a.m_tab[i][j], a.m_tab[i][j+1], a.m_tab[i][j+2], a.m_tab[i][j+3], a.m_tab[i][j+4], a.m_tab[i][j+5], a.m_tab[i][j+6], a.m_tab[i][j+7]);
                                __m256i r = _mm256_add_epi32(f, s);
                                T *res    = (T*) (&r);
                                tmp += res[7];
                                tmp += res[6];
                                tmp += res[5];
                                tmp += res[4];
                                tmp += res[3];
                                tmp += res[2];
                                tmp += res[1];
                                tmp += res[0];
                        }

Maintenant, ma multiplication avec des int est 2 fois plus rapide que celle naive avec les optimisations.

Mais avec des float,  elle est toujours environ 2 fois plus lente.

Merci de vos reponses

lmghs a écrit:

Il y a un (voire deux) problèmes d’inefficacités au niveau de la gestion de la mémoire dynamique avant d'aller chercher la réécriture avec une utilisation explicite de la vectorisation.

Si vraiment tu veux des perfs, je ne vois aucune raison de ne pas employer une bibliothèque dédiée (eigen/blaze/armadillo) et de plugger le noyau de calculs sur la MKL que l'on peut obtenir légalement et gratuitement aujourd'hui.

Salut, tu peux preciser svp que je puisse corriger ca.

-
Edité par Cencoremoi 3 janvier 2021 à 14:45:58

Tu as une allocation+libération en plein milieu d'un boucle.

Et tu vas générer une nouvelle matrice à chaque appel de fonction -- celle là, je n'ai pas vérifié si l'élimination de temporaires est simple à faire comme c'est le cas dès qu'il y a des sommes ou des vecteurs; je ne sais pas comment ils font dans blitz++, eigen, et cie pour les multiplications de matrices.

Je n'ai pas vérifié non plus si tu avais bien interverti les deux boucles les plus profondes, ni si tu travaillais bien sur des sous-blocs.

PS: sur les archi intel, la MKL est ce qu'il y a de plus rapide, et il est pratiquement impossible de faire mieux.

-
Edité par lmghs 3 janvier 2021 à 15:02:34

Tu parles de

                std::vector<T> colA(this->m_width);
                for(size_t j = 0u; j < this->m_width; j++)
                        colA[j] = a.m_tab[j][i];

C'est pour que tout soit aligné dans le cache, j'ai essayer par blocs mais comme je n'ai besoin de le faire que pour des matrices de taille 512 par 512(strassen qui s'arrete quand les matrices font cette taille.), cette methode est plus efficace(a moins que je m'y sois mal pris).

Pour l'instant avec des int j'obtient de tres bonnes perf(7.5 secondes pour des matrices de 8192 par 8192) en parallelisant et avec le SIMD.

Je parle de ça. Ca n'a rien à faire dans une boucle.

for(size_t i = 0u; i < a.m_width; i++) {
                std::vector<T> colA(this->m_width);

Et si tu veux que cela soit aligné, utilises un allocateur qui te garantisse l'alignement (ex. boost::aligned_allocator).

J’utiliserai aussi une lib portable et générique de vectorisation. On commence à en avoir plusieurs: xsimd, VcDevel, etc.

-
Edité par lmghs 3 janvier 2021 à 15:18:48

Ca marche merci de tes reponses, je vais essayer d'aller voir vers une générique de vectorisation comme tu me l'as conseiller.

En revanche pour l'allocation memoire, elle est indispensable, je vais aller chercher plusieurs fois les elements d'une meme colonne dans le cache ce qui va etre beaucoup plus rapide lorsque tout est aligné dans la memoire.

De plus, meme si je n'ai jamias utiliser boost::aligned_allocator, je ne suis pas sur que cela change grand chose si je dois parcourir les elements de la meme colonne rapidement.

"Sors l'alloc de la boucle, et réutilise le vecteur" est mon message

Quant à l'alignement, c'est pour pouvoir directement manipuler des séquences calées sur les tailles des caches sans avoir à avoir une première boucle de peeling qui peut autrement être nécessaire. Cela ne retirera pas un éventuel besoin de boucle finale à la main.

En effet, je n'avais pas vu la déclaration du vector avec un width dans le for.

Réduire absolument le nombre d'allocations gagnera beaucoup. Donc alloue tout ce que tu as à allouer avant les boucles. Mais je dis ce qui a été dit en fait...

Juste pour apporter ma pierre à l'édifice : je bosse dans une boite ou on fait des programmes qui calculent beaucoup, et dégager les allocations, ça optimise vraiment bien les zones critiques, c'est un fait.

-
Edité par Fvirtman 3 janvier 2021 à 21:54:14