c - 是什么导致我的数组地址在传递给函数时被损坏(更改)？

Question

我正在执行压缩稀疏原始矩阵 vector 乘法(CSR SPMV):这涉及将数组A划分为多个 block ，然后通过引用函数传递该 block ，但是只有数组的第一部分(A[0] 从数组开头开始的第一个 block )被修改。然而，从第二个循环A[0 + chunkIndex]开始，当函数读取子数组时，它会跳转并读取超出总数组地址范围的不同地址，尽管索引是正确的。

供引用:

<小时/>

SPMV 内核是:

void serial_matvec(size_t TS,  double *A, int *JA, int *IA,  double *X, double *Y)
{
    double sum;
    for (int i = 0; i < TS; ++i)
    {   
        sum = 0.0;
        for (int j = IA[i]; j < IA[i + 1]; ++j)
        {
                sum += A[j] * X[JA[j]]; // the error is here , the function reads diffrent 
                                        // address of A,  and JA, so the access 
                                       // will be out-of-bound
            }
            Y[i] = sum;
        }
    }

它的调用方式是这样的:

int chunkIndex = 0;
for(size_t k = 0; k < rows/TS; ++k)
{
    chunkIndex = IA[k * TS];
    serial_matvec(TS, &A[chunkIndex], &JA[chunkIndex], &IA[k*TS], &X[0], &Y[k*TS]);
}

假设我处理(8x8)矩阵，并且我处理每个 block 2行，因此循环k将为r ows/TS = 4 次循环，传递给函数的 chunkIndex 和数组如下:

chunkIndex: 0 --> loop k = 0, &A[0], &JA[0]

chunkIndex: --> loop k = 1, &A[16], &JA[16] //[ERROR here, function reads different address]

chunkIndex: --> loop k = 2, &A[32], &JA[32] //[ERROR here, function reads different address]

chunkIndex: --> loop k = 3, &A[48], &JA[48] //[ERROR here, function reads different address]

当我运行代码时，只有第一个 block 正确执行，其他 3 个 block 内存已损坏，并且数组指针跳入超出数组大小的边界。

我已经手动检查了所有参数的所有索引，它们都是正确的，但是当我打印地址时它们并不相同。 (已经调试了3天了)

我使用了valgrind，它报告:

Invalid read of size 8 and Use of uninitialised value of size 8 at the sum += A[j] * X[JA[j]]; line

我用-g -fsanitize=address编译它，我得到了

heap-buffer-overflow

我尝试在函数之外手动访问这些 block ，它们是正确的，那么什么会导致堆内存像这样被损坏？

代码是here ，这是我起码能做到的。

Answer 1

问题是我在索引传递给函数的数组部分 (chunk) 时使用了全局索引(main 内的索引)，因此输出 -越界问题。

解决方案是在每次函数调用时从 0 开始索引子数组，但我遇到了另一个问题。在每次函数调用时，我都会处理 TS 行，每行都有不同数量的非零。

举个例子，看图，chunk 1，抱歉我的笔迹不好，这样比较容易。正如您所看到的，我们需要 3 个索引，其中一个索引用于每个 block TS 行i ，另一个是因为每行有不同数量的非零 j ，第三个用于索引子数组传递了 l，这是最初的问题。

并且serial_matvec函数如下:

void serial_matvec(size_t TS, const double *A, const int *JA, const int *IA,
                   const double *X, double *Y) {
  int  l = 0;
  for (int i = 0; i < TS; ++i) {
    for (int j = 0; j < (IA[i + 1] - IA[i]); ++j) {
      Y[i] += A[l] * X[JA[l]];
      l++;
    }    
  }
}

经过测试的完整代码是here如果有人有更优雅的解决方案，我们非常欢迎。