Fuzzy Unsupervised C均值算法OpenMP并行化性能不升反降问题求助

我用OpenMP实现了模糊无监督C均值（FCM）算法的并行化，但发现当线程数从4/8提升至16/32时，性能不仅没提升反而下降，且线程数并未超过CPU核心数。我确定这是**竞争条件（Race Conditions）**导致的，但不知道该如何修复。以下是我的并行化代码：

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

#define MAX_DATA_POINTS 10000
#define MAX_CLUSTER 100
#define MAX_DATA_DIMENSION 5

int num_data_points;
int num_clusters;
int num_dimensions;
double low_high[MAX_DATA_DIMENSION][2];
double degree_of_memb[MAX_DATA_POINTS][MAX_CLUSTER];
double epsilon;
double fuzziness;
double data_point[MAX_DATA_POINTS][MAX_DATA_DIMENSION];
double cluster_centre[MAX_CLUSTER][MAX_DATA_DIMENSION];
double norms[MAX_DATA_POINTS][MAX_CLUSTER]; // Precomputed norms

// Initialize data and membership matrix
int init(char *fname) {
    int i, j, r, rval;
    FILE *f;
    double s;
    if ((f = fopen(fname, "r")) == NULL) {
        printf("Failed to open input file.");
        return -1;
    }
    fscanf(f, "%d %d %d", &num_data_points, &num_clusters, &num_dimensions);
    if (num_clusters > MAX_CLUSTER || num_data_points > MAX_DATA_POINTS || num_dimensions > MAX_DATA_DIMENSION) {
        printf("Input data exceeds defined limits.\n");
        fclose(f);
        exit(1);
    }

    fscanf(f, "%lf %lf", &fuzziness, &epsilon);
    if (fuzziness <= 1.0 || epsilon <= 0.0 || epsilon > 1.0) {
        printf("Invalid fuzziness or epsilon.\n");
        fclose(f);
        exit(1);
    }

    // Initialize data points and their range
    for (i = 0; i < num_dimensions; i++) {
        low_high[i][0] = __DBL_MAX__;
        low_high[i][1] = -__DBL_MAX__;
    }

    for (i = 0; i < num_data_points; i++) {
        for (j = 0; j < num_dimensions; j++) {
            fscanf(f, "%lf", &data_point[i][j]);
            if (data_point[i][j] < low_high[j][0])
                low_high[j][0] = data_point[i][j];
            if (data_point[i][j] > low_high[j][1])
                low_high[j][1] = data_point[i][j];
        }
    }

    // Initialize membership matrix randomly
    for (i = 0; i < num_data_points; i++) {
        s = 0.0;
        r = 100;
        for (j = 1; j < num_clusters; j++) {
            rval = rand() % (r + 1);
            r -= rval;
            degree_of_memb[i][j] = rval / 100.0;
            s += degree_of_memb[i][j];
        }
        degree_of_memb[i][0] = 1.0 - s;
    }

    fclose(f);
    return 0;
}

// Precompute norms to avoid redundant calculations
void precompute_norms() {
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < num_data_points; i++) {
        for (int j = 0; j < num_clusters; j++) {
            double sum = 0.0;
            for (int k = 0; k < num_dimensions; k++) {
                double diff = data_point[i][k] - cluster_centre[j][k];
                sum += diff * diff;
            }
            norms[i][j] = sqrt(sum);
        }
    }
}

// Calculate new cluster centers
int calculate_centre_vectors() {
    double t[MAX_DATA_POINTS][MAX_CLUSTER];

    #pragma omp parallel for collapse(2)
    for (int i = 0; i < num_data_points; i++) {
        for (int j = 0; j < num_clusters; j++) {
            t[i][j] = pow(degree_of_memb[i][j], fuzziness);
        }
    }

    #pragma omp parallel for collapse(2)
    for (int j = 0; j < num_clusters; j++) {
        for (int k = 0; k < num_dimensions; k++) {
            double numerator = 0.0;
            double denominator = 0.0;
            for (int i = 0; i < num_data_points; i++) {
                numerator += t[i][j] * data_point[i][k];
                denominator += t[i][j];
            }
            cluster_centre[j][k] = numerator / denominator;
        }
    }
    return 0;
}

// Update membership values
double update_degree_of_membership() {
    double max_diff = 0.0;

    // Precompute norms
    precompute_norms();

    #pragma omp parallel for reduction(max:max_diff) collapse(2)
    for (int i = 0; i < num_data_points; i++) {
        for (int j = 0; j < num_clusters; j++) {
            double sum = 0.0;
            double norm_ij = norms[i][j];
            for (int k = 0; k < num_clusters; k++) {
                sum += pow(norm_ij / norms[i][k], 2.0 / (fuzziness - 1));
            }
            double new_uij = 1.0 / sum;
            double diff = fabs(new_uij - degree_of_memb[i][j]);
            if (diff > max_diff) {
                max_diff = diff;
            }
            degree_of_memb[i][j] = new_uij;
        }
    }
    return max_diff;
}

// FCM clustering process
int fcm(char *fname) {
    double max_diff;
    if (init(fname) != 0) return -1;
    do {
        calculate_centre_vectors();
        max_diff = update_degree_of_membership();
    } while (max_diff > epsilon);
    return 0;
}

// Print membership matrix to a file or stdout
void print_membership_matrix(char *fname) {
    int i, j;
    FILE *f;
    if (fname == NULL)
        f = stdout;
    else if ((f = fopen(fname, "w")) == NULL) {
        printf("Cannot create output file.\n");
        exit(1);
    }

    fprintf(f, "Membership matrix Parallel:\n");
    for (i = 0; i < num_data_points; i++) {
        fprintf(f, "Data[%d]: ", i);
        for (j = 0; j < num_clusters; j++) {
            fprintf(f, "%lf ", degree_of_memb[i][j]);
        }
        fprintf(f, "\n");
    }
    if (fname != NULL)
        fclose(f);
}

// Main function
int main(int argc, char **argv) {
    if (argc != 2) {
        printf("USAGE: fcm <input file>\n");
        exit(1);
    }
    double start_time = omp_get_wtime();
    fcm(argv[1]);
    double end_time = omp_get_wtime();
    double execution_time = end_time - start_time;
    printf("Number of data points: %d\n", num_data_points);
    printf("Number of clusters: %d\n", num_clusters);
    printf("Number of data-point dimensions: %d\n", num_dimensions);
    printf("Accuracy margin: %lf\n", epsilon);
    print_membership_matrix("membership.matrix");
    printf("The program took %f seconds.\n", execution_time);
    return 0;
}

你的代码里并没有明显的写竞争，但性能下降的核心原因是伪共享（False Sharing）和缓存局部性差，这两种问题会随着线程数增加被放大，表现出类似竞争条件的性能衰退。以下是具体分析和修复步骤：

核心问题分析 #

1. 全局数组的伪共享：
   所有大数组（degree_of_memb、cluster_centre等）都是连续内存存储的，当多个线程访问相邻元素时，这些元素可能落在同一个缓存行中。一个线程修改缓存行内的元素会导致其他线程的缓存失效，引发频繁的缓存同步，线程数越多，同步开销越大。

2. calculate_centre_vectors的内存与并行低效：

   • 栈上声明的大数组t[MAX_DATA_POINTS][MAX_CLUSTER]会导致栈溢出，引发未定义行为；同时该数组的访问会加剧伪共享。
   • 内层遍历数据点的循环是串行的，没有利用多线程，且每个线程重复计算同一聚类中心的累加值，浪费资源。

3. precompute_norms的缓存命中率低：
   collapse(2)后的循环同时遍历数据点和聚类中心，导致cluster_centre的访问是跳跃式的，缓存无法有效复用，内存访问开销大。

修复步骤 #

1. 解决伪共享：内存对齐与调度优化 #

对全局数组进行缓存行对齐（64字节是主流CPU的缓存行大小），避免不同线程的访问落在同一缓存行：

// 修改全局数组声明
double degree_of_memb[MAX_DATA_POINTS][MAX_CLUSTER] __attribute__((aligned(64)));
double cluster_centre[MAX_CLUSTER][MAX_DATA_DIMENSION] __attribute__((aligned(64)));
double norms[MAX_DATA_POINTS][MAX_CLUSTER] __attribute__((aligned(64)));

同时调整并行循环的调度策略，让每个线程处理连续的大块数据，提升缓存局部性：

// 在precompute_norms等并行循环中添加调度参数
#pragma omp parallel for collapse(2) schedule(static, 64)

2. 重构calculate_centre_vectors函数 #

去掉栈上的大数组t，直接在累加时计算t_ij，同时并行化内层数据点遍历，用reduction保证累加的线程安全：

int calculate_centre_vectors() {
    #pragma omp parallel for collapse(2) schedule(static)
    for (int j = 0; j < num_clusters; j++) {
        for (int k = 0; k < num_dimensions; k++) {
            double numerator = 0.0;
            double denominator = 0.0;
            // 并行遍历数据点，累加分子分母
            #pragma omp parallel for reduction(+:numerator, denominator)
            for (int i = 0; i < num_data_points; i++) {
                double t_ij = pow(degree_of_memb[i][j], fuzziness);
                numerator += t_ij * data_point[i][k];
                denominator += t_ij;
            }
            cluster_centre[j][k] = numerator / denominator;
        }
    }
    return 0;
}

3. 优化precompute_norms的缓存局部性 #

调整循环顺序，先遍历聚类中心，将当前中心的所有维度加载到缓存后再遍历数据点，提升缓存复用率：

void precompute_norms() {
    #pragma omp parallel for schedule(static)
    for (int j = 0; j < num_clusters; j++) {
        // 预加载当前聚类中心到局部变量，避免重复访问全局数组
        double centre[MAX_DATA_DIMENSION];
        for (int k = 0; k < num_dimensions; k++) {
            centre[k] = cluster_centre[j][k];
        }
        for (int i = 0; i < num_data_points; i++) {
            double sum = 0.0;
            for (int k = 0; k < num_dimensions; k++) {
                double diff = data_point[i][k] - centre[k];
                sum += diff * diff;
            }
            norms[i][j] = sqrt(sum);
        }
    }
}

4. 调整update_degree_of_membership的并行粒度 #

去掉collapse(2)，改为每个线程处理完整的一个数据点的所有成员度更新，保证连续内存访问，避免伪共享：

double update_degree_of_membership() {
    double max_diff = 0.0;

    precompute_norms();

    #pragma omp parallel for reduction(max:max_diff) schedule(static)
    for (int i = 0; i < num_data_points; i++) {
        double diff_i = 0.0;
        // 每个线程处理一个数据点的所有聚类成员度更新
        for (int j = 0; j < num_clusters; j++) {
            double sum = 0.0;
            double norm_ij = norms[i][j];
            for (int k = 0; k < num_clusters; k++) {
                sum += pow(norm_ij / norms[i][k], 2.0 / (fuzziness - 1));
            }
            double new_uij = 1.0 / sum;
            double diff = fabs(new_uij - degree_of_memb[i][j]);
            if (diff > diff_i) {
                diff_i = diff;
            }
            degree_of_memb[i][j] = new_uij;
        }
        if (diff_i > max_diff) {
            max_diff = diff_i;
        }
    }
    return max_diff;
}

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