Оптимизация алгоритма Гаусса

@volver · Регистрация: 18.09.2021

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

можно ли оптимизировать данный код еще лучше, используя SIMD-инструкции?

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
#include <iostream>
#include <vector>
#include <immintrin.h>
#include <chrono>
 
void gaussEliminationSIMD(std::vector<std::vector<double>>& A, std::vector<double>& b) {
    int n = A.size();
 
    for (int i = 0; i < n; ++i) {
        // Приведение матрицы к верхнетреугольному виду
        for (int j = i + 1; j < n; ++j) {
            double factor = A[j][i] / A[i][i];
            b[j] -= factor * b[i];
            for (int k = i; k < n; ++k) {
                A[j][k] -= factor * A[i][k];
            }
        }
    }
 
    // Обратный проход для нахождения решения
    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        __m256d sum = _mm256_setzero_pd();
        for (int j = i + 1; j < n; ++j) {
            __m256d Av = _mm256_loadu_pd(&A[i][j]);
            __m256d xv = _mm256_broadcast_sd(&x[j]);
            sum = _mm256_add_pd(sum, _mm256_mul_pd(Av, xv));
        }
 
        x[i] = (b[i] - _mm256_cvtsd_f64(_mm256_hadd_pd(sum, sum))) / A[i][i];
    }
 
    // Вывод решения
    std::cout << "Решение x:\n";
    for (int i = 0; i < n; ++i) {
        std::cout << x[i] << " ";
    }
    std::cout << std::endl;
}
 
int main() {
    setlocale(LC_ALL, "rus");
    // Пример использования
    std::vector<std::vector<double>> A = { {1043000, 134300, 434399}, {123445, 453421, 112434}, {1213212, 5443232, 8746199} };
    std::vector<double> b = { 80009990, 126606090, 346060990 };
 
    // Измерение времени выполнения
    auto start = std::chrono::high_resolution_clock::now();
 
    gaussEliminationSIMD(A, b);
 
    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    std::cout << "Время выполнения: " << duration.count() << " секунд." << std::endl;
 
    return 0;
}

@VTsaregorodtsev · 30.11.2023, 15:34

Матрица 3*3? Смешно.
1000*1000 хотя-бы возьми - там даже скалярный (т.е. без векторизации компилятором/вручную) сишный код должен не дольше чем за пару секунд отрабатывать.

@volver · 30.11.2023, 17:04 **[ТС]**

VTsaregorodtsev, сделал 1000*1000
теперь возникает такая ошибка.

код теперь такой:

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
68
69
70
71
72
#include <iostream>
#include <vector>
#include <chrono>
#include <cstdlib>
#include <ctime>
#include <immintrin.h>
 
void gaussElimination(std::vector<std::vector<double>>& A, std::vector<double>& b) {
    int n = A.size();
 
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            double factor = A[j][i] / A[i][i];
 
            __m256d factor_vec = _mm256_set1_pd(factor);
            for (int k = i; k < n; k += 4) {
                __m256d A_row = _mm256_loadu_pd(&A[i][k]);
                __m256d factor_A = _mm256_mul_pd(factor_vec, A_row);
 
                __m256d B_row = _mm256_loadu_pd(&A[j][k]);
                B_row = _mm256_sub_pd(B_row, factor_A);
 
                _mm256_storeu_pd(&A[j][k], B_row);
            }
 
            b[j] -= factor * b[i];
        }
    }
 
    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        double sum = 0.0;
        for (int j = i + 1; j < n; ++j) {
            sum += A[i][j] * x[j];
        }
 
        x[i] = (b[i] - sum) / A[i][i];
    }
 
    std::cout << "Решение x:\n";
    for (int i = 0; i < n; ++i) {
        std::cout << x[i] << " ";
    }
    std::cout << std::endl;
}
 
int main() {
    setlocale(LC_ALL, "rus");
    int n = 1000;
 
    std::vector<std::vector<double>> A(n, std::vector<double>(n));
    std::vector<double> b(n);
 
    std::srand(std::time(0));
 
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            A[i][j] = 1.0 + static_cast<double>(rand()) / (RAND_MAX / 9.0); // Значения от 1.0 до 10.0
        }
        b[i] = 1.0 + static_cast<double>(rand()) / (RAND_MAX / 9.0); // Значения от 1.0 до 10.0
    }
 
    auto start = std::chrono::high_resolution_clock::now();
 
    gaussElimination(A, b);
 
    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    std::cout << "Время выполнения: " << duration.count() << " секунд." << std::endl;
 
    return 0;
}

@VTsaregorodtsev · 30.11.2023, 18:06

Сообщение от volver

теперь возникает такая ошибка.
код теперь такой:

Выход за пределы массива. На первой итерации циклов в строках 31 и 33 (там при i=n-1 у внешнего цикла - на внутреннем цикле j получает значение j=i+1=n-1+1=n, а элементы массива имеют индексы от 0 до n-1).
Одну ошибку я нашёл - дальше другие возможные ошибки мне лень искать.

@volver · 30.11.2023, 22:05 **[ТС]**

VTsaregorodtsev, а как тогда цикл должен выглядеть?

@Aledveu · 01.12.2023, 01:45

volver, прежде чем что-то оптимизировать надо выяснить на что тратится время.
То что вы поставили пару simd инструкций не значит что стало быстрее.
Вы проверили что генерирует компилятор если просто на С написать код ? со всеми оптимизациями.
соответственно запустите дизассемблер и сравните.
Вы проверили последовательность обращения к памяти - для наилучшего использования кэша ?
замеры производительности надо производить после первого прогона - когда всё уже в кеше и на временной базе когда погрешность измерения времени хотябы в 1000 раз меньше измеряемой величины.
но тогда от 30, 40-45 строки надо избавиться - иначе вы будете замерять ещё и выделение памяти и вывод на экран ))

@Alexis · 01.12.2023, 02:52

volver, Вы запрограммировали неустойчивый алгоритм, если на очередном шаге случится так, что A[i][i]==0,
то он потерпит неудачу, если же A[i][i] будет равно очень малому числу, то будет потеряна вся точность и решение будет неадекватным (этот вариант даже хуже, так как его можно вовремя и не обнаружить).

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

@volver · 01.12.2023, 03:16 **[ТС]**

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

Добавлено через 6 минут
Alexis333, я сейчас добавил такие строки:

C++
1
2
3
 for (int l = k; l < n; ++l) {
            A[j][l] -= factor * A[i][l];
        }

Тоесть просто проверку, что мы не выходим за пределы массива при использовании SIMD-инструкций для блока размером 4 элемента. Так правильнее будет?

Добавлено через 14 минут
Alexis333, прошу прощения, но есть еще один вопрос: я попробовал запустить на линуксе, и "оптимизированный" алгоритм работает примерно в 2 раза дольше

не знаете, в чем может быть проблема?

Добавлено через 14 секунд
Alexis333, прошу прощения, но есть еще один вопрос: я попробовал запустить на линуксе, и "оптимизированный" алгоритм работает примерно в 2 раза дольше

не знаете, в чем может быть проблема?

@Alexis · 01.12.2023, 04:34

Сообщение от volver

что лучше сделать, чтобы не терять точность?

на этот вопрос ответ известен - специально выбирать A[i][i] путём перестановки строк или столбцов, лучше всего делить на максимальный по модулю элемент необработанного фрагмента матрицы. Но это заметно снизит производительность. Лучше всего просто переставлять строки, ставя каждый раз вперёд строку с максимальным по модулю элементом в оставшейся части ведущего столбца. Такой алгоритм и быстрый и достаточно точный.

Сообщение от volver

с помощью SIMD-инструкций, что лучше сделать, чтобы не терять точность?

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

Ещё один вариант - использовать SIMD для получения приближенного решения, которое потом уточнить с помощью нескольких простой итераций. Вычислительная сложность алгоритма Гаусса, на сколько я помню 3n^3 сложений-умножений, а простой итерации - всего n^2 сложений-умножений. Так можно сразу повысить и точность и производительность.

@volver · 01.12.2023, 04:38 **[ТС]**

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

@Alexis · 01.12.2023, 05:25

volver, попробуйте использовать библиотеки blas

@igorrr37 · 10.12.2023, 15:13

Сообщение от volver

теперь возникает такая ошибка

Запись за пределы массива в строке 23: разница n-k может быть меньше чем 4 а ф-ция _mm256_storeu_pd записывает 4 дабла в &A[j][k]. Хвост массива надо досчитывать обычным циклом если его длина не кратна четырём. Если вместо double взять float то в регистр влезет 8 чисел вместо 4 что должно дать ускорение.

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <iostream>
#include <chrono>
#include <vector>
#include <cassert>
#include <immintrin.h>
 
void gaussElimination(double* A, double* b, int n)
{
    for (int i = 0; i < n; ++i)
    {
        for (int j = i + 1; j < n; ++j)
        {
            double factor = -A[j * n + i] / A[i * n + i];
            b[j] = b[j] + factor * b[i];
 
            __m256d factor256 = _mm256_broadcast_sd(&factor);
 
            int k = i;
            for (; k <= n - 4; k += 4) // 4 дабла за раз
            {
                _mm256_storeu_pd(A + j * n + k, _mm256_fmadd_pd(factor256, _mm256_loadu_pd(A + i * n + k), _mm256_loadu_pd(A + j * n + k)));
            }
            for (; k < n; ++k) // добиваем хвост если (n - i) не кратно 4
            {
                A[j * n + k] = A[j * n + k] + factor * A[i * n + k];
            }
        }
    }
 
    double* x = new double[n];
    for (int i = n - 1; i >= 0; --i)
    {
        __m256d sum256 = _mm256_setzero_pd();
        int j = i + 1;
        for (; j <= n - 4; j += 4) // 4 дабла за раз
        {
            sum256 = _mm256_fmadd_pd(_mm256_loadu_pd(A + i * n + j), _mm256_loadu_pd(x + j), sum256);
        }
        auto psum = (double*)&sum256;
        double sum = psum[0] + psum[1] + psum[2] + psum[3];
        for (; j < n; ++j) // добиваем хвост если (n - (i + 1)) не кратно 4
        {
            sum = sum + A[i * n + j] * x[j];
        }
        x[i] = (b[i] - sum) / A[i * n + i];
    }
 
    std::cout << "Решение x:\n";
    for (int i = n - 5; i < n; ++i) // выводим не весь результат
    {
        std::cout << x[i] << " ";
    }
    std::cout << std::endl;
 
    delete[] x;
}
 
int main()
{
    setlocale(LC_ALL, "rus");
    int n = 2000;
    double* A = new double[n * n]; // заменил матрицу на массив как в гайде на хабре, но похоже скорость от этого не улучшилась
    double* b = new double[n];
    for (int i = 0; i < n * n; ++i)
    {
        A[i] = 1.0 + static_cast<double>(rand()) / (RAND_MAX / 9.0);
    }
    for (int i = 0; i < n; ++i)
    {
        b[i] = 1.0 + static_cast<double>(rand()) / (RAND_MAX / 9.0);
    }
 
    auto start = std::chrono::high_resolution_clock::now();
 
    gaussElimination(A, b, n);
 
    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    std::cout << "Время выполнения: " << duration.count() << " секунд." << std::endl;
 
    delete[] A;
    delete[] b;
}

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): сборка C/C++ проекта из консоли 8Observer8 30.01.2026 Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .	Установка Emscripten SDK (emsdk) и CMake на Windows для сборки C и C++ приложений в WebAssembly (Wasm) 8Observer8 30.01.2026 Чтобы скачать Emscripten SDK (emsdk) необходимо сначало скачать и уставить Git: Install for Windows. Следуйте стандартной процедуре установки Git через установщик. Система контроля версиями Git. . .	Подключение Box2D v3 к SDL3 для Android: физика и отрисовка коллайдеров 8Observer8 29.01.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами. Версия v3 была полностью переписана на Си, в. . .	Инструменты COM: Сохранение данный из VARIANT в файл и загрузка из файла в VARIANT bedvit 28.01.2026 Сохранение базовых типов COM и массивов (одномерных или двухмерных) любой вложенности (деревья) в файл, с возможностью выбора алгоритмов сжатия и шифрования. Часть библиотеки BedvitCOM Использованы. . .
Загрузка PNG с альфа-каналом на SDL3 для Android: с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 28.01.2026 Содержание блога SDL3 имеет собственные средства для загрузки и отображения PNG-файлов с альфа-каналом и базовой работы с ними. В этой инструкции используется функция SDL_LoadPNG(), которая. . .	Загрузка PNG с альфа-каналом на SDL3 для Android: с помощью SDL3_image 8Observer8 27.01.2026 Содержание блога SDL3_image - это библиотека для загрузки и работы с изображениями. Эта пошаговая инструкция покажет, как загрузить и вывести на экран смартфона картинку с альфа-каналом, то есть с. . .	Влияние грибов на сукцессию anaschu 26.01.2026 Бифуркационные изменения массы гриба происходят тогда, когда мы уменьшаем массу компоста в 10 раз, а скорость прироста биомассы уменьшаем в три раза. Скорость прироста биомассы может уменьшаться за. . .	Воспроизведение звукового файла с помощью SDL3_mixer при касании экрана Android 8Observer8 26.01.2026 Содержание блога SDL3_mixer - это библиотека я для воспроизведения аудио. В отличие от инструкции по добавлению текста код по проигрыванию звука уже содержится в шаблоне примера. Нужно только. . .

@VTsaregorodtsev 2637 / 1648 / 267 Регистрация: 19.02.2010 Сообщений: 4,368
	30.11.2023, 15:34
	Матрица 33? Смешно. 10001000 хотя-бы возьми - там даже скалярный (т.е. без векторизации компилятором/вручную) сишный код должен не дольше чем за пару секунд отрабатывать. 0

@volver 2 / 2 / 0 Регистрация: 18.09.2021 Сообщений: 418
	30.11.2023, 22:05 [ТС]
	VTsaregorodtsev, а как тогда цикл должен выглядеть? 0

@Aledveu 267 / 199 / 30 Регистрация: 26.11.2022 Сообщений: 866
	01.12.2023, 01:45
	volver, прежде чем что-то оптимизировать надо выяснить на что тратится время. То что вы поставили пару simd инструкций не значит что стало быстрее. Вы проверили что генерирует компилятор если просто на С написать код ? со всеми оптимизациями. соответственно запустите дизассемблер и сравните. Вы проверили последовательность обращения к памяти - для наилучшего использования кэша ? замеры производительности надо производить после первого прогона - когда всё уже в кеше и на временной базе когда погрешность измерения времени хотябы в 1000 раз меньше измеряемой величины. но тогда от 30, 40-45 строки надо избавиться - иначе вы будете замерять ещё и выделение памяти и вывод на экран )) 0

@Alexis 524 / 512 / 129 Регистрация: 31.10.2016 Сообщений: 4,156
	01.12.2023, 02:52
	volver, Вы запрограммировали неустойчивый алгоритм, если на очередном шаге случится так, что A[i][i]==0, то он потерпит неудачу, если же A[i][i] будет равно очень малому числу, то будет потеряна вся точность и решение будет неадекватным (этот вариант даже хуже, так как его можно вовремя и не обнаружить). Ну и обратный проход, который Вы оптимизируете, занимает намного меньше вычислительных затрат, чем прямой проход. 1

@volver 2 / 2 / 0 Регистрация: 18.09.2021 Сообщений: 418
	01.12.2023, 04:38 [ТС]
	Alexis333, в общем simd инструкции не дали мне никакого увеличения в оптимизации. я попробовал использовать встроенную библиотеку OpenMP, так стало лучше. хотело бы конечно самому это оптимизировать, но что-то безрезультативно 0

@Alexis 524 / 512 / 129 Регистрация: 31.10.2016 Сообщений: 4,156
	01.12.2023, 05:25
	volver, попробуйте использовать библиотеки blas 0