Mozenje matrica 4x4 (SIMD, C)

[bmaxa]

Mnogo puta ste se susreli sa mnozenjem matrica, ali da li ste ikad razmislili
da to ubrzate sa SIMD instrukcijama?
SIMD je prilicno aktraktivan za ovu operaciju, zato sto moze da pomnozi i sabere
od 4 pa cak i do 16 brojeva paralelno (avx512).
Usredsredicu se na mnozenje matrica 4x4 zato sto je to cest slucaj u programiranju
grafike, gde se float (32 bitni) koristi za boje pixela.
Znaci nasa operacija je ovo:

void cmul4x4(float res[4][4],volatile float a[4][4],volatile float b[4][4] ){
    for(int i = 0;i<4;++i){
        for(int j= 0;j<4;++j) {
            float sum = 0.0;
            for(int k=0;k<4;++k) {
                sum += a[i][k]*b[k][j];
            }
            res[i][j]=sum;
        }
    }
}

Dakle svaki element rezultata se dobija tako sto se izvede dot proizvod vrste i kolone.
Kako sad to ubrzati sa SIMD?
Pa, evo ovako: treba postaviti u jedan simd registar vrednost prvog elementa iz vrste
pa onda sa njim pomnoziti sve elemente vrste druge matrice.
Tako imamo 4 rezultata odjednom. Onda prelazimo na sledeci element u vrsti prve
a sledecu vrstu druge. I tako ponavljamo sve do kraja te sabiramo rezultate.
I sve to ponovimo za sve vrste prve.
Znaci umesto da radimo kako je u gore navedenom algoritmu, ovde zapravo idemo
u paralelno mnozenje i sabiranje.
To izgleda ovako:

void mul4x4(float res[4][4],volatile float a[4][4],volatile float b[4][4] ) {
    for (int i=0;i<4;++i) {
        __m128 sum = _mm_set1_ps(0.0);
        for(int j=0;j<4;++j){
            __m128 tmp = _mm_set1_ps(a[i][j]);
            __m128 prod = _mm_mul_ps(tmp,*(__m128*)b[j]);
            sum = _mm_add_ps(sum,prod);
        }
        _mm_storeu_ps(res[i],sum);
    }
}

Dakle, prakticno, citavu jednu petlju stedimo
A sad kompletan kod i merenje vremena:

#include <immintrin.h>
#include <stdio.h>
#include <time.h>
const int NITER = 10000000;
void mul4x4(float res[4][4],volatile float a[4][4],volatile float b[4][4] ) {
    for (int i=0;i<4;++i) {
        __m128 sum = _mm_set1_ps(0.0);
        for(int j=0;j<4;++j){
            __m128 tmp = _mm_set1_ps(a[i][j]);
            __m128 prod = _mm_mul_ps(tmp,*(__m128*)b[j]);
            sum = _mm_add_ps(sum,prod);
        }
        _mm_storeu_ps(res[i],sum);
    }
}
void cmul4x4(float res[4][4],volatile float a[4][4],volatile float b[4][4] ){
    for(int i = 0;i<4;++i){
        for(int j= 0;j<4;++j) {
            float sum = 0.0;
            for(int k=0;k<4;++k) {
                sum += a[i][k]*b[k][j];
            }
            res[i][j]=sum;
        }
    }
}
void print4x4(float a[4][4]) {
    for (int i=0;i<4;++i){
        for (int j=0;j<4;++j) {
            printf(" %f",a[i][j]);
        }
        puts("");
    }
}

int main(void) {
    float res[4][4];
    volatile float a[4][4]={{1,1,1,1},{2,2,2,2},{3,3,3,3},{4,4,4,4}},
                   b[4][4]={{16,15,14,13},{12,11,10,9},{8,7,6,5},{4,3,2,1}};
    clock_t start = clock();
    for(int i = 0;i<NITER;++i)
        mul4x4 (res,a,b);
    clock_t end = clock();
    print4x4(res);
    printf("%f\n",(end-start)/(double)CLOCKS_PER_SEC);
    start = clock();
    for(int i = 0;i<NITER;++i)
        cmul4x4 (res,a,b);
    end = clock();
    print4x4(res);
    printf("%f\n",(end-start)/(double)CLOCKS_PER_SEC);
}

kod mene rezultat izgleda ovako:

 40.000000 36.000000 32.000000 28.000000
 80.000000 72.000000 64.000000 56.000000
 120.000000 108.000000 96.000000 84.000000
 160.000000 144.000000 128.000000 112.000000
0.032028
 40.000000 36.000000 32.000000 28.000000
 80.000000 72.000000 64.000000 56.000000
 120.000000 108.000000 96.000000 84.000000
 160.000000 144.000000 128.000000 112.000000
0.374586

Volatile naravno ne treba, ali u ovom slucaju je koristan da se spreci optimizacija
petlji, tj da se kompajler natera da zaista prodje NITER puta petlju
Dakle, SIMD kod mene 10 puta ubrzava racunicu!
Dalje, nakon sto smo ovo uradili sledeci korak bi bio da se ubrza mnozenje matrica
mxn na ovaj nacin
To ostavljam vama u razmisljanje