Optimierungsalgorithmus für Matrix- und Vektormultiplikation

Question

Angenommen, die Dimensionen sind sehr groß (bis zu 1 Milliarde Elemente in einer Matrix). Wie würde ich einen Cache-Observious-Algorithmus für Matrix-Vektor-Produkte implementieren? Basierend auf wikipedia muss ich rekursiv teilen und erobern, aber ich habe das Gefühl, dass es eine Menge Overhead geben würde. Wäre es effizient, dies zu tun?

Follow-up-Frage und Antwort: OpenMP with matrices and vectors

Answer 1

Also die Antwort auf die Frage: „Wie kann ich diesen grundlegenden lineare Algebra Betrieb mache schnell“, ist immer und überall zu finden und ein Link zu einer abgestimmten BLAS-Bibliothek für Ihre Plattform. ZB GotoBLAS (deren Arbeit wird in OpenBLAS fortgesetzt), oder die langsamer Autotuned ATLAS, oder kommerzielle Pakete wie Intels MKL. Lineare Algebra ist für so viele andere Operationen so grundlegend, dass enorme Anstrengungen unternommen werden, um diese Pakete für verschiedene Plattformen zu optimieren, und es gibt einfach keine Chance, dass Sie sich in ein paar Nachmittagsarbeiten etwas einfallen lassen, das konkurrieren wird. Die bestimmten Subroutinenaufrufe, nach denen Sie nach einer allgemeinen dichten Matrix-Vektor-Multiplikation suchen, sind SGEMV/DGEMV/CGEMV/ZGEMV.

Cache-vergesslich Algorithmen oder Auto-Tuning, sind für, wenn Sie nicht-Tuning für die spezifische Cache-Architektur des Systems gestört werden - was in Ordnung sein könnte, in der Regel, aber da die Menschen sind bereit, das für BLAS zu tun Routinen, und stellen Sie dann die abgestimmten Ergebnisse zur Verfügung, bedeutet, dass Sie am besten diese Routinen verwenden.

Das Speicherzugriffsmuster für GEMV ist so einfach, dass Sie Divide and Conquer nicht wirklich brauchen (dasselbe gilt für den Standardfall der Matrixtransponierung) - Sie finden nur die Cache-Blockgröße und verwenden sie. In GEMV (y = Ax) müssen Sie immer noch einmal die gesamte Matrix durchsehen, damit Sie nichts für die Wiederverwendung (und damit die effektive Cache-Nutzung) tun können, aber Sie können versuchen, x so oft wie möglich wiederzuverwenden, um es zu laden einmal statt (Anzahl der Zeilen) mal - und du willst immer noch Zugriff auf A um cachefreundlich zu sein. So ist die offensichtliche Cache-Blockierung Sache zu tun, um Blöcke zu brechen:

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

Und Sie können das rekursiv tun. Aber eine naive Implementierung zu tun, dann ist es langsamer als der einfache Doppel-Loop und Weise langsamer als ein richtiger SGEMV Bibliotheksaufruf:

$ ./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 

Der Code folgt:

#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; 
}