我有一个在 32x32x32 float32 矩阵上运行的 7 点 3d 模板操作。矩阵使用运行之间“相同”的一些值进行初始化,因此不存在随机分量。 下面的原始实现使用深度复制方法,每次时间步迭代都会将 B 矩阵中的值复制到 A 矩阵,这样我们最终会在每次迭代之间得到 2 个相同的矩阵。注意:
arr_t
只是
***float
的类型定义。void stencil_3d_7point(arr_t A, arr_t B, const int nx, const int ny, const int nz)
{
int i, j, k;
for (int timestep = 0; timestep < MAX_TIME; ++timestep)
{
for (i = 1; i < nx - 1; i++)
for (j = 1; j < ny - 1; j++)
for (k = 1; k < nz - 1; k++)
B[i][j][k] = (A[i - 1][j][k] +
A[i][j - 1][k] +
A[i][j][k - 1] +
A[i][j][k] +
A[i][j][k + 1] +
A[i][j + 1][k] +
A[i + 1][j][k]) /
7.0;
for (i = 1; i < nx - 1; i++)
for (j = 1; j < ny - 1; j++)
for (k = 1; k < nz - 1; k++)
A[i][j][k] = B[i][j][k]; // Why not just swap the pointers?
}
我的想法是在迭代之间交换指针 A 和 B,以便我们覆盖 A 中我们不再感兴趣的旧值,同时读取和操作 B 中的新值。请参阅下面的实现。
for (int timestep = 0; timestep < MAX_TIME; ++timestep)
{
for (i = 1; i < nx - 1; i++)
for (j = 1; j < ny - 1; j++)
for (k = 1; k < nz - 1; k++)
B[i][j][k] = (A[i - 1][j][k] +
A[i][j - 1][k] +
A[i][j][k - 1] +
A[i][j][k] +
A[i][j][k + 1] +
A[i][j + 1][k] +
A[i + 1][j][k]) /
7.0;
// Swap A and B
float ***temp = A;
A = B;
B = temp;
}
然而,这两种方法会产生不同的结果。我通过打印 3d 矩阵中间的 10 个值来检查这一点。
这是我初始化并调用模板函数的主要功能的一部分:
int sz = 32;
int sx = 32;
int sy = 32;
float ***A = (float ***)malloc(sx * sizeof(float **));
for (i = 0; i < sy; i++)
{
A[i] = (float **)malloc(sy * sizeof(float *));
for (j = 0; j < sz; j++)
{
A[i][j] = (float *)malloc(sz * sizeof(float));
}
}
float ***B = (float ***)malloc(sx * sizeof(float **));
for (i = 0; i < sy; i++)
{
B[i] = (float **)malloc(sy * sizeof(float *));
for (j = 0; j < sz; j++)
{
B[i][j] = (float *)malloc(sz * sizeof(float));
}
}
// initialize the A matrix
initValues(A, sx, sy, sz, INNER, OUTER);
stencil_3d_7point(A, B, sz, sy, sz);
initValues
函数用-1和1值填充矩阵,这在执行之间是相同的。同样,初始值没有随机成分。
深层复制版本每次执行都会打印此内容: 1.000000 0.126466 -0.502245 -0.830534 -0.956082 -0.991445 -0.998770 -0.999873 -0.999991 -0.999999指针交换版本每次执行都会打印此内容: 1.000000 -0.192474 -0.617592 -0.869564 -0.965497 -0.993137 -0.998986 -0.999892 -0.999991 -0.999999
如您所见,这些并不相同。
这样做的全部目的是避免昂贵的深度复制操作以节省执行周期。
使用的编译器和标志如下:
CC = aarch64-linux-gnu-gcc
CFLAGS = -O2 -static -march=armv8-a+sve
我编译为armv8以在gem5模拟器中运行它。
我尝试在 A 和 B 中打印相同的位置来检查实际结果是否位于另一个矩阵中,但这也与深度复制结果不同。
我还尝试在模板函数末尾打印行,而不是在主函数中调用之后打印行,但这会打印相同的数据。我只是尝试了一下,看看范围是否有效果。
这就是矩阵的初始化方式:
void setElement(float ***array, float value, int x, int y, int z)
{
array[x][y][z] = value;
}
void initValues(float ***array, int sx, int sy, int sz, float inner_temp, float outer_temp)
{
for (int i = 1; i < (sx - 1); i++)
{
for (int j = 1; j < (sy - 1); j++)
{
for (int k = 1; k < (sz - 1); k++)
{
setElement(array, inner_temp, i, j, k);
}
}
}
for (int j = 0; j < sy; j++)
{
for (int k = 0; k < sz; k++)
{
setElement(array, outer_temp, 0, j, k);
setElement(array, outer_temp, sx - 1, j, k);
}
}
for (int i = 0; i < sx; i++)
{
for (int k = 0; k < sz; k++)
{
setElement(array, outer_temp, i, 0, k);
setElement(array, outer_temp, i, sy - 1, k);
}
}
for (int i = 0; i < sx; i++)
{
for (int j = 0; j < sy; j++)
{
setElement(array, outer_temp, i, j, 0);
setElement(array, outer_temp, i, j, sz - 1);
}
}
}
outer_temp
和
inner_temp
分别为 -1 和 1。以下是我打印值的方法:
// Print 10 first elements of the result
for (i = 9; i < 10; i++)
{
for (j = 9; j < 10; j++)
{
for (k = 0; k < 10; k++)
{
printf("%f ", A[i][j][k]);
}
printf("\n");
for (k = 0; k < 10; k++) {
printf("%f ", B[i][j][k]);
}
printf("\n");
}
printf("\n");
}
在分配局部变量时将引用传递给原始指针。,
void stencil_3d_7point(arr_t *pA, arr_t *pB, const int nx, const int ny, const int nz)
{ arr_t A = *pA;
arr_t B = *pB;
/* ... */
*pA = B;
*pB = A;
}