Boost::multi_array -- 引用太慢

Question

我必须通过引用或指针将数组传递给其他函数，只要它工作得快，我不在乎。这就是我开始使用 boost 库的原因。我是通过以下方式做到的：

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);

显然，我做错了，因为代码的运行速度比静态版本慢，静态版本不涉及任何引用/指针，只是将数组从一个函数传递到另一个函数。

我能做些什么来加速这个过程？

感谢您的任何帮助..

编辑：我正在发布这些代码的作用，来自 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];
            }

Answer 1

在您对 multi_aray 进行零初始化的地方，您可以尝试改用 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]

试试

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，您基本上是以速度换取灵活性和便利性。

Answer 2

见rodrigob's answer here。此外，使用具有相同编译器优化的 Blaze DynamicMatrix 可以提供几乎额外的因子 2 改进。