Алгоритм оптимизации матриц и векторов - PullRequest
Новая
       29

Алгоритм оптимизации матриц и векторов

1 голос
/ 28 марта 2012

Предположим, что размеры очень велики (до 1 миллиарда элементов в матрице). Как бы я реализовал алгоритм кеширования для матрично-векторного произведения? Основываясь на википедии, мне нужно будет рекурсивно разделять и завоевывать, но я чувствую, что будет много накладных расходов ... Было бы эффективно сделать это?

Последующие вопросы и ответы: OpenMP с матрицами и векторами

1 Ответ

3 голосов
/ 28 марта 2012

Таким образом, ответ на вопрос «как сделать эту базовую операцию линейной алгебры быстрым» всегда и везде можно найти и связать с настроенной библиотекой BLAS для вашей платформы.Например, GotoBLAS (работа которого продолжается в OpenBLAS ), или более медленная автонастройка ATLAS , или коммерческие пакеты, такие как Intel MKL.Линейная алгебра настолько фундаментальна для многих других операций, что огромные усилия направляются на оптимизацию этих пакетов для различных платформ, и нет никаких шансов, что вы придете к чему-то за несколько дней работы, которая будет конкурировать.Конкретные вызовы подпрограмм, которые вы ищете для общего умножения матрицы-вектора: SGEMV / DGEMV / CGEMV / ZGEMV.

Кэш-забывающие алгоритмы или автонастройка предназначены для случаев, когда вам не нужно настраиваться наконкретная архитектура кеша вашей системы - обычно это может быть хорошо, но поскольку люди готовы сделать это для подпрограмм BLAS, а затем сделать настроенные результаты доступными, это означает, что вам лучше всего просто использоватьэти процедуры.

Шаблон доступа к памяти для GEMV достаточно прост, так что вам не нужно делить и захватывать (то же самое для стандартного случая транспонирования матрицы) - вы просто находите размер блокировки кеша и используете его.В GEMV (y = Ax) вам все равно придется сканировать всю матрицу один раз, поэтому для повторного использования (и, следовательно, для эффективного использования кэша) там ничего не поделаешь, но вы можете попытаться использовать x как можно больше, чтобы загрузить егоодин раз вместо (количество строк) - и вы все еще хотите, чтобы доступ к A был дружественным к кешу.Таким образом, очевидная вещь, блокирующая кэш, состоит в том, чтобы разбивать блоки: 

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

И вы, конечно, можете сделать это рекурсивно.Но, делая простую реализацию, она медленнее, чем простой двойной цикл, и намного медленнее, чем правильный вызов библиотеки SGEMV: 

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

Код выглядит следующим образом: 

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