Boost :: multi_array - ссылка слишком медленная - PullRequest
Новая
       17

Boost :: multi_array - ссылка слишком медленная

3 голосов
/ 12 марта 2011

Я должен передать массивы в другие функции через ссылки или указатели, мне все равно, пока это работает быстро. Вот почему я начал использовать библиотеку наддува. Я сделал это следующим образом: 

using namespace boost;

typedef  multi_array<long double, 4> array_type;
typedef  multi_array<long double, 2> twod_array_type;
typedef  multi_array<long double, 1> vec_type;

как функции: 

void pde_3d_7_stencil_discretization(array_type& A, vec_type& b, vec_type& x,const int& xdim, const int& ydim,const int& zdim)

void gmressolver3d(array_type& A, vec_type& x, vec_type& rhs,const int& KrylovDim,const int& xdim,const int& ydim,const int& zdim,const int& COP, const int& threeDStencil)

и в основной функции: 

  array_type A(extents[threeDimStencil][COP][COP][xdim*ydim*zdim]);
  vec_type b(extents[xdim*ydim*zdim*COP]);
  vec_type x(extents[xdim*ydim*zdim*COP]);

  pde_3d_7_stencil_discretization(A,b,x,xdim,ydim,zdim);
  gmressolver3d(A,x,b,KrylovDim,xdim,ydim,zdim,COP,threeDimStencil);

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

Что я могу сделать, чтобы ускорить это?

Спасибо за любую помощь ..

edit: я публикую то, что делают эти коды, последовательность из решателя GMRES: все массивы в нем были инициализированы также с использованием Boost, например: 

vec_type pp(extents[zdim*xdim*ydim*COP]);
vec_type ppp(extents[zdim*xdim*ydim*COP]);
vec_type w(extents[zdim*xdim*ydim*COP]);
vec_type y(extents[KrylovDim]);
vec_type vv(extents[zdim*xdim*ydim*COP]);
vec_type b(extents[KrylovDim+1]);
vec_type ro(extents[zdim*xdim*ydim*COP]);
vec_type out1(extents[xdim*zdim*ydim*COP]);
vec_type m_jac(extents[xdim*zdim*ydim*COP]);
twod_array_type h(extents[KrylovDim+1][KrylovDim]);
twod_array_type v(extents[zdim*xdim*ydim*COP][KrylovDim]);
twod_array_type hess(extents[KrylovDim+1][KrylovDim]);
array_type maa(extents[threeDStencil][COP][COP][zdim*xdim*ydim]);
array_type maaa(extents[threeDStencil][COP][COP][zdim*xdim*ydim]);

for (i=0;i<m+1;i++){
            b[i] = 0;
            for(k=0;k<m;k++){
                h[i][k] = 0.0;
            }
        }

        for (i=0;i<n;i++){
            v[i][0] = ro[i]/r;
        }
        for(j=0;j<m;j++){
            b[0] = r;
            vector_zero_fill(n,ppp);
            for(i=0;i<n;i++){
                vv[i]=v[i][j];
            }
            //********************MATRIX FREE********************
            matrix_vector_product_heptadiagonal_discret(A,vv,pp,xdim,ydim,zdim);
            //two_vector_dot_product(n,pp,m_jac);
    //      if(isPrec)
    //      forback(A,pp);
            //********************MATRIX FREE********************
            //pretty fast**
            for(i=0;i<=j;i++){
                for(k=0;k<n;k++){
                    h[i][j] = h[i][j] + pp[k]*v[k][i];
                }
            }

            for(i=0;i<=j;i++){
                for(k=0;k<n;k++){
                    ppp[k] = ppp[k] + h[i][j]*v[k][i];
                }
            }
            p=0.0;

            for(i=0;i<n;i++){
                w[i] = pp[i] - ppp[i];
                p = p + pow(w[i],2);
            }

            h[j+1][j] = sqrt(p);

            for(i=0;i<=j+1;i++){
                for(k=0;k<=j;k++){
                    hess[i][k] = h[i][k];
                }
            }
            for(i=0;i<j+1;i++){
                c = hess[i][i]/sqrt(pow(hess[i][i],2)+pow(hess[i+1][i],2));
                s = hess[i+1][i]/sqrt(pow(hess[i][i],2)+pow(hess[i+1][i],2));
                for (k=0;k<=j;k++){
                    inner1=c*hess[i][k]+s*hess[i+1][k];
                    inner2=(-s)*hess[i][k]+c*hess[i+1][k];
                    hess[i][k] = inner1;
                    hess[i+1][k] = inner2;
                }
                b[i+1] = -s*b[i];
                b[i] = c*b[i];
            }

Ответы [ 2 ]

4 голосов
/ 12 марта 2011

Там, где вы инициализируете ваши multi_arays с нуля, вы можете попробовать использовать std::memset.Например, 

std::memset(b.data(), 0, size_of_b_in_bytes);

В вашем коде есть несколько мест, где вы можете индексировать один и тот же элемент multi_array более одного раза.Например, вместо 

h[i][j] = h[i][j] + pp[k]*v[k][i]

try 

h[i][j] += pp[k]*v[k][i]

Обычно оптимизатор автоматически делает такие замены для вас, но, возможно, это невозможно с multi_array.

Я также заметил два цикла for, которые можно объединить в один, чтобы избежать многократного индексирования одного и того же элемента multi_array: 

/*
for(i=0; i<=j; i++)
{
    for(k=0; k<n; k++)
    {
        h[i][j] = h[i][j] + pp[k]*v[k][i];
    }
}

for(i=0; i<=j; i++)
{
    for(k=0; k<n; k++)
    {
        ppp[k] = ppp[k] + h[i][j]*v[k][i];
    }
}
*/

for(i=0; i<=j; i++)
{
    for(k=0; k<n; k++)
    {
        long double& h_elem = h[i][j];
        long double v_elem = v[k][i];
        h_elem += pp[k]*v_elem;
        ppp[k] += h_elem*v_elem;
    }
}

Возможно, таких будет больше.Обратите внимание на использование ссылок и переменных для «запоминания» элемента и избежания необходимости пересчитывать его позицию в multi_array.

В последнем цикле for из вашего кода вы можете избежать многократного пересчета индексов multi_array с помощьюиспользуя временные переменные и ссылки: 

/*
for(i=0;i<j+1;i++){
    c = hess[i][i]/sqrt(pow(hess[i][i],2)+pow(hess[i+1][i],2));
    s = hess[i+1][i]/sqrt(pow(hess[i][i],2)+pow(hess[i+1][i],2));
    for (k=0;k<=j;k++){
        inner1=c*hess[i][k]+s*hess[i+1][k];
        inner2=(-s)*hess[i][k]+c*hess[i+1][k];
        hess[i][k] = inner1;
        hess[i+1][k] = inner2;
    }
    b[i+1] = -s*b[i];
    b[i] = c*b[i];
}
*/

for(i=0;i<j+1;i++){
    long double hess_i_i = hess[i][i];
    long double hess_ip1_i = hess[i+1][i];
    long double temp = sqrt(pow(hess_i_i,2)+pow(hess_ip1_i,2));
    c = hess_i_i/temp;
    s = hess_ip1_i/temp;
    for (k=0;k<=j;k++){
        long double& hess_i_k = hess[i][k];
        long double& hess_ip1_k = hess[i+1][k];
        inner1=c*hess_i_k+s*hess_ip1_k;
        inner2=(-s)*hess_i_k+c*hess_ip1_k;
        hess_i_k = inner1;
        hess_ip1_k = inner2;
    }
    long double b_i& = b[i];
    b[i+1] = -s*b_i;
    b_i = c*b_i;
}

Дважды проверьте мою работу - наверняка я где-то допустил ошибку.Обратите внимание, что я сохранил sqrt(pow(hess_i_i,2)+pow(hess_ip1_i,2)) в переменной, чтобы он не вычислялся без необходимости дважды.

Я сомневаюсь, что эти незначительные изменения приведут к уменьшению времени выполнения до 5 секунд.Проблема с multi_array заключается в том, что размеры массива известны только во время выполнения.Поддержка упорядочения по главной / основной колонке, вероятно, также вызывает некоторые накладные расходы.

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

Используя Boost multi_arrays, вы в основном теряете скорость ради гибкости и удобства.

0 голосов
/ 19 июня 2019

См. ответ Родригоба здесь .Кроме того, использование Blaze DynamicMatrix с той же оптимизацией компилятора может дать почти дополнительное улучшение в 2 раза.

...