Как ускорить мой решатель разреженной матрицы?

Думаю, мне удалось его оптимизировать, вот код, создать новый проект на VC ++, добавить этот код и просто скомпилировать в разделе «Выпуск».

#include <iostream>
#include <cstdlib>
#include <cstring>

#define _WIN32_WINNT 0x0400
#define WIN32_LEAN_AND_MEAN
#include <windows.h>

#include <conio.h>

using namespace std;

size_t d_nx = 128, d_ny = 128;
float *d_x, *d_b, *d_w, *d_e, *d_s, *d_n;

void step_original() {
    size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
    for (size_t y = 1; y < d_ny - 1; ++y) {
        for (size_t x = 1; x < d_nx - 1; ++x) {
            d_x[ic] = d_b[ic]
                - d_w[ic] * d_x[iw] - d_e[ic] * d_x[ie]
                - d_s[ic] * d_x[is] - d_n[ic] * d_x[in];
            ++ic; ++iw; ++ie; ++is; ++in;
        }
        ic += 2; iw += 2; ie += 2; is += 2; in += 2;
    }
}
void step_new() {
    //size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
    float
        *d_b_ic,
        *d_w_ic,
        *d_e_ic,
        *d_x_ic,
        *d_x_iw,
        *d_x_ie,
        *d_x_is,
        *d_x_in,
        *d_n_ic,
        *d_s_ic;

    d_b_ic = d_b;
    d_w_ic = d_w;
    d_e_ic = d_e;
    d_x_ic = d_x;
    d_x_iw = d_x;
    d_x_ie = d_x;
    d_x_is = d_x;
    d_x_in = d_x;
    d_n_ic = d_n;
    d_s_ic = d_s;

    for (size_t y = 1; y < d_ny - 1; ++y)
    {
        for (size_t x = 1; x < d_nx - 1; ++x)
        {
            /*d_x[ic] = d_b[ic]
                - d_w[ic] * d_x[iw] - d_e[ic] * d_x[ie]
                - d_s[ic] * d_x[is] - d_n[ic] * d_x[in];*/
            *d_x_ic = *d_b_ic
                - *d_w_ic * *d_x_iw - *d_e_ic * *d_x_ie
                - *d_s_ic * *d_x_is - *d_n_ic * *d_x_in;
            //++ic; ++iw; ++ie; ++is; ++in;
            d_b_ic++;
            d_w_ic++;
            d_e_ic++;
            d_x_ic++;
            d_x_iw++;
            d_x_ie++;
            d_x_is++;
            d_x_in++;
            d_n_ic++;
            d_s_ic++;
        }
        //ic += 2; iw += 2; ie += 2; is += 2; in += 2;
        d_b_ic += 2;
        d_w_ic += 2;
        d_e_ic += 2;
        d_x_ic += 2;
        d_x_iw += 2;
        d_x_ie += 2;
        d_x_is += 2;
        d_x_in += 2;
        d_n_ic += 2;
        d_s_ic += 2;
    }
}

void solve_original(size_t iters) {
    for (size_t i = 0; i < iters; ++i) {
        step_original();
    }
}
void solve_new(size_t iters) {
    for (size_t i = 0; i < iters; ++i) {
        step_new();
    }
}

void clear(float *a) {
    memset(a, 0, d_nx * d_ny * sizeof(float));
}

int main(int argc, char **argv) {
    size_t n = d_nx * d_ny;
    d_x = new float[n]; clear(d_x);
    d_b = new float[n]; clear(d_b);
    d_w = new float[n]; clear(d_w);
    d_e = new float[n]; clear(d_e);
    d_s = new float[n]; clear(d_s);
    d_n = new float[n]; clear(d_n);

    if(argc < 3)
        printf("app.exe (x)iters (o/n)algo\n");

    bool bOriginalStep = (argv[2][0] == 'o');
    size_t iters = atoi(argv[1]);

    /*printf("Press any key to start!");
    _getch();
    printf(" Running speed test..\n");*/

    __int64 freq, start, end, diff;
    if(!::QueryPerformanceFrequency((LARGE_INTEGER*)&freq))
        throw "Not supported!";
    freq /= 1000000; // microseconds!
    {
        ::QueryPerformanceCounter((LARGE_INTEGER*)&start);
        if(bOriginalStep)
            solve_original(iters);
        else
            solve_new(iters);
        ::QueryPerformanceCounter((LARGE_INTEGER*)&end);
        diff = (end - start) / freq;
    }
    printf("Speed (%s)\t\t: %u\n", (bOriginalStep ? "original" : "new"), diff);
    //_getch();


    //cout << d_x[0] << endl; // prevent the thing from being optimized away
}

Запустите его так:

app.exe 10000 o

app.exe 10000 n

«o» означает старый код, ваш.

«n» - мое, новое.

Мои результаты: Скорость (исходная):

1515028

1523171

1495988

Скорость (новая):

966012

984110

1006045

Улучшение около 30%.

Логика, лежащая в основе: Вы использовали счетчики индекса для доступа / управления. Я использую указатели. Во время работы точка останова в определенной строке кода вычисления в VC ++ отладчик и нажмите F8. Вы получите окно дизассемблера. Вы увидите созданные коды операций (код сборки).

В любом случае, посмотрите:

int * x = ...;

x [3] = 123;

Это указывает компьютеру поместить указатель x в регистр (скажем, EAX). {{ 1}} Добавить (3 * sizeof (int)). Только тогда установите значение 123.

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

Надеюсь, это поможет.

Замечание сотрудникам stackoverflow.com: Отличный веб-сайт, надеюсь, я слышал о нем давно!

Во-первых, здесь, похоже, есть проблема конвейерной обработки. Цикл читает из значения в d_x , которое только что было записано, но, очевидно, он должен дождаться завершения этой записи. Простое изменение порядка вычислений, выполнение чего-то полезного во время ожидания делает его почти вдвое быстрее:

d_x[ic] = d_b[ic]
                        - d_e[ic] * d_x[ie]
    - d_s[ic] * d_x[is] - d_n[ic] * d_x[in]
    - d_w[ic] * d_x[iw] /* d_x[iw] has just been written to, process this last */;

Это Имон Нербон понял это. Большое ему спасибо! Никогда бы не догадался.

Ответ Пони мне кажется правильным.

Я просто хочу отметить, что в задачах этого типа вы часто получаете преимущества от локальности памяти. Прямо сейчас массивы b, w, e, s, n находятся в разных местах памяти. Если бы вы могли не уместить проблему в кэше L3 (в основном в L2), это было бы плохо, и решение такого рода было бы полезным:

size_t d_nx = 128, d_ny = 128;
float *d_x;

struct D { float b,w,e,s,n; };
D *d;

void step() {
    size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
    for (size_t y = 1; y < d_ny - 1; ++y) {
        for (size_t x = 1; x < d_nx - 1; ++x) {
            d_x[ic] = d[ic].b
                - d[ic].w * d_x[iw] - d[ic].e * d_x[ie]
                - d[ic].s * d_x[is] - d[ic].n * d_x[in];
            ++ic; ++iw; ++ie; ++is; ++in;
        }
        ic += 2; iw += 2; ie += 2; is += 2; in += 2;
    }
}
void solve(size_t iters) { for (size_t i = 0; i < iters; ++i) step(); }
void clear(float *a) { memset(a, 0, d_nx * d_ny * sizeof(float)); }

int main(int argc, char **argv) {
    size_t n = d_nx * d_ny;
    d_x = new float[n]; clear(d_x);
    d = new D[n]; memset(d,0,n * sizeof(D));
    solve(atoi(argv[1]));
    cout << d_x[0] << endl; // prevent the thing from being optimized away
}

Например, это решение с разрешением 1280x1280 является чуть менее чем в 2 раза быстрее, чем решение Пони (13 с против 23 в моем тесте - тогда ваша исходная реализация составляет 22 с), а при 128x128 оно на 30% медленнее (7 с против 10 с - -ваше оригинал 10с).

(Итерации были увеличены до 80000 для базового случая и до 800 для более крупного случая 1280x1280 в 100 раз.)

Думаю, вы правы в том, что память является узким местом. Это довольно простой цикл с простой арифметикой на итерацию. ic, iw, т.е. есть, и in index, кажется, находятся на противоположных сторонах матрицы, поэтому я предполагаю, что там есть куча промахов кеша.

Я предлагаю добавить несколько операторов предварительной выборки, а также исследовать «дизайн, ориентированный на данные»:

void step_original() {
    size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
    float dw_ic, dx_ic, db_ic, de_ic, dn_ic, ds_ic;
    float dx_iw, dx_is, dx_ie, dx_in, de_ic, db_ic;
    for (size_t y = 1; y < d_ny - 1; ++y) {
        for (size_t x = 1; x < d_nx - 1; ++x) {
// Perform the prefetch
// Sorting these statements by array may increase speed;
//    although sorting by index name may increase speed too.
            db_ic = d_b[ic];
            dw_ic = d_w[ic];
            dx_iw = d_x[iw];
            de_ic = d_e[ic];
            dx_ie = d_x[ie];
            ds_ic = d_s[ic];
            dx_is = d_x[is];
            dn_ic = d_n[ic];
            dx_in = d_x[in];
// Calculate
            d_x[ic] = db_ic
                - dw_ic * dx_iw - de_ic * dx_ie
                - ds_ic * dx_is - dn_ic * dx_in;
            ++ic; ++iw; ++ie; ++is; ++in;
        }
        ic += 2; iw += 2; ie += 2; is += 2; in += 2;
    }
}

Это отличается от вашего второй метод, поскольку значения копируются в локальные временные переменные перед выполнением вычисления.

Я не специалист по этому вопросу, но я видел, что есть несколько научных статей об улучшении использования кеша методом Гаусса-Зейделя.

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