Метод Гаусса, valarray, c-array и sse2 оптимизация, скорость решения

@mdt::Vladimir · Регистрация: 23.09.2011

Студворк — интернет-сервис помощи студентам

Не знаю, в правильном ли месте задаю вопрос. Дело в том, что есть одна процедура в большой программе, на выполнение которой уходит 85% времени работы программы.
На всякий случай напишу про программу - это система моделирования систем электроснабжения (СЭС), в ней составляется система линейных дифференциальных уравнений, которые решаются методом Рунге-Кутты, на каждой итерации составляется система уравнений, которая приводится к задаче Коши, т.е. все производные выносятся в левую часть уравнений, для этого требуется решить систему линейных уравнений, чтобы выразить производные. И на каждом шаге происходит решение, обычно это около 50 уравнений в системе. За всё время работы решается несколько сотен тысяч таких систем, а в дальнейшем может потребоваться и такое применение программы, что понадобится и сотни миллионов раз решить. Естественно, алгоритм решения захотелось сделать максимально эффективным.
Вот код, который используется для решения:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
dvalarray linear_solve_c(const dvalarray &M, const size_t n)
{
    //копируем матрицу и создаем матрицу результатов
    dvalarray R(n);
    dvalarray A(M);
    double *a = &A[0];
    const size_t m = n+1;
    slice diag(0,n,m+1); //диагональ матрицы
    //приводим матрицу к треугольной
    for (size_t i=0; i<n-1; ++i)
    {
        slice s0(i*(m+1), m-i, 1);
        if (A[slice_index(diag, i)] == 0.0)
        {
            slice t(slice_col_c(n, m, i));
            size_t k;
            for (k=i+1; k<n; ++k)
                if (A[slice_index(t, k)] != 0)
                {
                    slice z(m*k+i, m-i, 1);
                    dvalarray tmp(A[z]);
                    A[z] = A[s0];
                    A[s0] = tmp;
                    break;
                }
                if (k == n)
                    break;
        }
#define MARRAY
#ifdef MARRAY
        double mul;
        size_t d = m+1, i0 = d*i, j0 = i0 + m, jd = m - i, jn = j0+jd-m, j;
#ifdef MP
        #pragma omp parallel for private(mul, j, jn) num_threads(8)
#endif
        for (j = j0; j<n*m; j += m) //на этот цикл уходит 80% всего времени
        {
#ifdef MP
            jn = j+jd;
#else
            jn += m;
#endif
            mul = a[j] / a[i0];
            double *a0 = a + i0, *ai = a + j;
            for (size_t k = 0; k<jd; ++k)
                ai[k] -= a0[k] * mul;
        }
#else
        for (size_t j=i+1; j<n; ++j)
        {
            slice s1(m*j+i, m-i, 1);
 
            double mul = A[slice_index(s1, 0)] / A[slice_index(diag, i)];
            A[s1] -= dvalarray(A[s0]) * mul;
        }
#endif
    }
    //вычисляем результаты
    slice right(slice_col_c(n, m, n));
    for (int i = n-1; i>=0; --i)
    {
        slice cur(slice_col_c(n, m, i));
        R[i] = A[slice_index(right, i)] / A[slice_index(cur, i)];
        A[right] -= dvalarray(A[cur]) * R[i];
    }
    return R;
}

dvalarray - это typedef valarray<double>
Параметры функции: 1 - сама система, представленная в виде матрицы, 2 - количество уравнений. valarray содержит n*(n+1) чисел, строки матрицы расположены в массиве последовательно.
Вообще самое узкое место, на которое уходит практически всё время - это то, что заключено в директивы #ifdef MARRAY - #else - #endif, то есть один цикл. Просто проверял, как работает в различных реализациях.
Использую gcc для компиляции, и вопрос собственно состоит из нескольких частей:
1. При определённом MARRAY решение происходит в два раза быстрее, т.е. при использовании обычного массива быстрее, чем valarray, хотя по сути происходит одно и то же.
2. При использовании автоматической оптимизации sse2 скорость нисколько не увеличивается, хотя обычно в других реализациях, например умножение матриц, я всегда замечал прирост скорости ровно в 2 раза.
3. Совсем необязательный, да и тут может быть влияют другие части программы, но при мереходе с Core 2 Duo 3 ГГц на Core i7 4.7 ГГц скорость увеличилась в 60 раз. Правда в целом у всей программы, а она многопоточная, т.е. на каждое ядро увеличение скорости в 30 раз, хотя частота возросла в 1.5 раза.
Кто-нибудь может объяснить, как такое может быть?
Ну также я пробовал использовать OpenMP, но даже на четырёхъядерном процессоре скорость решения возрастает только при более 100 уравнений в системе (да и пока не очень хорошо оптимизировал использование многопоточности здесь). Но и это не вариант, т.к. в целом в программе используется оптимизация СЭС при помощи генетического алгоритма, а там как раз несколько СЭС моделируются параллельно.

Недостающие функции

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
inline size_t slice_index(const slice &s, size_t index)
{
    return s.start() + s.stride() * index;
}
 
inline slice slice_row_c(size_t rows, size_t cols, size_t index)
{
    rows = rows; //disable warning
    return slice(index*cols, cols, 1);
}
 
inline slice slice_col_c(size_t rows, size_t cols, size_t index)
{
    return slice(index, rows, cols);
}

Новые блоги и статьи Все статьи Все блоги /
PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .
Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .