基于OpenMP实现无额外内存分配的3D点云变换并行化

我正在开发一款C程序，使用3×3变换矩阵T对3D点云执行线性变换，每个点对应n×3矩阵中的一行。我的目标是就地更新点云，但在并行更新过程中遇到了数据依赖问题，且因内存限制无法为结果分配额外矩阵。

原始代码 #

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char *argv[]) {
    int num_of_threads = atoi(argv[1]);
    int n = atoi(argv[2]);
    int points[n][3], T[3][3], i, j, k;
    srand(time(NULL));

    omp_set_num_threads(num_of_threads);
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            T[i][j] = rand() % 100;
        }
    }

    // Initialize points
    #pragma omp parallel for collapse(2)
    for (i = 0; i < n; i++) {
        for (j = 0; j < 3; j++) {
            points[i][j] = rand() % 1000;
        }
    }

    // Apply transformation
    #pragma omp parallel for private(i, j, k, temp)
    for (i = 0; i < n; i++) {
        for (j = 0; j < 3; j++) {
            int temp = 0;
            for (k = 0; k < 3; k++) {
                temp += points[i][k] * T[j][k];
            }
            points[i][j] = temp; // Direct update introduces data dependency issues
        }
    }

    return 0;
}

面临挑战 #

• 直接更新points[i][j]会引发数据依赖问题，因为不同坐标的变换并非相互独立。
• 使用临界区会严重影响性能，且无法解决读写同一位置的问题。

问题 #

是否存在更优方法，可在不分配额外内存存储变换后点云的前提下，利用OpenMP高效处理依赖并实现该操作的并行化？

修改后的代码 #

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double my_rand(long *seed)
{
    long next_rnd = (1366 * (*seed) + 150889) % 714025;
    *seed = next_rnd;
    return ((double)next_rnd / (double)714025);
}

int main(int argc, char *argv[])
{

    int num_of_threads = atoi(argv[1]);
    int n = atoi(argv[2]);
    int T[3][3], i, j, k;
    double points[n][3];
    long seed;
    srand(time(NULL));

    omp_set_num_threads(num_of_threads);
    // matrix T
    for (i = 0; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
        {
            T[i][j] = rand() % 100;
        }
    }

    // points
#pragma omp parallel for collapse(2) private(i, j, seed) shared(points,n)
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < 3; j++)
        {
            seed = omp_get_thread_num() + 123456789;
            points[i][j] = my_rand(&seed);
        }
    }

//  transformation
#pragma omp parallel for private(i, j, k) shared(points, n)
    for (i = 0; i < n; i++)
    {
        double temp[3] = {0, 0, 0}; // temporary local storage for transformed coordinates; private by default
        for (j = 0; j < 3; j++)
        {
            for (k = 0; k < 3; k++)
            {
                temp[j] += points[i][k] * T[j][k];
            }
        }
        for (j = 0; j < 3; j++)
        {
            points[i][j] = temp[j]; // write transformed points back to points array
        }
    }

    return 0;
}

内容的提问来源于stack exchange，提问作者untitled43