Geek asks and answers

Boost :: multi_array - link too slow

I need to pass arrays to other functions using a link or pointer, I don't care as long as it runs fast. This is why I started using the boost library. I did it as follows:

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;

as functions:

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)

and in the main function:

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

Obviously, I am doing something wrong, because the code works much slower than the static version, which does not contain references / pointers, just passing arrays from one function to another.

What can I do to speed this up?

Thanks for any help ..

edit: I send what these codes do, a sequence from the GMRES solver: all arrays in it were also initialized using Boost, for example:

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];
            }
+3
source share
2

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]


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

, , , multi_array.

, , 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.

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, .

+4

. , Blaze DynamicMatrix 2 .

0

More articles:


All Articles