le meilleur pari sera d'utiliser des instructions SSE ps qui fonctionnent sur 4 flotteurs à la fois. Vous pouvez profiter du fait qu'un float 0.0 est tout 0 bits pour utiliser une instruction andps pour masquer les éléments indésirables:

#include <stdint.h>
#include <xmmintrin.h>

union {
    uint32_t i[4];
    __m128   xmm;
} mask[16] = {
 {  0,  0,  0,  0 },
 { ~0,  0,  0,  0 },
 {  0, ~0,  0,  0 },
 { ~0, ~0,  0,  0 },
 {  0,  0, ~0,  0 },
 { ~0,  0, ~0,  0 },
 {  0, ~0, ~0,  0 },
 { ~0, ~0, ~0,  0 },
 {  0,  0,  0, ~0 },
 { ~0,  0,  0, ~0 },
 {  0, ~0,  0, ~0 },
 { ~0, ~0,  0, ~0 },
 {  0,  0, ~0, ~0 },
 { ~0,  0, ~0, ~0 },
 {  0, ~0, ~0, ~0 },
 { ~0, ~0, ~0, ~0 },
};

float dot(__m128 *a, uint64_t b) {
    __m128 sum = { 0.0 };
    for (int i = 0; i < 16; i++, b>>=4)
        sum += _mm_and_ps(a[i], mask[b&0xf].xmm);
    return sum[0] + sum[1] + sum[2] + sum[3];
}

Si vous vous attendez à être beaucoup de 0 dans le masque, il pourrait être plus rapide à court cicruit les 0:

for (int i = 0; b; i++, b >>= 4)
    if (b & 0xf)
        sum += _mm_and_ps(a[i], mask[b&0xf].xmm);

mais si b est aléatoire, ce sera plus lent.

pour développer la réponse D'Aki Suihkonen, le remodelage de la corde est utile pour le déplacement conditionnel des flotteurs. Dans la solution ci-dessous, une permutation de bits en deux étapes à l'aide des instructions SSE PMOVMASKB et PSHUFB, plus L'instruction BLENDVPS a été utilisée pour obtenir 1,25 éléments manipulés/cycle sur un noyau 2 Duo 2,26 GHz, soit 20 fois la vitesse de mon code C de référence.

[EDIT: Un AVX2 mise en œuvre a été ajouté. La Performance est inconnue parce que je ne peux pas la tester moi-même, mais est devrait être le double de la vitesse. ]

Voici mon implémentation et Test bench, l'explication suit.

SSE4.1 (ancien, pour AVX2 voir ci-dessous)

Code

/* Includes */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <smmintrin.h>  /* SSE 4.1 */
#include <time.h>



/* Defines */
#define ALIGNTO(n) __attribute__((aligned(n)))
#define USE_PINSRW 1
#define NUM_ITERS  2260000



/**
 * Bit mask shuffle.
 * 
 * This version uses a loop to store eight u16 and reloads them as one __m128i.
 */

__m128 bitMaskShuffleStoreAndReload(__m128i mask){
    const __m128i perm    ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4,
                                                     11, 3, 10, 2, 9,  1, 8,  0);
    int i;
    uint16_t interMask[8] ALIGNTO(16);

    /* Shuffle bitmask */
        /* Stage 1 */
    for(i=7;i>=0;i--){
        interMask[i] = _mm_movemask_epi8(mask);
        mask = _mm_slli_epi32(mask, 1);
    }

        /* Stage 2 */
    return _mm_castsi128_ps(
                    _mm_shuffle_epi8(
                        _mm_load_si128((const __m128i*)interMask),
                        perm)
                );
}

/**
 * Bit mask shuffle.
 * 
 * This version uses the PINSTRW instruction.
 */

__m128 bitMaskShufflePINSRW(__m128i mask){
    const __m128i perm    ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4,
                                                     11, 3, 10, 2, 9,  1, 8,  0);
    __m128i  imask ALIGNTO(16);

    /* Shuffle bitmask */
        /* Stage 1 */
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 7);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 6);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 5);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 4);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 3);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 2);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 1);
    mask = _mm_slli_epi16(mask, 1);
    imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 0);

        /* Stage 2 */
    return _mm_castsi128_ps(
                    _mm_shuffle_epi8(
                        imask,
                        perm)
                );
}

/**
 * SSE 4.1 implementation.
 */

float dotSSE41(__m128 f[32], unsigned char maskArg[16]){
    int i, j, k;
    __m128i  mask      ALIGNTO(16) = _mm_load_si128((const __m128i*)maskArg);
    __m128   shufdMask ALIGNTO(16);
    __m128   zblended  ALIGNTO(16);
    __m128   sums      ALIGNTO(16) = _mm_setzero_ps();
    float    sumsf[4]  ALIGNTO(16);

    /* Shuffle bitmask */
#if USE_PINSRW
    shufdMask = bitMaskShufflePINSRW(mask);
#else
    shufdMask = bitMaskShuffleStoreAndReload(mask);
#endif

    /* Dot product */
    for(i=1;i>=0;i--){
        for(j=1;j>=0;j--){
            for(k=7;k>=0;k--){
                zblended  = _mm_setzero_ps();
                zblended  = _mm_blendv_ps(zblended, f[i*16+j+k*2], shufdMask);
                sums      = _mm_add_ps(sums, zblended);
                shufdMask = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(shufdMask), 1));
            }
        }
    }

    /* Final Summation */
    _mm_store_ps(sumsf, sums);
    return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3];
}

/**
 * Reference C implementation
 */


float dotRefC(float f[128], unsigned char mask[16]){
    float sum = 0.0;
    int i;

    for(i=0;i<128;i++){
        sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0;
    }

    return sum;
}

/**
 * Main
 */

int main(void){
    /* Variables */
        /* Loop Counter */
    int           i;

        /* Data to process */
    float         data[128] ALIGNTO(16);
    unsigned char mask[16]  ALIGNTO(16);
    float         refCSum, sseSum;

        /* Time tracking */
    clock_t       t1, t2, t3;
    double        refCTime, sseTime;



    /* Initialize mask and float arrays with some random data. */
    for(i=0;i<128;i++){
        if(i<16)
            mask[i]=rand();

        data[i] = random();
    }


    /* RUN TESTS */
    t1 = clock();
    for(i=0;i<NUM_ITERS;i++){
        refCSum = dotRefC(data, mask);
    }
    t2 = clock();
    for(i=0;i<NUM_ITERS;i++){
        sseSum  = dotSSE41((__m128*)data, mask);
    }
    t3 = clock();


    /* Compute time taken */
    refCTime = (double)(t2-t1)/CLOCKS_PER_SEC;
    sseTime  = (double)(t3-t2)/CLOCKS_PER_SEC;


    /* Print out results */
    printf("Results:\n"
           "RefC: Time: %f    Value: %f\n"
           "SSE:  Time: %f    Value: %f\n",
           refCTime, refCSum,
           sseTime,  sseSum);

    return 0;
}

Explication

BLENDVPS utilise le bit supérieur dans les quatre voies 32 bits du registre xmm0 128 bits pour déterminer s'il faut déplacer ou non la valeur dans la voie correspondante de son opérande source dans son opérande destination. Lors du chargement des données avec MOVAPS, on obtient 4 flotteurs consécutifs: par exemple, les 8ème, 9ème, 10ème et 11ème flotteurs. Bien sûr, leur sélection ou désélection doit être contrôlée par l'ensemble correspondant de bits: Par exemple, les 8e, 9e, 10e et 11e bits dans la chaîne de bits.

le problème est que lorsque le masque est chargé pour la première fois, les morceaux de ces ensembles sont juste à côté les uns des autres (dans les positions 8, 9, 10 et 11), alors qu'en fait ils devraient être à 32 bits l'un de l'autre; rappelez-vous, à un moment donné ils devront occuper la position du 31ème bit de chaque voie (les 31ème, 63ème, 95ème et 127ème positions dans le registre XMM0).

ce qui se produirait idéalement est une transposition un peu qui apporte des bits 0, 4, 8, 12,... sur la voie zéro, bits 1, 5, 9, 13, ... dans la voie 1, bits 2, 6, 10, 14,... ligne 2 et bits 3, 7, 11, 15,... dans le couloir de trois. Ainsi, tous les ensembles de 4 bits qui étaient auparavant contigus sont maintenant rayés à 32 bits l'un de l'autre, un dans chacune des quatre voies de 32 bits. Alors tout ce qu'il faut est un boucle qui itère 32 fois, chaque fois en déplaçant dans la position du bit supérieur de chaque voie un nouvel ensemble de 4 bits.

malheureusement x86 n'est pas doté de bonnes instructions de manipulation de bits, donc, faute d'une manière propre de faire une transposition parfaite, un compromis raisonnable est celui-ci.

dans le masque, les 128 bits

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
 16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31
 32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47
 48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63
 64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79
 80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
 96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

sont permutés, par huit PMOVMASKB et huit instructions PSLLW, d'abord à

  0   8  16  24  32  40  48  56  64  72  80  88  96 104 112 120
  1   9  17  25  33  41  49  57  65  73  81  89  97 105 113 121
  2  10  18  26  34  42  50  58  66  74  82  90  98 106 114 122
  3  11  19  27  35  43  51  59  67  75  83  91  99 107 115 123
  4  12  20  28  36  44  52  60  68  76  84  92 100 108 116 124
  5  13  21  29  37  45  53  61  69  77  85  93 101 109 117 125
  6  14  22  30  38  46  54  62  70  78  86  94 102 110 118 126
  7  15  23  31  39  47  55  63  71  79  87  95 103 111 119 127

et par un seul RSSFB instruction à

  0   8  16  24  32  40  48  56   4  12  20  28  36  44  52  60
 64  72  80  88  96 104 112 120  68  76  84  92 100 108 116 124
  1   9  17  25  33  41  49  57   5  13  21  29  37  45  53  61
 65  73  81  89  97 105 113 121  69  77  85  93 101 109 117 125
  2  10  18  26  34  42  50  58   6  14  22  30  38  46  54  62
 66  74  82  90  98 106 114 122  70  78  86  94 102 110 118 126
  3  11  19  27  35  43  51  59   7  15  23  31  39  47  55  63
 67  75  83  91  99 107 115 123  71  79  87  95 103 111 119 127

. Nous itérons maintenant sur quatre "runs", dont chacun contient huit ensembles de quatre bits étalés à des intervalles de 32 bits (comme nous l'avons souhaité), en utilisant ces ensembles comme la commande de masque pour BLENDVPS. La maladresse inhérente du brassage de bits conduit à la boucle maladroite à triple emboîtement dans dotSSE41(), mais

clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -msse4.1 -mssse3 dot.c -o dottest

les boucles sont de toute façon déroulées. Les itérations de la boucle interne se composent de 16 répétitions

blendvps    0x90(%rsi),%xmm1
addps       %xmm4,%xmm1
pslld       x1,%xmm2
movdqa      %xmm2,%xmm0
xorps       %xmm4,%xmm4

.

à part cela, je n'ai pas été en mesure de préciser laquelle de mes deux versions de bit shuffle était la plus rapide, donc j'ai donné les deux implémentations dans ma réponse.

AVX2 (de nouveau, mais non testé)

Code

/* Includes */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>  /* AVX2 */
#include <time.h>



/* Defines */
#define ALIGNTO(n) __attribute__((aligned(n)))
#define NUM_ITERS  2260000


/**
 * Bit mask shuffle.
 * 
 * This version uses the PINSTRW instruction.
 */

__m256 bitMaskShufflePINSRW(__m256i mask){
    __m256i  imask ALIGNTO(32);

    /* Shuffle bitmask */
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 7);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 6);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 5);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 4);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 3);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 2);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 1);
    mask = _mm256_slli_epi32(mask, 1);
    imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 0);

    /* Return bitmask */
    return _mm256_castsi256_ps(imask);
}

/**
 * AVX2 implementation.
 */

float dotAVX2(__m256 f[16], unsigned char maskArg[16]){
    int i, j, k;
    /* Use _mm_loadu_si128 */
    __m256i  mask      ALIGNTO(32) = _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg));
    __m256   shufdMask ALIGNTO(32);
    __m256   zblended  ALIGNTO(32);
    __m256   sums      ALIGNTO(32) = _mm256_setzero_ps();
    float    sumsf[8]  ALIGNTO(32);

    /* Shuffle bitmask */
    shufdMask = bitMaskShufflePINSRW(mask);
    shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 16));

    /* Dot product */
    for(i=15;i>=0;i--){
        zblended  = _mm256_setzero_ps();
        /* Replace f[i] with _mm256_loadu_ps((float*)&f[i]) */
        zblended  = _mm256_blendv_ps(zblended, f[i], shufdMask);
        sums      = _mm256_add_ps(sums, zblended);
        shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 1));
    }

    /* Final Summation */
    _mm256_store_ps(sumsf, sums);
    return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3] + sumsf[4] + sumsf[5] + sumsf[6] + sumsf[7];
}

/**
 * Reference C implementation
 */

float dotRefC(float f[128], unsigned char mask[16]){
    float sum = 0.0;
    int i;

    for(i=0;i<128;i++){
        sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0;
    }

    return sum;
}

/**
 * Main
 */

int main(void){
    /* Variables */
        /* Loop Counter */
    int           i;

        /* Data to process */
    float         data[128] ALIGNTO(32);
    unsigned char mask[16]  ALIGNTO(32);
    float         refCSum, sseSum;

        /* Time tracking */
    clock_t       t1, t2, t3;
    double        refCTime, sseTime;



    /* Initialize mask and float arrays with some random data. */
    for(i=0;i<128;i++){
        if(i<16)
            mask[i]=rand();

        data[i] = random();
    }


    /* RUN TESTS */
    t1 = clock();
    for(i=0;i<NUM_ITERS;i++){
        refCSum = dotRefC(data, mask);
    }
    t2 = clock();
    for(i=0;i<NUM_ITERS;i++){
        sseSum  = dotAVX2((__m256*)data, mask);
    }
    t3 = clock();


    /* Compute time taken */
    refCTime = (double)(t2-t1)/CLOCKS_PER_SEC;
    sseTime  = (double)(t3-t2)/CLOCKS_PER_SEC;


    /* Print out results */
    printf("Results:\n"
           "RefC: Time: %f    Value: %f\n"
           "SSE:  Time: %f    Value: %f\n",
           refCTime, refCSum,
           sseTime,  sseSum);

    return 0;
}

Explication

le même concept que pour SSE4.1 est utilisé. La différence est que maintenant nous traitons 8 flotteurs à la fois et en utilisant les registres 256 bits D'AVX2 et PMOVMASKB de ymm registres (qui rassemblent 256/8 = 32 bits). Pour cette raison, nous avons maintenant un mélange plus simple, et une boucle Beaucoup plus simple.

dans le masque, les 256 bits

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
 16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31
 32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47
 48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63
 64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79
 80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
 96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255

sont permutés, en utilisant 8 PMOVMASKB et 8 instructions PSLLW, à

  0   8  16  24  32  40  48  56  64  72  80  88  96 104 112 120
128 136 144 152 160 168 176 184 192 200 208 216 224 232 240 248
  1   9  17  25  33  41  49  57  65  73  81  89  97 105 113 121
129 137 145 153 161 169 177 185 193 201 209 217 225 233 241 249
  2  10  18  26  34  42  50  58  66  74  82  90  98 106 114 122
130 138 146 154 162 170 178 186 194 202 210 218 226 234 242 250
  3  11  19  27  35  43  51  59  67  75  83  91  99 107 115 123
131 139 147 155 163 171 179 187 195 203 211 219 227 235 243 251
  4  12  20  28  36  44  52  60  68  76  84  92 100 108 116 124
132 140 148 156 164 172 180 188 196 204 212 220 228 236 244 252
  5  13  21  29  37  45  53  61  69  77  85  93 101 109 117 125
133 141 149 157 165 173 181 189 197 205 213 221 229 237 245 253
  6  14  22  30  38  46  54  62  70  78  86  94 102 110 118 126
134 142 150 158 166 174 182 190 198 206 214 222 230 238 246 254
  7  15  23  31  39  47  55  63  71  79  87  95 103 111 119 127
135 143 151 159 167 175 183 191 199 207 215 223 231 239 247 255

. 128 élément bits-avec-float dot-produits, nous avons ensuite itérer en parallèle sur huit ensembles de 16 éléments. Cette mise en œuvre peut facilement être étendue à 256 DPs en itérant sur des ensembles de 32 éléments. Un seul boucle est nécessaire maintenant.

en particulier, pour changer cela pour travailler pour 256-elment dot produits, vous

• Double la taille des arguments de fonction. __m256 f[32], unsigned char maskArg[32].
• changez la charge du masque (= _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg));) avec = _mm256_load_si256((const __m256i*)maskArg);.
• supprimer le décalage compensateur laissé par 16 juste en dessous de l'appel à bitMaskShufflePINSRW.
• Exécuter la boucle vers le bas à partir de la configuration 31 au lieu de 15: for(i=31;i>=0;i--)

je ne peux ni essai ni même lancer le code comme mon CPU est SSE4.1, mais Clang avec

clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -mavx2 -msse4.1 -mssse3 dotavx2.c -o dottest

compilé proprement( vous pouvez obtenir du code plus rapide sans le dérouler), produisant ceci:

(gdb) disas dotAVX2
Dump of assembler code for function dotAVX2:
0x0000000100001730 <+0>:    push   %rbp
0x0000000100001731 <+1>:    mov    %rsp,%rbp
0x0000000100001734 <+4>:    vmovdqa (%rsi),%xmm0
0x0000000100001738 <+8>:    vpslld x1,%ymm0,%ymm1
0x000000010000173d <+13>:   vpslld x1,%ymm1,%ymm2
0x0000000100001742 <+18>:   vpmovmskb %ymm2,%eax
0x0000000100001746 <+22>:   vpslld x1,%ymm2,%ymm2
0x000000010000174b <+27>:   vpmovmskb %ymm2,%ecx
0x000000010000174f <+31>:   vxorps %ymm3,%ymm3,%ymm3
0x0000000100001753 <+35>:   vmovd  %ecx,%xmm4
0x0000000100001757 <+39>:   vpinsrd x1,%eax,%xmm4,%xmm4
0x000000010000175d <+45>:   vpmovmskb %ymm1,%eax
0x0000000100001761 <+49>:   vpinsrd x2,%eax,%xmm4,%xmm1
0x0000000100001767 <+55>:   vpslld x1,%ymm2,%ymm2
0x000000010000176c <+60>:   vpslld x1,%ymm2,%ymm4
0x0000000100001771 <+65>:   vpslld x1,%ymm4,%ymm5
0x0000000100001776 <+70>:   vpmovmskb %ymm0,%eax
0x000000010000177a <+74>:   vpinsrd x3,%eax,%xmm1,%xmm0
0x0000000100001780 <+80>:   vpmovmskb %ymm5,%eax
0x0000000100001784 <+84>:   vpslld x1,%ymm5,%ymm1
0x0000000100001789 <+89>:   vpmovmskb %ymm1,%ecx
0x000000010000178d <+93>:   vmovd  %ecx,%xmm1
0x0000000100001791 <+97>:   vpinsrd x1,%eax,%xmm1,%xmm1
0x0000000100001797 <+103>:  vpmovmskb %ymm4,%eax
0x000000010000179b <+107>:  vpinsrd x2,%eax,%xmm1,%xmm1
0x00000001000017a1 <+113>:  vpmovmskb %ymm2,%eax
0x00000001000017a5 <+117>:  vpinsrd x3,%eax,%xmm1,%xmm1
0x00000001000017ab <+123>:  vinserti128 x1,%xmm0,%ymm1,%ymm0
0x00000001000017b1 <+129>:  vpslld x10,%ymm0,%ymm0
0x00000001000017b6 <+134>:  vblendvps %ymm0,0x1e0(%rdi),%ymm3,%ymm1
0x00000001000017c0 <+144>:  vpslld x1,%ymm0,%ymm0
0x00000001000017c5 <+149>:  vblendvps %ymm0,0x1c0(%rdi),%ymm3,%ymm2
0x00000001000017cf <+159>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000017d3 <+163>:  vpslld x1,%ymm0,%ymm0
0x00000001000017d8 <+168>:  vblendvps %ymm0,0x1a0(%rdi),%ymm3,%ymm2
0x00000001000017e2 <+178>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000017e6 <+182>:  vpslld x1,%ymm0,%ymm0
0x00000001000017eb <+187>:  vblendvps %ymm0,0x180(%rdi),%ymm3,%ymm2
0x00000001000017f5 <+197>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000017f9 <+201>:  vpslld x1,%ymm0,%ymm0
0x00000001000017fe <+206>:  vblendvps %ymm0,0x160(%rdi),%ymm3,%ymm2
0x0000000100001808 <+216>:  vaddps %ymm2,%ymm1,%ymm1
0x000000010000180c <+220>:  vpslld x1,%ymm0,%ymm0
0x0000000100001811 <+225>:  vblendvps %ymm0,0x140(%rdi),%ymm3,%ymm2
0x000000010000181b <+235>:  vaddps %ymm2,%ymm1,%ymm1
0x000000010000181f <+239>:  vpslld x1,%ymm0,%ymm0
0x0000000100001824 <+244>:  vblendvps %ymm0,0x120(%rdi),%ymm3,%ymm2
0x000000010000182e <+254>:  vaddps %ymm2,%ymm1,%ymm1
0x0000000100001832 <+258>:  vpslld x1,%ymm0,%ymm0
0x0000000100001837 <+263>:  vblendvps %ymm0,0x100(%rdi),%ymm3,%ymm2
0x0000000100001841 <+273>:  vaddps %ymm2,%ymm1,%ymm1
0x0000000100001845 <+277>:  vpslld x1,%ymm0,%ymm0
0x000000010000184a <+282>:  vblendvps %ymm0,0xe0(%rdi),%ymm3,%ymm2
0x0000000100001854 <+292>:  vaddps %ymm2,%ymm1,%ymm1
0x0000000100001858 <+296>:  vpslld x1,%ymm0,%ymm0
0x000000010000185d <+301>:  vblendvps %ymm0,0xc0(%rdi),%ymm3,%ymm2
0x0000000100001867 <+311>:  vaddps %ymm2,%ymm1,%ymm1
0x000000010000186b <+315>:  vpslld x1,%ymm0,%ymm0
0x0000000100001870 <+320>:  vblendvps %ymm0,0xa0(%rdi),%ymm3,%ymm2
0x000000010000187a <+330>:  vaddps %ymm2,%ymm1,%ymm1
0x000000010000187e <+334>:  vpslld x1,%ymm0,%ymm0
0x0000000100001883 <+339>:  vblendvps %ymm0,0x80(%rdi),%ymm3,%ymm2
0x000000010000188d <+349>:  vaddps %ymm2,%ymm1,%ymm1
0x0000000100001891 <+353>:  vpslld x1,%ymm0,%ymm0
0x0000000100001896 <+358>:  vblendvps %ymm0,0x60(%rdi),%ymm3,%ymm2
0x000000010000189d <+365>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000018a1 <+369>:  vpslld x1,%ymm0,%ymm0
0x00000001000018a6 <+374>:  vblendvps %ymm0,0x40(%rdi),%ymm3,%ymm2
0x00000001000018ad <+381>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000018b1 <+385>:  vpslld x1,%ymm0,%ymm0
0x00000001000018b6 <+390>:  vblendvps %ymm0,0x20(%rdi),%ymm3,%ymm2
0x00000001000018bd <+397>:  vaddps %ymm2,%ymm1,%ymm1
0x00000001000018c1 <+401>:  vpslld x1,%ymm0,%ymm0
0x00000001000018c6 <+406>:  vblendvps %ymm0,(%rdi),%ymm3,%ymm0
0x00000001000018cc <+412>:  vaddps %ymm0,%ymm1,%ymm0
0x00000001000018d0 <+416>:  vpshufd x1,%xmm0,%xmm1
0x00000001000018d5 <+421>:  vaddss %xmm1,%xmm0,%xmm1
0x00000001000018d9 <+425>:  vmovhlps %xmm0,%xmm0,%xmm2
0x00000001000018dd <+429>:  vaddss %xmm1,%xmm2,%xmm1
0x00000001000018e1 <+433>:  vpshufd x3,%xmm0,%xmm2
0x00000001000018e6 <+438>:  vaddss %xmm1,%xmm2,%xmm1
0x00000001000018ea <+442>:  vextracti128 x1,%ymm0,%xmm0
0x00000001000018f0 <+448>:  vaddss %xmm1,%xmm0,%xmm1
0x00000001000018f4 <+452>:  vpshufd x1,%xmm0,%xmm2
0x00000001000018f9 <+457>:  vaddss %xmm1,%xmm2,%xmm1
0x00000001000018fd <+461>:  vpshufd x3,%xmm0,%xmm2
0x0000000100001902 <+466>:  vmovhlps %xmm0,%xmm0,%xmm0
0x0000000100001906 <+470>:  vaddss %xmm1,%xmm0,%xmm0
0x000000010000190a <+474>:  vaddss %xmm0,%xmm2,%xmm0
0x000000010000190e <+478>:  pop    %rbp
0x000000010000190f <+479>:  vzeroupper 
0x0000000100001912 <+482>:  retq   
End of assembler dump.

comme nous pouvons le voir, le noyau contient 3 instructions (vblendvps, vaddps, vpslld) maintenant.

si on permet une permutation légèrement différente dans données float data[128], ou la permutation correspondante dans le bitmask pour le __m128 mask;, on peut améliorer légèrement l'algorithme suggéré par Chris Dodd ci-dessus. (Sans compter le temps requis pour permuter le masque, cette implémentation (+ overhead) est environ 25% plus rapide). Bien sûr, ce n'est qu'une ébauche de mon idée fournie dans les commentaires.

union {
    unsigned int i[4];
    float f[4];
    __m128   xmm;
} mask = { 0xFF00FF00, 0xF0F0F0F0, 0xCCCCCCCC, 0xAAAAAAAA };

float dot2(__m128 *a, __m128 mask);
// 20M times 1.161s

float dotref(__m128 *a, unsigned int *mask) // 20M times 8.174s
{
    float z=0.0f;
    int i;
    for (i=0;i<32;i++) {
       if (mask[0] & (0x80000000U >> i)) z+= a[i][0];
       if (mask[1] & (0x80000000U >> i)) z+= a[i][1];
       if (mask[2] & (0x80000000U >> i)) z+= a[i][2];
       if (mask[3] & (0x80000000U >> i)) z+= a[i][3];       
    }
   return z;
}

la mise en oeuvre de l'assembleur correspondant serait:

dot2:
    // warm up stage: fill in initial data and
    // set up registers
    pxor %xmm1, %xmm1      ;; // clear partial sum1
    pxor %xmm2, %xmm2      ;; // clear partial sum2
    movaps (%rdi), %xmm3   ;; // register warm up stage1
    movaps 16(%rdi), %xmm4 ;; // next 4 values
    pxor %xmm5, %xmm5
    pxor %xmm6, %xmm6
    lea 32(%rdi), %rdi
    movl , %ecx            ;; // process 2x4 items per iteration (total=128)
a:  ;; // inner loop -- 2 independent data paths
    blendvps %xmm3, %xmm5
    pslld , %xmm0
    movaps (%rdi), %xmm3   
    blendvps %xmm4, %xmm6
    pslld , %xmm0
    movaps 16(%rdi), %xmm4
    addps %xmm5, %xmm1
    pxor  %xmm5, %xmm5
    addps %xmm6, %xmm2
    pxor  %xmm6, %xmm6
    lea 32(%rdi), %rdi
    loop a
 ;; // cool down stage: gather results (xmm0 = xmm1+xmm2)
 ;; // in beautiful world this stage is interleaved
 ;; // with the warm up stage of the next block
    addps %xmm2, %xmm1
    movaps  %xmm1, %xmm2
    movaps  %xmm1, %xmm0
    shufps  , %xmm1, %xmm2
    addss   %xmm2, %xmm0
    movaps  %xmm1, %xmm2
    unpckhps %xmm1, %xmm2
    shufps  5, %xmm1, %xmm1
    addss   %xmm2, %xmm0
    addss   %xmm1, %xmm0
    ret

Voici quelques choses à essayer.

essayez de faire en sorte que le compilateur utilise CMOV au lieu d'une branche. ( notez que l'utilisation d'une union de cette façon est bien définie en C11 mais non définie en C++11.)

union {
    int i;
    float f;
} u;
u.i = 0;
if (b & 1) {
    u.f = a[i];
}
sum += u.f;

Utiliser une multiplication au lieu d'une branche.

sum += (b & 1) * a[i];

garder plusieurs sommes et les ajouter à la fin pour réduire les dépendances de flux de données. (Vous pouvez combiner l'une des suggestions ci-dessus avec celui-ci.)

float sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
for (int i = 0; i < 64; i += 4; b >>= 4) {
    if (b & 1) sum0 += a[i];
    if (b & 2) sum1 += a[i+1];
    if (b & 4) sum2 += a[i+2];
    if (b & 8) sum3 += a[i+3];
}
return sum0 + sum1 + sum2 + sum3;

Réduire le nombre de branches en traitant plusieurs bits à la fois:

for (int i = 0; i < 64; i += 4, b >>= 4) {
    switch (b & 0xf) {
        case 0:
            break;
        case 1:
            sum += a[i];
            break;
        case 2:
            sum += a[i + 1];
            break;
        case 3:
            sum += a[i] + a[i+1];
            break;
        case 4:
            sum += a[i+2];
            break;
        // etc. for cases up to and including 15
    }
}

Vous pouvez garder plusieurs sommes et pour chaque somme traiter plusieurs bits à la fois. Dans ce cas, vous voudrez probablement utiliser une macro ou une fonction inline et à l'appeler quatre fois.

j'ai trouvé que l'assemblage généré pour le code de Chris Dodd est très fortement dépendant du compilateur; clang le transforme en boucle pendant que icc (12.x et 13.x) déroulez la boucle à la place. Néanmoins, on peut réduire les dépendances (le besoin d'attendre lesum +=) en le transformant en carte-réduire,

float dot(__m128 *a, uint64_t b) {
    __m128 sum[8];
    int i;

    for (i = 0; i < 8; i++) {
            sum[i] = _mm_add_ps(
                    _mm_and_ps(a[2*i], mask[b & 0xf].xmm),
                    _mm_and_ps(a[2*i+1], mask[(b & 0xf0) >> 4].xmm));
            b >>= 8;
    }
    for (i = 0; i < 4; i++) {
            sum[i] = _mm_add_ps(sum[2*i], sum[2*i+1]);
    }
    sum[0] = _mm_add_ps(sum[0], sum[1]);
    sum[2] = _mm_add_ps(sum[2], sum[3]);

    sum[0] = _mm_add_ps(sum[0], sum[2]);

    sum[0] = _mm_hadd_ps(sum[0], sum[0]);
    sum[0] = _mm_hadd_ps(sum[0], sum[0]);

    i = _mm_extract_ps(sum[0], 0);
    return *((float*)(&i));
}

Cela crée distinctement inférieur assemblée clang (qui stocke les sum[] sur la pile) mais mieux code (pas de dépendances sur d'autres addps) avec gcc et icc. Fait intéressant, seulement gcc a l'idée à la fin que le bas floatsum[0] peut être retourné en place...

bel exercice sur la façon de modifier pour les compilateurs spécifiques ...

Si vous avez i7, alors vous avez SSE4.1 et vous pouvez utiliser le _mm_dp_ps intrinsèque, qui fait un produit de point. Avec votre exemple de code, il ressemblerait à

#include <stdint.h>
#include <immintrin.h>

const float fltmask[][4] =
  {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {0, 0, 1, 1},
   {0, 1, 0, 0}, {0, 1, 0, 1}, {0, 1, 1, 0}, {0, 1, 1, 1},
   {1, 0, 0, 0}, {1, 0, 0, 1}, {1, 0, 1, 0}, {1, 0, 1, 1},
   {1, 1, 0, 0}, {1, 1, 0, 1}, {1, 1, 1, 0}, {1, 1, 1, 1}};

// a has 64 elements. b is a bitvector of 64 dimensions.
float dot(float * restrict a, uint64_t b) {
     int i;
    float sum = 0;
    for(i=0; b && i<64; i+=4,b>>=4) {
        __m128 t0 = _mm_load_ps (a);
        a += 4;
        __m128 t1 = _mm_load_ps (fltmask[b & 15]);
        sum += _mm_cvtss_f32 (_mm_dp_ps (t0, t1, 15));
    }
    return sum;
}

PS. Les tableaux a et fltmask mieux vaut être aligné sur 16 octets!

PPS. Lorsqu'il est compilé avec gcc -std=c99 -msse4 -O2 la boucle ressemble à:

.L3:
    movq    %rdx, %rax
    movaps  (%rcx), %xmm1
    shrq    , %rdx
    andl    , %eax
    addq    , %rcx
    addl    , %r8d
    salq    , %rax
    testq   %rdx, %rdx
    dpps    , (%r9,%rax), %xmm1
    addss   %xmm1, %xmm0
    jne .L13

et -O3 son déroulé, bien sûr.

Vous pouvez éliminer la branche comme ceci:

for(int i=0; b && i<64; i++, b>>=1)
    sum += a[i] * (b & 1);

bien que cela ajoutera un mult supplémentaire, au moins il ne sera pas baiser votre pipeline.

une autre façon de contrôler la branche est de l'utiliser comme vous l'êtes, mais aussi d'utiliser une macro de compilateur. Je suppose qu'en gcc c'est le likely(if ...) macro. Vous utiliserez la branche, mais de cette façon vous direz au compilateur que la branche sera exécutée plus que souvent, et gcc en optimisera davantage.

une Autre optimisation qui peut être fait est la "mise en cache" le produit scalaire. Alors au lieu d'avoir une fonction pour calculer le produit scalaire, vous auriez une variable contenant le produit initialisé à 0 au début. Et chaque fois que vous insérez/supprimez/mettez à jour un élément du vecteur, alors vous mettez également à jour la variable contenant le résultat.