Оптимизировать код

[Nick]

Всем привет.
Очень хотелось бы понять как оптимизировать код:

Код: (C++)

    #include <vector>

    using namespace std;  

    double R=330.3;
    double A=400;
    vector<double> tmp_const(M), pieces(M), x(M), y(M), z(M), r(M);
    vector<double> global(N), first_summand(N), summands(N);
    vector<vector<double> > tmp_matrix(M, std::vector<double>(M));  
    vector<double> params(N),result(N);


    double vconst = 4/3 * M_PI;
    double V = vconst * pow(R,3);

    int k, j, i;
//const
    for( j = 0; j < M; j++)
    {
        x[j] = data[4 * j];
        y[j] = data[1 + j * 4];
        z[j] = data[2 + j * 4];
        r[j] = data[3 + j * 4];      
        pieces[j]= vconst * pow(r[j], 3);
        tmp_const[j] = sqrt(pow(x[j], 2) + pow(y[j], 2)+ pow(z[j], 2)) / M_PI;
    }

   
//matr
    for( j = 0; j < M; j++)
    {
        for(i = 0; i < M; i++)
        {
            tmp_matrix[j][i] = sqrt(pow(x[j]-x[i], 2) + pow(y[j]-y[i], 2) + pow(z[j]-z[i], 2)) / M_PI;
        }
    }
 
//
    for(k = 0; k < N; k++)
    {  
        summands[k] = 0.0;
        result[k] = 0.0;
        first_summand[k] = 0.0;
        global[k] = cos(R * params[k]) * V;

        for(j = 0; j < M; j++)
        {
            first_summand[k] += pieces[j] * cos(r[j] * params[k]) * sin(params[k] * tmp_const[j]);
            for(i = 0; i < M; i++)
            {
                summands[k] += pieces[i] * pieces[j] * cos(r[j] * params[k]) * cos(r[i] * params[k]) * sin(params[k] * tmp_matrix[j][i]);
             }
        }
        result[k] = summands[k] + pow(A ,2) * pow(global[k], 2) - 2 * A * first_summand[k] * global[k];
    }
 

Критический момент то, что при достаточно больших M(>40000)
матрица tmp_matrix достигает такого размера что памяти тупо перестает хватать.

Планируется переписать этот код на чистый C. Соответственно, буду признателен за любые советы по оптимизации.


Всем СПАСИБО!

[darkelf]

Самый простой вариант переделки Вашего варианта, если надо влезть в память - не вычислять tmp_matrix:

Код: (C++)

    #include <vector>

    using namespace std;  

    double R=330.3;
    double A=400;
    vector<double> tmp_const(M), pieces(M), x(M), y(M), z(M), r(M);
    vector<double> global(N), first_summand(N), summands(N);
    vector<double> params(N),result(N);

    double vconst = 4/3 * M_PI;
    double V = vconst * pow(R,3);
    double tmp_value;

    int k, j, i;
//const
    for( j = 0; j < M; j++)
    {
        x[j] = data[4 * j];
        y[j] = data[1 + j * 4];
        z[j] = data[2 + j * 4];
        r[j] = data[3 + j * 4];      
        pieces[j]= vconst * pow(r[j], 3);
        tmp_const[j] = sqrt(pow(x[j], 2) + pow(y[j], 2)+ pow(z[j], 2)) / M_PI;
    }
//
    for(k = 0; k < N; k++)
    {  
        summands[k] = 0.0;
        result[k] = 0.0;
        first_summand[k] = 0.0;
        global[k] = cos(R * params[k]) * V;

        for(j = 0; j < M; j++)
        {
            first_summand[k] += pieces[j] * cos(r[j] * params[k]) * sin(params[k] * tmp_const[j]);
            for(i = 0; i < M; i++)
            {
                tmp_value = sqrt(pow(x[j]-x[i], 2) + pow(y[j]-y[i], 2) + pow(z[j]-z[i], 2)) / M_PI;
                summands[k] += pieces[i] * pieces[j] * cos(r[j] * params[k]) * cos(r[i] * params[k]) * sin(params[k] * tmp_value);
            }
        }
        result[k] = summands[k] + pow(A ,2) * pow(global[k], 2) - 2 * A * first_summand[k] * global[k];
    }
 

т.е. придётся в M раз чаще вычислять значение ячейки tmp_matrix[ j ][ i ], но зато не выделять M*M*sizeof(double) байтов памяти.

Критический момент то, что при достаточно больших M(>40000)
матрица tmp_matrix достигает такого размера что памяти тупо перестает хватать.

Планируется переписать этот код на чистый C. Соответственно, буду признателен за любые советы по оптимизации.

Кстати, что в данном случае даёт использование vector<double>xxx(M), нельзя-ли обойтись обычными массивами  double xxx[M]?

[Dale]

1. Однозначно заменить массив tmp_matrix[j] функцией tmp_matrix(i, j).
2. Раз уж код планируется портировать код на C - победить лень и заменить все вызовы pow() на явное возведение в квадрат или куб.
3. Не увлекаться делением целого на целое вроде 4/3. Результат может неприятно удивить.
4. Избавиться от ненужных одномерных массивов. Например, забивать целый массив значениями x[j] = data[4 * j] не есть разумное использование памяти, эти значения вычисляются достаточно эффективно на ходу. Также, если summands[k], first_summand[k] и т.п. используются только для вычисления result[k], их незачем хранить все - вместо массива достаточно одной временной переменной.
5. Воспользовавшись фактом, что матрица tmp_matrix симметрична, т.е. tmp_matrix[i][j] == tmp_matrix[j][i], можно сократить объем вычислений почти вдвое.
6. Повыносить инварианты из циклов. Например, во фрагменте

Код: (C)

for (i = 0; i < M; i++)
{
  summands[k] += pieces[i] * pieces[j] * cos(r[j] * params[k]) * cos(r[i] * params[k]) * sin(params[k] * tmp_matrix[j][i]);
}

можно вычислить инвариантное подвыражение pieces[j] * cos(r[j] * params[k]) до входа в цикл. Аналогично и в других циклах.
Конечно, правильный оптимизирующий компилятор должен в принципе выполнить эту оптимизацию самостоятельно, но на компилятор надейся...

[Dale]

P.S. Еще немаловажный момент совсем позабыл:
7. Раскидать вычисления на несколько потоков. Сегодня вряд ли можно найти процессор с менее чем двумя ядрами (в телефонах их и то по четыре), а в исходном варианте кода задействовано будет лишь одно из них.

[Nick]

2darkelf,  Этот код вызывался в scipy.weave, vector было удобнее инициализировать переменной приходящей из python.


2Dale Спасибо, за исчерпывающий ответ  .

Все указания ясны, кроме  пункта 5 при организации функции tmp_matrix().

Такой вопрос, в целом насколько адекватно вычисление result[], т.е. я имею ввиду не имеет ли смысла поменять циклы по j,k,i местами и "перекомпановать"  их тела?

[Dale]

Все указания ясны, кроме  пункта 5 при организации функции tmp_matrix().

Посмотрите внимательно на формулу для вычисления элемента матрицы:

Код: (C)

tmp_matrix[j][i] = sqrt(pow(x[j]-x[i], 2) + pow(y[j]-y[i], 2) + pow(z[j]-z[i], 2)) / M_PI;

Очевидно, что матрица симметрична относительно главной диагонали, т.е. tmp_matrix[j][i] == tmp_matrix[i][j].

Этим можно воспользоваться. Например, если уже вычислено значение tmp_matrix[2][3], можно не считать tmp_matrix[3][2] - будет то же самое. Таким образом, можно не считать элементы выше (или ниже, как больше нравится) главной диагонали матрицы.

[Dale]

Такой вопрос, в целом насколько адекватно вычисление result[]

Довольно грязная арифметика, ее долго можно причесывать до приличного вида. Вот пример (фрагмент выражения):

Код:

pow(A ,2) * pow(global[k], 2) - 2 * A * first_summand[k] * global[k];

Для ясности обозначу global[k] просто g, first_summand[k] - f. Перепишем выражение:

A²*g² - 2*A*f*g

Считаем операции:
(A*A*g*g) - (2*A*f*g)
Итого 6 умножений и одно вычитание.

А теперь вынесем за скобки A*g:

(A*g)*(A*g - 2*f)

Произведение A*g можно вычислить один раз и сохранить во временной переменной; тогда останется еще 2 умножения и одно вычитание. Одним махом ускорили вычисление вдвое - вместо 6 умножений осталось лишь 3.

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

[Nick]

Ok, в целом понятно.

Еще раз спасибо, будем работать над арифметикой;).

[Вад]

Я бы ещё косинусы/синусы предвычислял: они там M раз для каждой позиции во вложенном цикле повторно вычисляются (я про cos( r[j] * params[k] ) и т.п.).
А ещё лучше так: всё, что зависит только от внешних параметров, считать до внутреннего цикла (по i), т.е. значения выражений

Код: (C++)

pieces[j] * cos(r[j] * params[k])

занести заранее в массив - для параметра i есть ровно такой же множитель  pieces[i] * cos( r[i] * params[k] ), его из той же таблицы и брать.

Кроме того, умножение на этот член внутри цикла по i можно не делать (вообще, предполагая произвольный диапазон x/y/z, даже и лучше не делать - точность вычислений, скорее всего, куда выше будет), а сделать его в конце, после накопления суммы: он же для всех слагаемых одинаковый.