犰狳中的快速阵列置换（广义张量转置）（C）

我有一个项目涉及3D阵列的大量排列（ arma::Cube<cx_double> ） . 特别地，所需的排列是通过切片交换列 . 在Matlab中，这是由 permute(cube,[1,3,2]) 有效计算的，在Python中由 numpy.transpose(cube,axis=[0,2,1]) 有效计算 .

不幸的是，Armadillo本身没有 permute 功能 . 我尝试了不同的方法，但与Matlab相比，它们都相当慢 . I would like to know what's the faster way to permute (rather large) Cubes in Armadillo . 使用 gprof 对代码进行分析，大部分时间用于我在下面尝试的置换函数，而在Matlab中，对于相同的移植项目，大部分时间用于SVD或QR矩阵分解（重塑和置换都很快）在matlab） .

我想了解哪种方法是在犰狳中进行这种排列的最快方法，以及为什么有些方法比其他方法效果更好 .

Option 1: Raw permutation （最快的选择）（有更快的方法吗？）

输入多维数据集的元素分配给输出多维数据集 .

template <typename T>
static Cube<T> permute (Cube<T>& cube){

uword D1=cube.n_rows;
uword D2=cube.n_cols;
uword D3=cube.n_slices;

Cube<T> output(D1,D3,D2);

for (uword s = 0; s < D3; ++s){
for (uword c = 0; c < D2; ++c){
for (uword r = 0; r < D1; ++r){
    output.at(r, s, c) = cube.at(r, c, s);
    // output[ D1*D3*c + D1*s+ r ] = cube[ D1*D2*s + D1*c + r ]; 

}
}
}

return output;
}

Option 2: Filling slices （很慢）

通过非连续的 subcube 视图填充输出多维数据集的切片 .

template <typename T>
static Cube<T> permute (Cube<T>& cube_in){

    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    Cube<T> output; 
    output.zeros(D1, D3, D2);

    for (uword c=0; c<D2; ++c) {
        output.slice(c) = cube_in.subcube( span(0,D1-1),span(c),span(0,D3-1) );
    }

    return output;
}

Option 3: Transposing layers （比原始排列慢但可比较）

我们可以迭代输入多维数据集的图层（固定行）并转置它们 .

template <typename T>
static Cube<T> permute (Cube<T>& cube_in){
    // in a cube, permute {1,3,2} (permute slices by columns)

    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    if(D3 > D2){
        cube_in.resize(D1,D3,D3);
    } else if (D2 > D3) {
        cube_in.resize(D1,D2,D2);
    }

    for (uword r=0; r<D1; ++r) {        
        static cmat layer = cmat(cube_in.rows(r,r)); 
        inplace_strans(layer);
        cube_in.rows(r,r)=layer;               
    }

    cube_in.resize(D1,D3,D2);

    return cube_in;
}

Option 4: Look-up table 通过读取向量中的indeces获取非连续访问 .

template <typename T>
arma::Cube<T> permuteCS (arma::Cube<T> cube_in){
    // in a cube, permute {1,3,2} (permute slices by columns)
    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    cx_vec onedcube = cube_in.elem(gen_trans_idx(cube_in));            

    return arma::Cube<T>(onedcube.memptr(), D1, D3, D2, true ) ;
}

其中 gen_trans_idx 是一个生成置换立方体的indeces的函数：

template <typename T>
uvec gen_trans_idx(Cube<T>& cube){

    uword D1 = cube.n_rows;
    uword D2 = cube.n_cols;
    uword D3 = cube.n_slices;

    uvec perm132(D1*D2*D3);    
    uword ii = 0;                
    for (int c = 0; c < D2; ++c){
    for (int s = 0; s < D3; ++s){
    for (int r = 0; r < D1; ++r){
        perm132.at(ii) = sub2ind(size(cube), r, c, s);
        ii=ii+1;
    }}} 

    return perm132;
}

理想情况下，如果预先确定立方体尺寸，则可以预先计算这些查找表 .

Option 5 （就地换位）非常慢，内存效率高

// Option: In-place transpose
template <typename T>
arma::Cube<T> permuteCS (arma::Cube<T> cube_in, uvec permlist ){    

    T* Qpoint = cube_in.memptr(); // pointer to first element of cube_in

    uvec updateidx = find(permlist - arma::linspace<uvec>(0,cube_in.n_elem-1,cube_in.n_elem)); // index of elements that change position in memory           
    uvec skiplist(updateidx.n_elem,fill::zeros);

    uword rr = 0; // aux index for updatelix
    for(uword jj=0;jj<updateidx.n_elem;++jj){                

        if(any(updateidx[jj] == skiplist)){ // if element jj has already been updated
            // do nothing
        } else {

            uword scope = updateidx[jj];
            T target = *(Qpoint+permlist[scope]);        // store the value of the target element             

            while(any(scope==skiplist)-1){              // while wareyou has not been updated

                T local  = *(Qpoint+scope);                  // store local value                                                                
                *(Qpoint+scope) = target;       
                skiplist[rr]=scope;
                ++rr;            
                uvec wareyou = find(permlist==scope);        // find where the local value will appear                
                scope = wareyou[0];
                target = local;                            
            }

        }
    }
   cube_in.reshape(cube_in.n_rows,cube_in.n_slices,cube_in.n_cols);

    return cub

e_in;
}