c++ - C++如何在内存中连续(连贯)存储三维数组？

Question

这是我第一次使用 SO。对不起我的英语，但我会尽力而为，并完整地描述我的大学任务问题，以便您更好地理解。我确定我的问题更多是与代码和数组相关，而不是数学。

我正在使用 MPI、OpenMP 以及大尺寸和网格大小解决立方体中的数值 3D 波问题。现在我正在编写代码以在我们大学的 super 计算机上执行它。

所以一些小的解释和公式可以更好地理解问题:

1. 我有一个立方体，每个轴上有 N+1 个点(数字从 0 到 N )
2. 让我有 s = MPI_size 排名为 的进程总数0 到 s-1
3. 我正在做 s 平行区域，制作垂直于 Y 轴的切片(具有固定的 j 数)。我的意思是索引 i 与往常一样对应于 x 轴，j 对应于 Y，k 对应于 Z。
4. 所以现在每个等级从 0 到 s-1 的进程负责我的 Y 轴上的 p=(N+1)/s 个点:从 rank*p 到 rank*p+p-1
5. 现在我在创建三个大小为 [N+1][p][N+1] 的数组时遇到了问题。三是因为在计算的每一步中，我都使用前一时间步和“双前(前一个)”的值(抱歉，我不知道怎么说对:))

我试着简单地做到这一点:

    for(int i = 0; i < N_p+1; i++) {
      u[i] = new double*[N_p+1];
      u_prev[i] = new double*[N_p+1];
      u_prev_prev[i] = new double*[N_p+1];

      for(int j = 0; j < N_p+1; j++) {
          u[i][j] = new double[N_p+1];
          u_prev[i][j] = new double[N_p+1];
          u_prev_prev[i][j] = new double[N_p+1];

          for(int k = 0; k < N_p+1; k++) { 
              u[i][j][k] = 0.0;
              u_prev[i][j][k] = 0.0;
              u_prev_prev[i][j][k] = 0.0;
          }
      }
}

但是我已经完成这项任务的 friend 告诉我，当我需要向其他进程发送消息或从其他进程接收消息时，我会遇到问题(因为我需要发送大小为 (N+1)^2 的完整层 垂直于 Y 轴。在内存中放置这样的数组将是一个大问题。

此外，我们不需要像我一开始那样在 Y 轴上设置所有 N+1 点，每个过程只需要 p 个点。因为内存问题是非常实际的(在某些测试中最大 N 将在 1536 左右。而我们用于测试的资源非常少)

所以他提议这样做:(他做了垂直于 X 轴的平行切片，而不是像我一样垂直于 Y 轴和从 0 到 N-1 的计数而不是从 0 到 N 对我来说，这不是什么大问题。但是这段代码对我来说是某种魔法，反正我也不完全理解它。)

float* buffer = new float[N * N * (N/s) * 3];
float ***u;

for (int i = left; i < right; ++i) { 
  u[i] = new float*[N]; 
  for (int j = 0; j < N; ++j) { 
    u[i][j] = buffer + (i - left)*N*N + j*N; 
  } 
}

所以当我的大脑着火并且很快就会爆炸时，我尝试做一些类比:

float* buffer = new float[(N+1) * p * (N+1) * 3];
float ***u;

for (int i = 0; i <= N; ++i) {
u[i] = new float*[p];
//p = right-left+1
  for (int j = left; j <= right; ++j) {
    u[i][j] = buffer + ???; //SOS
  }
 }

请有人尝试理解这种方法并解释放置什么而不是“???”。或者任何其他更好的解决方案。

我的队友也告诉我可以在代码中简单地写 u[i][j][k] 但是我不确定我是否会有 k 具有此类定义的索引以及我是否需要更改所有计算语法。

很抱歉这么大的解释和问题。真不想自己弄明白，赶紧解决。但现在我卡住了。

还粘贴了一个立方体的小图像，并在其中用红色标记了图层以进行一些简单的可视化。

Answer 1

我更喜欢将此类问题封装在类中(OOP 对吗？)。它更容易理解、可重用并且开销很小。任何类型的 3D 矩阵都是:

template<typename T>
class Matrix {
private:
    int _dimX, _dimY, _dimZ;
    T *_storage;
public:
    inline Matrix(int dimX, int dimY, int dimZ) {
        _dimX = dimX; _dimY = dimY; _dimZ = dimZ;
        _storage = new T[_dimX * _dimY * _dimZ];
    }
    inline ~Matrix() { delete _storage; }
    inline T* getStorage() const { return _storage; }
    inline T &operator()(int x, int y, int z) {
        return _storage[x * _dimX * _dimY + y * _dimY + z];
    }
};

并使用它:

#include <iostream>
using std::cout;
using std::endl;

int main() {
    int xlen = 10, ylen = 5, zlen = 2;
    // Creating matrix of doubles
    Matrix<double> matrix(xlen, ylen, zlen);

    // Filling using contiguous memory array (to export)
    double *d = matrix.getStorage();
    int count = 0;
    for (int i = 0; i < xlen * ylen * zlen; i++)
        *d++ = ++count; // same as d[i] = ++count;

    // Using with indexes
    cout << "matrix(1,0,0) = " << matrix(1, 0, 0) << endl; // prints 51
    matrix(1, 0, 0) = 34;
    cout << "matrix(1,0,0) = " << matrix(1, 0, 0) << endl; // prints 34
    matrix(1, 0, 0)++;
    cout << "matrix(1,0,0) = " << matrix(1, 0, 0) << endl; // prints 35

    return 0;
}