Algorithme d'optimisation de la multiplication matricielle et vectorielle

Question

Supposons que les dimensions soient très grandes (jusqu'à 1 milliard d'éléments dans une matrice). Comment pourrais-je implémenter un algorithme inconsistant de cache pour le produit de matrice-vecteur? Basé sur wikipedia, je vais devoir diviser et conquérir récursivement, mais je pense qu'il y aurait beaucoup de frais généraux. Serait-il efficace de le faire?

Suivi de questions et réponses: OpenMP with matrices and vectors

Answer 1

La réponse à la question, « comment puis-je faire cette opération d'algèbre linéaire de base rapide », est toujours et partout pour trouver et un lien vers une bibliothèque BLAS à l'écoute pour votre Plate-forme. Par exemple, GotoBLAS (dont le travail se poursuit dans OpenBLAS), ou le plus lent autotuned ATLAS, ou des paquets commerciaux comme Intel MKL. L'algèbre linéaire est si fondamentale pour tant d'autres opérations que d'énormes efforts sont consacrés à l'optimisation de ces paquets pour diverses plates-formes, et il n'y a aucune chance que vous trouviez quelque chose dans un travail de quelques après-midi. Les appels de sous-programme particuliers que vous recherchez pour la multiplicaiton matrice-vecteur dense générale sont SGEMV/DGEMV/CGEMV/ZGEMV.

algorithmes de cache-inconscients, ou autoréglage, sont quand vous ne pouvez pas être pris la peine mise au point pour l'architecture spécifique de cache de votre système - qui pourrait bien, normalement, mais puisque les gens sont prêts à le faire pour BLAS routines, et ensuite rendre les résultats à l'écoute disponibles, signifie que vous êtes le meilleur simplement en utilisant ces routines. Le modèle d'accès mémoire pour GEMV est assez simple pour que vous n'ayez pas vraiment besoin de diviser pour régner (même chose pour le cas standard de la transposition matricielle) - vous trouvez simplement la taille de blocage du cache et l'utilisez. Dans GEMV (y = Axe), vous devez toujours parcourir la matrice entière une seule fois, donc rien à faire pour la réutilisation (et donc l'utilisation efficace du cache), mais vous pouvez essayer de réutiliser x autant que possible pour le charger une fois au lieu de (nombre de lignes) fois - et vous voulez toujours que l'accès à A soit compatible avec le cache. Le cache évident blocage chose à faire est de briser le long des blocs:

A x -> [ A11 | A12 ] | x1 | = | A11 x1 + A12 x2 | 
     [ A21 | A22 ] | x2 | | A21 x1 + A22 x2 | 

Et vous pouvez certainement le faire récursive. Mais faire une mise en œuvre naïve, il est plus lent que le simple double boucle, et le mode plus lent qu'un appel bibliothèque SGEMV appropriée:

$ ./gemv 
Testing for N=4096 
Double Loop: time = 0.024995, error = 0.000000 
Divide and conquer: time = 0.299945, error = 0.000000 
SGEMV: time = 0.013998, error = 0.000000 

Le code suit:

#include <stdio.h> 
#include <stdlib.h> 
#include <sys/time.h> 
#include "mkl.h" 

float **alloc2d(int n, int m) { 
    float *data = malloc(n*m*sizeof(float)); 
    float **array = malloc(n*sizeof(float *)); 
    for (int i=0; i<n; i++) 
     array[i] = &(data[i*m]); 
    return array; 
} 

void tick(struct timeval *t) { 
    gettimeofday(t, NULL); 
} 

/* returns time in seconds from now to time described by t */ 
double tock(struct timeval *t) { 
    struct timeval now; 
    gettimeofday(&now, NULL); 
    return (double)(now.tv_sec - t->tv_sec) + ((double)(now.tv_usec - t->tv_usec)/1000000.); 
} 

float checkans(float *y, int n) { 
    float err = 0.; 
    for (int i=0; i<n; i++) 
     err += (y[i] - 1.*i)*(y[i] - 1.*i); 
    return err; 
} 

/* assume square matrix */ 
void divConquerGEMV(float **a, float *x, float *y, int n, 
        int startr, int endr, int startc, int endc) { 

    int nr = endr - startr + 1; 
    int nc = endc - startc + 1; 

    if (nr == 1 && nc == 1) { 
     y[startc] += a[startr][startc] * x[startr]; 
    } else { 
     int midr = (endr + startr+1)/2; 
     int midc = (endc + startc+1)/2; 
     divConquerGEMV(a, x, y, n, startr, midr-1, startc, midc-1); 
     divConquerGEMV(a, x, y, n, midr, endr, startc, midc-1); 
     divConquerGEMV(a, x, y, n, startr, midr-1, midc, endc); 
     divConquerGEMV(a, x, y, n, midr, endr, midc, endc); 
    } 
} 
int main(int argc, char **argv) { 
    const int n=4096; 
    float **a = alloc2d(n,n); 
    float *x = malloc(n*sizeof(float)); 
    float *y = malloc(n*sizeof(float)); 
    struct timeval clock; 
    double eltime; 

    printf("Testing for N=%d\n", n); 

    for (int i=0; i<n; i++) { 
     x[i] = 1.*i; 
     for (int j=0; j<n; j++) 
      a[i][j] = 0.; 
     a[i][i] = 1.; 
    } 

    /* naive double loop */ 
    tick(&clock); 
    for (int i=0; i<n; i++) { 
     y[i] = 0.; 
     for (int j=0; j<n; j++) { 
      y[i] += a[i][j]*x[j]; 
     } 
    } 
    eltime = tock(&clock); 
    printf("Double Loop: time = %lf, error = %f\n", eltime, checkans(y,n)); 

    for (int i=0; i<n; i++) y[i] = 0.; 

    /* naive divide and conquer */ 
    tick(&clock); 
    divConquerGEMV(a, x, y, n, 0, n-1, 0, n-1); 
    eltime = tock(&clock); 
    printf("Divide and conquer: time = %lf, error = %f\n", eltime, checkans(y,n)); 

    /* decent GEMV implementation */ 
    tick(&clock); 

    float alpha = 1.; 
    float beta = 0.; 
    int incrx=1; 
    int incry=1; 
    char trans='N'; 

    sgemv(&trans,&n,&n,&alpha,&(a[0][0]),&n,x,&incrx,&beta,y,&incry); 
    eltime = tock(&clock); 
    printf("SGEMV: time = %lf, error = %f\n", eltime, checkans(y,n)); 

    return 0; 
}