Vektorilaskenta

Johdanto

• Ohjelmoinnin peruskursseilla lähetään liikkeelle siitä, että dataa käsitellään yksittäisinä data-alkioina eli skalaareina (float, double, int, jne)
• Skalaarimuuttujien väliset alkeisoperaatiot saattavat olla esimerkiksi yksinkertaisia aritmeettisia laskutoimituksia:

int a = 1;
int b = 55;
int c = -3;
int r = a * b + c

• Oleellinen asia on se, että kukin alkeisoperaatio ottaa argumenttinaan joukon skalaarimuuttujia ja paluuarvo on myöskin skalaarimuuttuja

• Sama pätee myöskin vertailuoperaatioiden tapauksessa:

if(r < 34) {
    ...
} else {
    ...
}

• Vertailuoperaation osapuolet ovat skalaareja ja tulos on bool-tyyppinen skalaari eli voitaisiin kirjoittaa:

bool ret = r < 34; // false, koska c = 1 * 55 + (-3) = 52
if(ret) { ...

• Vektorilaskennan ideana on käsitellä dataa skalaarien sijasta vektoreina
• Matemaattisesti ajateltuna kyseessä on suure jolla on sekä suunta, että suuruus (pituus)
• Ohjelmoinnin kannalta vektori voidaan ajatella muuttujaksi, joka sisältää indeksoidun joukon skalaareja
• Vektorit eroavat taulukoista siinä, että taulukkojen tapauksessa alkiot käsitellään tyypillisesti erikseen skalaareina, mutta vektoreiden tapauksessa alkiot on tarkoitus käsitellä yhdessä

• Vektorimuuttuja voisi sisältää esimerkiksi neljä kappaletta yksinkertaisen tarkkuuden liukulukuja:

float4 va = { 4.6f, 3.4f, -6.9f, 1.0f };
float4 vb = { 5.3f, -7.6f, 5.9f, 0.5f };

• Yllä esitelty float4-tietotyyppi ei ole millään tavalla standardi vaan eri vektoriohjelmontirajapinnat määrittelevät tyypillisesti omat tietotyyppinsä ja alkeisoperaationsa

• Useimman vektorilaskentaohjelmointirajapinnat tarjoavat mahdollisuuden käsitellä vektorimuuttujien alkioita yksittäisinä skalaarimuuttujina:

float d1 = va.s0; // = vektorimuuttujan va ensimmäinen alkio
float d2 = va.s1; // = vektorimuuttujan va toinen alkio
float d3 = d1 * d2 + vb.s3;
vb.s1 = d3; // vektorimuuttujan vb toinen alkio = d3

• Vektorimuuttujia ei ole kuitenkaan järkevää käsitellä tällä tavalla vaan meidän tulisi pyrkiä käyttämään vektoroituja alkeisoperaatioita.
• Lähdetään seuraavaksi tarkastelemaan esimerkkien kautta kuinka tämä olisi mahdollista.

Esimerkki

• Oletetaan, että meillä on taulukollinen yksinkertaisen tarkkuuden liukulukuja \[ \left( x_1, x_2, \dots, x_{n-1}, x_n \right) \] ja haluamme muodostaa uuden taulukon \[ \left( f(x_1), f(x_2), \dots, f(x_{n-1}), f(x_n) \right), \] jossa \[ \begin{align} f : \mathbb{R} \to \mathbb{R}^+, \; f(x) = \begin{cases} |x|, & 1 \leq |x| \\ x^2, & |x| < 1. \end{cases} \end{align} \]

Skalaariversio

• C-ohjelmointikielellä taulukon muodostaminen voitaisiin toteuttaa esimerkiksi seuraavasti:

void function(float *buff, int n) {
  
    for(int i = 0; i < n; i++) {
        float val = buff[i];
        float abs_val = fabs(val);
      
        float ret;
        if(abs_val < 1.0f)
            ret = val*val;
        else
            ret = abs_val;
 
        buff[i] = ret;
    }
}

• Useimmissa tilanteissa kääntäjä tuottaa ohjelmakoodista skalaariversion, jossa ehtolause on toteutettu hyppykäskyllä.

• Testiohjelman vektorointi-esimerkki/scalar.cpp tulostus tietokoneella, josta löytyy Intel Core i5 4690K -prosessori:

Time: 0.0583539 s
Flops: 0.666179 GFlops
Max error: 0

Advanced Vector Extensions

• Lähdetään seuraavaksi tarkastelemaan, että miten algoritmi voitaisiin vektoroida manuaalisesti käyttäen moderneista x86-prosessoreista löytyvää AVX-käskykanta-laajennusta (Advanced Vector Extensions) ja siihen liittyvää C-kielen käärettä (wrapper).

Silmukan muokaaminen ja datan hakeminen muistista

• AVX-käskykantalaajennuksen avulla 8 yksinkertaisen tarkkuuden liukulukua voidaan pakata samaan __m256-tyyppiseen vektorimuuttujaan.
• Laajennus tukee myöskin kokonaislukujen ja kaksinkertaisen tarkkuuden liukulukujen käsittelyä.
• Muokkaamme aluksi silmukkaa siten että siinä käsitellään 8 liukulukua yhdessä:

for(int i = 0; i < n; i += 8) {
    ...
}

• Kaikki 8 liukulukua voidaan ladata buff-taulukosta _mm256_load_ps-aliohjelmakutsulla:

 __m256 val = _mm256_load_ps(buff+i); // vmovaps ymm, m256

• Tässä kohtaa oletataan, että n on jaollinen 8:lla ja, että osoitinmuuttujan buff arvo on jaollinen 32:lla sillä _mm256_load_ps-aliohjelma (vmovaps-käsky) heittää laitteistotason poikkeuksen mikäli sille annettu osoite ei jaollinen 32:lla.

Itseisarvon ja neliön laskeminen

• Seuraavaksi haluamme laskea kaikille vektorin val komponenteille itseisarvot. Tätä varten olemme ennen silmukkaa määritelleet seuraavat muuttujat:

const unsigned int abs_coef = 0x7fffffff; 

// vbroadcastss ymm, m32
const __m256 abs_mask = _mm256_broadcast_ss((float*) &abs_coef);

• Muuttuja abs_coef pitää sisällään bittimaskin, jonka avulla yksinkertaisen tarkkuuden liukuluvun merkkibitti voidaan nollata eli luvusta otetaan itseisarvo.
• _mm256_broadcast_ss-aliohjelma palauttaa vektorin (abs_mask), jonka jokaisen komponentin arvoksi on asetettu argumenttina annetusta muistiosoitteesta löytyvä liukuluku.
• Tässä tapauksessa emme kuitenkaan tallenna abs_mask-vektorin komponentteihin oikeaa liukulukua (0x7fffffff = NaN (Not a Number)) vaan käytämme vektorin binääriesitystä bittimaskina.

• Nyt itseisarvot voidaan laskea ottamalla komponenteittainen bittitason ja-operaatio (logical AND) vektoreista val ja abs_mask:

__m256 abs_val = _mm256_and_ps(val, abs_mask); // vandps ymm, ymm, ymm

• Seuraavaksi laskemme kaikille vektorin val komponenteille neliöt:

__m256 squ_val = _mm256_mul_ps(val, val); // vmulps ymm, ymm, ymm

Vertailuoperaatio

• Vertailuoperaatio voidaan toteuttaa vektoroidusti _mm256_cmp_ps-aliohjelmakutsulla.
• Aliohjelma ottaa argumenttinaan vertailtavat vektorit ja vertailuoperaattorin. Aliohjelma palauttaa vektorin, jonka kukin komponentti pitää sisällään sitä vastaavan vertailuoperaation tuloksen eli tässä tapauksessa joko 0x00000000 tai 0xffffffff.

• Kuten itseisarvon laskemisenkin yhteydessä, olemme etukäteen määritelleet seuraavat muuttujat:

const float one_coef = 1.0f;
const __m256 one_mask = _mm256_broadcast_ss(&one_coef);

• Eli vektorin one_mask jokaisesta komponentista löytyy liukuluku 1.0f.

Nyt tarvitsemamme "pienempi kuin" -vertailu onnistuu esimerkiksi _CMP_LT_OS-vertailuoperaattorilla:

// vcmpps ymm, ymm, ymm, imm
__m256 pre = _mm256_cmp_ps(abs_val, one_mask, _CMP_LT_OS);

Lopullisen tuloksen muodostaminen "blend"-komennolla

• Seuraavaksi muodostamme vastausvektorin valitsemalla abs_var ja squ_val vektoreista sopivat komponentit.
• Yhdistäminen tapahtuu _mm256_blendv_ps-aliohjelmalla, joka ottaa argumentteinaan kaksi vektoria ja ohjausvektorin, joka määrittelee millä tavalla nämä kaksi vektoria tulisi yhdistää uudeksi vektoriksi.
• Mikäli ohjausvektorin komponentti on nolla, valitaan uuden vektorin vastaavaksi komponentiksi ensimmäisen argumenttina annetun vektorin vastaava komponentti. Muissa tapauksissa valitaan toisen argumenttina annetun vektorin vastaava komponentti.

• Tässä tapauksessa saamme muodotettua haluavamme vastausvektorin käyttämällä edellä muodostettua pre-vektoria ohjausvektorina:

// vblendvps ymm, ymm, ymm, ymm
_m256 ret = _mm256_blendv_ps(abs_val, squ_val, pre); 

Viimeistely

• Lopuksi tallennamme ret-vektorin takaisin buff taulukkoon:

_mm256_store_ps(buff+i, ret); // vmovaps m256, ymm

Lopullinen aliohjelma

• Tässä vielä koko algoritmi kokonaisuudessaan:

const unsigned int abs_coef = 0x7fffffff;
const float one_coef = 1.0f;

void function(float *buff, int n) {
  
    const __m256 abs_mask = _mm256_broadcast_ss((float*) &abs_coef);
    const __m256 one_mask = _mm256_broadcast_ss(&one_coef);
    
    for(int i = 0; i < n; i += 8) {
        __m256 val = _mm256_load_ps(buff+i);
        __m256 abs_val = _mm256_and_ps(val, abs_mask);
        __m256 squ_val = _mm256_mul_ps(val, val);
        __m256 pre = _mm256_cmp_ps(abs_val, one_mask, _CMP_LT_OS);
        __m256 ret = _mm256_blendv_ps(abs_val, squ_val, pre);
        _mm256_store_ps(buff+i, ret);
    }
}

• Testiohjelman vektorointi-esimerkki/vector.cpp tulostus tietokoneella, josta löytyy Intel Core i5 4690K -prosessori:

Time: 0.012953 s
Flops: 3.00116 GFlops (estimate)
Max error: 0

• Vektoroitu ohjelmakoodi on siis noin 4,5 kertaa nopeampi kuin skalaariversio

Vertailuoperaation hyödyntäminen ehtolauseessa

• Edellisen luvun esimerkissä silmukan sisällä ei esiintynyt yhtään ehtolausetta vaan kävimme läpi molemmat mahdollisuudet (\(|x_i| < 1\) ja \(1 \leq |x_i|\)) ja valitsimme sen jälkeen oikean vaihtoehdon "blend"-komennolla hyödyntäen pre-vertailuvektoria.
• Voisimme kuitenkin lisätä ehtolauseen, jolloin välttyisimme neliön laskemiselta ja "blend"-komennolta niissä tilanteissa, joissa kaikki vektorin komponentit toteuttaisivat ehdon \(1 \leq |x_i|\).
• Tämä onnistuu parhaiten _mm256_movemask_ps-aliohjelmalla, joka palauttaa int-tyyppisen muuttujan, jonka 8 vähiten merkitsevää bittiä määräytyvät argumenttina annetun vektorin komponenttien eniten merkitsevien bittien mukaan.

• Tällöin voisimme antaa vertailuvektorin pre argumenttina kyseiselle aliohjelmalle ja käyttää paluuarvoa ehtolauseessa:

int cond = _mm256_movemask_ps(pre); // vmovmskps r32, ymm
if(cond) {
    __m256 squ_val = _mm256_mul_ps(val, val);
    __m256 ret = _mm256_blendv_ps(abs_val, squ_val, pre);
    _mm256_store_ps(buff+i, ret);
} else {
    _mm256_store_ps(buff+i, abs_val);
}

• Seuraava kuva sisältä kaksi esimerkkitapausta _mm256_movemask_ps-aliohjelma käytöstä:

• Ylemmän kaltaiseen tapaukseen törmätään silloin kun val-vektori sisältää komponentteja, joissa \(|x_i| < 1\) eli osa pre-vertailuvektorin komponenteista on saanut arvon 0xffffffff. Tällöin näitä elementtejä vastaavat bitit cond-muuttujassa saavat arvon 1 ja muuttujan cond totuusarvo vertailulauseessa on true.
• Jälkimmäisessä tapauksessa kaikki vektorin val komponentit toteuttavat ehdon \(1 \leq |x_i|\), joten kaikki pre-vertailuvektorin komponentit ovat nollia ja siten cond-muuttuja saa arvon false.

• Emme ota tässä kohtaa sen tarkemmin kantaa siihen, että onko ehtolauseen arvon palauttamisesta mitään oleellista hyötyä yleisessä tapauksessa.
• Voidaan kuitenkin mainita sen verran, että mikäli yksittäinen vektorin komponentti toteuttaa ehdon \(|x_i| < 1\) todennäköisyydellä \(p\) niin tällöin ainakin yksi vektorin komponenteista toteuttaa ehdon todennäköisyydellä \(1-(1-p)^8\).

• Esimerkiksi tapaus, jossa pelkästään 1% luvuista kuuluu puoliavoimelle välille \((-1,1)\) johtaa tilanteeseen, jossa cond-muuttujalla on lähes 8% todennäköisyys saada nollasta eroava arvo. Yhdeksän prosentin todennäköisyys puolestaan johtaa jo lähes 53% todennäköisyyteen.
• Kannattaa myöskin huomioida, että hyppykäskyillä on yleisesti ottaen negatiivinen vaikutus ohjelman suorituskykyyn sillä modernit prosessorit pyrkivät spekuloimaan ehtolauseiden lopputuloksilla jo etukäteen ja epäonnistunut ennuste johtaa pahimmillaan koko liukuhihnan mitätöimiseen.

• Testiohjelman vektorointi-esimerkki/vector2.cpp tulostus tietokoneella, josta löytyy Intel Core i5 4690K -prosessori:

Time: 0.0164099 s
Flops: 2.36895 GFlops (estimate)
Max error: 0

• Ehtolauseen lisääminen hidasti koodia kuten edellä spekulointiinkin

Vektorin komponenttien permutointi

• Joissain tilanteissa saattaisimme haluta permutoida vektorin komponentit keskenään.
• Lähdetään esimerkin vuoksi tarkastelemaan tapausta, jossa meille on annettu kaksi __m256-tyyppistä vektoria --- a ja b --- ja haluaisimme muodostaa uuden vektorin, joka sisältää ne vektoreiden a ja b komponentit, joiden indeksinumerot ovat parittomia.

• Tätä varten määrittelemme ensiksi permutaatiovektorin, joka siirtää vektorin parittomasti indeksoidut komponentit vektorin alkuun ja parillisesti indeksoidut komponentit vektorin loppuun:

union perm_vector {
    __m256i m;
      __m256 f;
    int i[8];
};

const union perm_vector perm1 = { .i = { 1, 3, 5, 7, 0, 2, 4, 6 } };

• Nyt voimme permutoida vektorin komponentit _mm256_permutevar8x32_ps-aliohjelmalla (vaatii AVX2-käskykantalaajennuksen), joka ottaa argumentteinaan permutoitavan vektorin ja __m256i-tyyppisen (8 etumerkillistä 32-bittistä kokonaislukua) permutaatiovektorin. Paluuarvona saadaan permutoitu vektori.
• Permutaatiovektorin komponentti määrää sen permutoitavan vektorin komponentin, jonka arvo kopioidaan permutoidun vektorin vastaavaan komponenttiin.

Voimme nyt soveltaa perm1-permutaatiovektoria vektoriin a:

 __m256 aa = _mm256_permutevar8x32_ps(a, perm1.m);

• Samalla tavalla voimme määritellä permutaatiovektorin, joka siirtää vektorin parillisesti indeksoidut komponentit vektorin alkuun ja parittomasti indeksoidut komponentit vektorin loppuun.
• Sen jälkeen voimme soveltaa kyseistä permutaatiovektoria vektoriin b:

const union perm_vector perm2 = { .i = { 0, 2, 4, 6, 1, 3, 5, 7 } };
__m256 bb = _mm256_permutevar8x32_ps(b, perm2.m);

• Lopuksi voimme yhdistää permutoidut vektorit aa ja bb siten että lopullisen vektorin neljä ensimmäistä komponenttia valitaan vektorin aa alusta ja neljä jälkimmäistä komponenttia vektorin bb lopusta. Tämä onnistuu tutulla _mm256_blendv_ps-aliohjelmalla:

const union perm_vector blend = { .i = { 0, 0, 0, 0, 1, 1, 1, 1 } };
return _mm256_blendv_ps(aa, bb, blend.f);

• Lopuksi vielä koko ohjelmakoodi aliohjelman muodossa:

#include "immintrin.h"

union perm_vector {
    __m256i m;
      __m256 f;
    int i[8];
};

const union perm_vector perm1 = { .i = { 1, 3, 5, 7, 0, 2, 4, 6 } };
const union perm_vector perm2 = { .i = { 0, 2, 4, 6, 1, 3, 5, 7 } };
const union perm_vector blend = { .i = { 0, 0, 0, 0, 1, 1, 1, 1 } };

__m256 even_elements(__m256 a, __m256 b) {
    __m256 aa = _mm256_permutevar8x32_ps(a, perm1.m);
    __m256 bb = _mm256_permutevar8x32_ps(b, perm2.m);
    return _mm256_blendv_ps(aa, bb, blend.f);
}

Käänteisluvun laskeminen vektoroidusti

• Luentomonisteen alussa esitetty esimerkkikoodi, jossa laskettiin lukujen käänteislukuja, voidaan myöskin vektoroida varsin helposti.
• Tällä kertaa käytämme __m256d-tyyppisiä vektorimuuttujia, joihin voidaan tallentaa neljä kaksinkertaisen tarkkuuden liukulukua.
• Sen lisäksi hyödynnämme _mm256_fmadd_pd-aliohjelmaa, joka ottaa argumenttinaan kolme vektoria ja laskee komponentittaisen fma-operaation (Fused multiply-add) eli ensimmäisen argumentin komponentti kerrotaan toisen argumentin vastaavalla komponentilla ja tulokseen lisätään kolmannen argumentin vastaava komponentti.

#include "immintrin.h"

void process(double *a, double *b, int n) {
  
    // Vektoriarvoiset vakiot alkuarvauksen ja Newtonin iteraation laskemiseen
    const __m256d coefA = { 48.0/17.0, 48.0/17.0, 48.0/17.0, 48.0/17.0 };
    const __m256d coefB = { 32.0/17.0, 32.0/17.0, 32.0/17.0, 32.0/17.0 };
    const __m256d coefC = { 1.0, 1.0, 1.0, 1.0 };

    // Muokattu silmukka
    for(int i = 0; i < n; i += 4) {

        __m256d md = - _mm256_load_pd(b+i);

        // Alkuarvaus, y_j = 32/17 * (md)_j + 48/17
        __m256d y = _mm256_fmadd_pd(coefB, md, coefA); 

        __m256d t;

        // Yksittainen Newtonin iteraatio on jaettu kahteen fma-operaatioon
        t = _mm256_fmadd_pd(md, y, coefC); // (t)_j = (md)_j * (y)_j + 1.0
        y = _mm256_fmadd_pd(y, t, y);      // (y)_j = (y)_j * (t)_j + (y)_j

        t = _mm256_fmadd_pd(md, y, coefC);
        y = _mm256_fmadd_pd(y, t, y);

        t = _mm256_fmadd_pd(md, y, coefC);
        y = _mm256_fmadd_pd(y, t, y);

        t = _mm256_fmadd_pd(md, y, coefC);
        y = _mm256_fmadd_pd(y, t, y);

        _mm256_store_pd(b+i, y);
    }
}