如何加快这个用于选择子矩阵的 GSL 代码？

Question

我在 GSL 中写了一个非常简单的函数，select 来自结构中现有矩阵的子矩阵。

编辑：我的计时非常不正确，没有注意到 front.Still 中零的数量发生了变化，我希望这可以加快速度

对于 10000x10000 矩阵的 100x100 子矩阵，需要 1.2E-5 秒。因此，重复 1E4 次，比我对角化 100x100 矩阵所需的时间长 50 倍。编辑：我意识到，即使我注释掉 return(0) 以外的所有内容，它也会发生；因此，我推测，它一定是关于 struct TOWER 的。这是 TOWER 的样子：

struct TOWER 
{
    int array_level[TOWERSIZE];
    int array_window[TOWERSIZE];
    gsl_matrix *matrix_ordered_covariance;
    gsl_matrix *matrix_peano_covariance;

    double array_angle_tw[XISTEP];
    double array_correl_tw[XISTEP]; 
    gsl_interp_accel *acc_correl;   // interpolating for correlation
    gsl_spline *spline_correl;

    double array_all_eigenvalues[TOWERSIZE]; //contains all eiv. of whole matrix

    std::vector< std::vector<double> > cropped_peano_covariance, peano_mask;

};

下面是我的函数！

/* --- --- */
int monolevelsubmatrix(int i, int j, struct TOWER *tower, gsl_matrix *result)  //relying on spline!! //must addd auto vanishing
{
    int firstrow, firstcol,mu,nu,a,b;
    double aux, correl;

    firstrow = helix*i;
    firstcol = helix*j;

    gsl_matrix_view Xi = gsl_matrix_submatrix (tower ->matrix_ordered_covariance, firstrow, firstcol, helix, helix);
    gsl_matrix_memcpy (result, &(Xi.matrix));

    return(0);  
}
/* --- --- */

Answer 1

问题几乎可以肯定是gls_matric_memcpy。其来源在 copy_source.c 中，其中：

    const size_t src_tda = src->tda ;
    const size_t dest_tda = dest->tda ;
    size_t i, j;

    for (i = 0; i < src_size1 ; i++)
      {
        for (j = 0; j < MULTIPLICITY * src_size2; j++)
          {
            dest->data[MULTIPLICITY * dest_tda * i + j] 
              = src->data[MULTIPLICITY * src_tda * i + j];
          }
      }

这会很慢。请注意 gls_matrix_memcpy returns a GLS_ERROR 如果矩阵大小不同，那么很可能可以在 dest 和 src 的数据成员上使用 CRT memcpy 为数据成员提供服务。

这个循环很慢。每个单元格都通过数据成员的 dest 和 src 结构取消引用，然后进行索引。

您可以选择编写库的替代品，或者编写您自己的此矩阵副本的个人版本，例如（此处未测试的建议代码）：

unsigned int cellsize = sizeof( src->data[0] ); // just psuedocode here

memcpy( dest->data, src->data, cellsize * src_size1 * src_size2 * MULTIPLICITY )

请注意 MULTIPLICITY 是一个定义，通常是 1 或 2，可能取决于库配置 - 可能不适用于您的使用（如果它是 1 ）

现在，重要的警告....如果源矩阵是一个子视图，那么你必须按行...也就是说，在 i 中的行循环，其中 crt 的 memcpy 一次仅限于行，而不是我上面显示的整个矩阵。

换句话说，您确实必须考虑从中获取子视图的源矩阵几何形状...这可能就是为什么它们为每个单元格编制索引（使其变得简单）。

但是，如果您了解几何结构，则很可能可以通过这种方式优化您所看到的性能。

如果您所做的只是去掉 src/dest 取消引用，您会看到一些性能提升，如：

        const size_t src_tda = src->tda ;
        const size_t dest_tda = dest->tda ;
        size_t i, j;

        float * dest_data = dest->data; // psuedocode here
        float * src_data  = src->data; // psuedocode here

        for (i = 0; i < src_size1 ; i++)
          {
            for (j = 0; j < MULTIPLICITY * src_size2; j++)
              {
                dest_data[MULTIPLICITY * dest_tda * i + j] 
                  = src_data[MULTIPLICITY * src_tda * i + j];
              }
          }

我们希望编译器无论如何都能识别出这一点，但是...有时...

如何加快这个用于选择子矩阵的 GSL 代码？

How to speed up this GSL code for selecting a submatrix?

c++

gsl