如何改进子矩阵-向量乘法c代码

Question

我有用于计算通用子矩阵向量（gesubmv）乘法的 C 代码 y(rinds) = A(rinds, cinds)*x(cinds) + y(rinds)，其中 rinds 和 cinds 分别是行和列（不一定排序）索引的向量。对于通用性，rinds = rindsA_data(rindsA_1:rindsA_m) 和 cindsA_data(cindsA_1:cindsA_n).

我的代码如下：

void gesubmv(const double A_data[], const int A_size[2],
             const int rindsA_data[], int rindsA_1, int rindsA_m,
             const int cindsA_data[], int cindsA_1, int cindsA_n,
             const double x_data[], double y_data[]) {
    int i;
    int j;
    for (i = rindsA_1; i <= rindsA_m; i++) {
        int k;
        double tmp_sum;
        tmp_sum = 0.0F;
        for (j = cindsA_1; j <= cindsA_n; j++) {
            k = cindsA_data[j - 1] - 1;
            tmp_sum +=
                x_data[k] * A_data[k + (A_size[1] * (rindsA_data[i - 1] - 1))];
        }
        k = rindsA_data[i - 1] - 1;
        y_data[k] += tmp_sum;
    }
}

问题是计算 y(:) = A(:, cinds)*x(cinds) + y(:) 然后选择 y(rinds) 似乎更有效。我怀疑这可能与指令的执行方式有关 (https://godbolt.org/)。或者我可以做一些（更明显的）错误。有什么想法吗？

备选方案 y(:) = A(:, cinds)*x(cinds) + y(:) 的代码如下。

void gesubmv2(const double A_data[], const int A_size[2],
             const int cindsA_data[], int cindsA_1, int cindsA_n,
             const double x_data[], double y_data[]) {
    int i;
    int j;
    for (i = 0; i < A_size[1]; i++) {
        int k;
        double tmp_sum;
        tmp_sum = 0.0F;
        for (j = cindsA_1; j <= cindsA_n; j++) {
            k = cindsA_data[j - 1] - 1;
            tmp_sum += x_data[k] * A_data[k + (A_size[1] * i)];
        }
        y_data[i] += tmp_sum;
    }
}

Answer 1

您可以将 (A_size[1] * (rindsA_data[i - 1] - 1)) 表达式从 j 循环中取出。优化器可能已经在执行此操作，因此您可能看不到任何性能差异。同样，在第二个片段中，(A_size[1] * i) 不需要位于 j 循环内。如果您的优化器没有按照我的建议进行，这可能就是您看到性能差异的原因。

此外，您使用 for 语句递增 j，然后使用 j-1 作为索引，因此您可以通过在循环之前和计算 tmp_sum 之后分别放置一次 k assignmet 来节省一个操作。例如

    kk = (A_size[1] * (rindsA_data[i - 1] - 1));
    k = cindsA_data[cindsA_1 - 1] - 1;
    for (j = cindsA_1; j < cindsA_n; j++) {
        tmp_sum +=
            x_data[k] * A_data[k + kk];
        k = cindsA_data[j] - 1;
    }
    tmp_sum += x_data[k] * A_data[k + kk];