Accelerare hardware: același algoritm, de zece ori mai repede

#ifndef COMUN_H
#define COMUN_H

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

/* Imagine de test: 1920 x 1080, trei octeti pe pixel. */
#define LATIME   1920
#define INALTIME 1080
#define PIXELI   ((size_t)(LATIME) * (INALTIME))

/* Ceas monoton: nu sare inapoi cand se sincronizeaza ora sistemului. */
static inline double acum(void) {
    struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return t.tv_sec + t.tv_nsec * 1e-9;
}

/* Imagine determinista: aceleasi date la fiecare rulare, deci comparabile. */
static void umple(uint8_t* rgb, size_t pixeli) {
    for (size_t i = 0; i < pixeli; i++) {
        rgb[i * 3 + 0] = (uint8_t)((i * 7) & 0xFF);
        rgb[i * 3 + 1] = (uint8_t)((i * 13) & 0xFF);
        rgb[i * 3 + 2] = (uint8_t)((i * 29) & 0xFF);
    }
}

/* Verificarea corectitudinii: o varianta rapida dar gresita nu valoreaza nimic. */
static int compara(const uint8_t* a, const uint8_t* b, size_t n, const char* nume) {
    for (size_t i = 0; i < n; i++) {
        if (a[i] != b[i]) {
            printf("  EROARE in %s la pixelul %zu: %u != %u\n", nume, i, a[i], b[i]);
            return 0;
        }
    }
    return 1;
}

/* Masurare cu incalzire si repetari: intoarce timpul MEDIAN, nu media.
   Medianul nu este influentat de o rulare izolata perturbata de sistem. */
static int cmp_double(const void* x, const void* y) {
    double a = *(const double*)x, b = *(const double*)y;
    return (a > b) - (a < b);
}

#define REPETARI 11

static double masoara(void (*functie)(const uint8_t*, uint8_t*, size_t),
                      const uint8_t* intrare, uint8_t* iesire, size_t pixeli) {
    /* Incalzire: aduce datele in memoria cache si stabilizeaza frecventa. */
    for (int i = 0; i < 3; i++) functie(intrare, iesire, pixeli);

    double timpi[REPETARI];
    for (int i = 0; i < REPETARI; i++) {
        double t0 = acum();
        functie(intrare, iesire, pixeli);
        timpi[i] = acum() - t0;
    }
    qsort(timpi, REPETARI, sizeof(double), cmp_double);
    return timpi[REPETARI / 2];
}

#endif

#include "comun.h"

void gri_scalar(const uint8_t* rgb, uint8_t* gri, size_t pixeli) {
    for (size_t i = 0; i < pixeli; i++) {
        uint32_t r = rgb[i * 3 + 0];
        uint32_t g = rgb[i * 3 + 1];
        uint32_t b = rgb[i * 3 + 2];
        gri[i] = (uint8_t)((77 * r + 150 * g + 29 * b) >> 8);
    }
}

#include "comun.h"
#include <arm_neon.h>       /* functiile intrinseci NEON */

void gri_neon(const uint8_t* rgb, uint8_t* gri, size_t pixeli) {
    /* Ponderile, replicate in toate cele 8 benzi ale registrului. */
    const uint8x8_t p_r = vdup_n_u8(77);
    const uint8x8_t p_g = vdup_n_u8(150);
    const uint8x8_t p_b = vdup_n_u8(29);

    size_t i = 0;
    for (; i + 8 <= pixeli; i += 8) {
        /* vld3_u8 incarca 24 de octeti si ii DESPARTE automat pe cele trei
           canale: 8 valori R, 8 valori G, 8 valori B. Exact ce ne trebuie
           pentru date intretesute - o instructiune pe care un compilator
           rareori o alege singur. */
        uint8x8x3_t px = vld3_u8(rgb + i * 3);

        /* Inmultire cu largire: 8 biti x 8 biti -> 16 biti, ca sa nu depaseasca.
           77 * 255 = 19635, deci 8 biti nu ar fi ajuns. */
        uint16x8_t suma = vmull_u8(px.val[0], p_r);

        /* Inmultire cu largire SI acumulare, intr-o singura instructiune. */
        suma = vmlal_u8(suma, px.val[1], p_g);
        suma = vmlal_u8(suma, px.val[2], p_b);

        /* Deplasare la dreapta cu 8 SI ingustare inapoi la 8 biti,
           tot printr-o singura instructiune. */
        vst1_u8(gri + i, vshrn_n_u16(suma, 8));
    }

    /* Coada: pixelii ramasi, cand numarul lor nu e multiplu de 8.
       Aceasta parte se uita foarte usor si produce o dunga gresita in imagine. */
    for (; i < pixeli; i++) {
        uint32_t r = rgb[i * 3 + 0], g = rgb[i * 3 + 1], b = rgb[i * 3 + 2];
        gri[i] = (uint8_t)((77 * r + 150 * g + 29 * b) >> 8);
    }
}

#include "comun.h"
#include <omp.h>

void gri_neon(const uint8_t* rgb, uint8_t* gri, size_t pixeli);   /* varianta 3 */

void gri_paralel(const uint8_t* rgb, uint8_t* gri, size_t pixeli) {
    int fire = omp_get_max_threads();

    /* Impartim in bucati multiplu de 8, ca fiecare fir sa poata folosi
       calea NEON completa, fara coada proprie. */
    size_t bucata = ((pixeli / fire) / 8) * 8;

    #pragma omp parallel for schedule(static)
    for (int f = 0; f < fire; f++) {
        size_t start = f * bucata;
        size_t cate = (f == fire - 1) ? (pixeli - start) : bucata;
        gri_neon(rgb + start * 3, gri + start, cate);
    }
}

#include "comun.h"

void gri_scalar(const uint8_t*, uint8_t*, size_t);
void gri_neon(const uint8_t*, uint8_t*, size_t);
void gri_paralel(const uint8_t*, uint8_t*, size_t);

typedef struct {
    const char* nume;
    void (*functie)(const uint8_t*, uint8_t*, size_t);
} Varianta;

int main(void) {
    uint8_t* rgb = malloc(PIXELI * 3);
    uint8_t* iesire = malloc(PIXELI);
    uint8_t* referinta = malloc(PIXELI);
    if (!rgb || !iesire || !referinta) return 1;

    umple(rgb, PIXELI);

    /* Referinta de corectitudine, calculata o singura data. */
    gri_scalar(rgb, referinta, PIXELI);

    Varianta variante[] = {
        { "scalar",        gri_scalar  },
        { "NEON",          gri_neon    },
        { "NEON + 4 fire", gri_paralel },
    };

    printf("imagine %d x %d = %.1f milioane de pixeli\n\n",
           LATIME, INALTIME, PIXELI / 1e6);
    printf("%-16s %10s %12s %10s %8s\n",
           "varianta", "timp", "Mpixeli/s", "accelerare", "corect");
    printf("%s\n", "------------------------------------------------------------");

    double t_referinta = 0;
    for (size_t v = 0; v < sizeof(variante) / sizeof(variante[0]); v++) {
        memset(iesire, 0, PIXELI);
        double t = masoara(variante[v].functie, rgb, iesire, PIXELI);
        if (v == 0) t_referinta = t;

        int corect = compara(iesire, referinta, PIXELI, variante[v].nume);

        printf("%-16s %8.2f ms %12.1f %9.2fx %8s\n",
               variante[v].nume, t * 1000.0, PIXELI / t / 1e6,
               t_referinta / t, corect ? "da" : "NU");
    }

    /* Cate cadre pe secunda inseamna asta, in practica? */
    printf("\nla 30 de cadre pe secunda, bugetul este de %.1f ms pe cadru\n",
           1000.0 / 30.0);

    free(rgb); free(iesire); free(referinta);
    return 0;
}