Форум программистов, компьютерный форум CyberForum.ru

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

Восстановить пароль Регистрация
 
mdt::Vladimir
7 / 7 / 0
Регистрация: 23.09.2011
Сообщений: 32
18.11.2011, 19:01     Метод Гаусса, valarray, c-array и sse2 оптимизация, скорость решения #1
Не знаю, в правильном ли месте задаю вопрос. Дело в том, что есть одна процедура в большой программе, на выполнение которой уходит 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);
}
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
18.11.2011, 19:01     Метод Гаусса, valarray, c-array и sse2 оптимизация, скорость решения
Посмотрите здесь:

C++ метод гаусса для решения линейних уравнений
Метод Гаусса для решения СЛАУ с использованием одномерного массива C++
Оптимизация решения. C++
C++ Разработать программу для решения СЛАУ методом Гаусса.
Метод деления отрезка пополам для решения нелинейных уравнений (метод дихотомии) C++
C++ Метод Гаусса решения СЛАУ с полным выбором. C++
C++ Решения СЛАУ методом Гаусса по шагу
C++ Решения матрицы методом Гаусса

Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Ответ Создать тему
Опции темы

Текущее время: 13:48. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru