C语言pthreads多线程n体模拟出现堆损坏及线程异常问题

我正在用C语言编写n个天体的引力运动模拟程序，单线程版本运行正常，现在尝试用POSIX pthreads库实现多线程版本。程序把天体初始数据存在全局指针数组data中，启动12个线程（基于6核处理器），每个线程负责一组天体的相互作用计算，包括力、加速度、速度、位置更新，以及碰撞合并逻辑。

线程函数代码 #

void* calculate_step(void*index_val) {
    int index = * (int *)index_val;
    long double x_dist;
    long double y_dist;
    long double distance;
    long double force;

    for (int i = 0; i < (rows/nthreads); ++i) {  //遍历当前线程分配的所有天体
        data[i+index][X_FORCE] = 0; //重置所有力为0
        data[i+index][Y_FORCE] = 0;
        data[i+index][X_ACCEL] = 0;
        data[i+index][X_ACCEL] = 0; //重复赋值X_ACCEL，Y_ACCEL未初始化
        for (int j = 0; j < rows; ++j) { //遍历所有可能的天体对
            if (i != j && data[j][DELETED] != 1 && data[i+index][DELETED] != 1) { //排除自身和已删除的天体
                x_dist = data[j][X_POS] - data[i+index][X_POS];
                y_dist = data[j][Y_POS] - data[i+index][X_POS]; //错误：应为Y_POS减Y_POS
                distance = sqrtl(powl(x_dist, 2) + powl(y_dist, 2));
                if (distance > data[i+index][RAD] + data[j][RAD]) {
                    force = G * data[i+index][MASS] * data[j][MASS] /
                            powl(distance, 2); //计算非碰撞天体对的加速度、速度、位置
                    data[i+index][X_FORCE] += force * (x_dist / distance);
                    data[i+index][Y_FORCE] += force * (y_dist / distance);
                    data[i+index][X_ACCEL] = data[i+index][X_FORCE]/data[i+index][MASS];
                    data[i+index][X_VEL] += data[i+index][X_ACCEL]*dt;
                    data[i+index][X_POS] += data[i+index][X_VEL]*dt;
                    data[i+index][Y_ACCEL] = data[i+index][Y_FORCE]/data[i+index][MASS];
                    data[i+index][Y_VEL] += data[i+index][Y_ACCEL]*dt;
                    data[i+index][Y_POS] += data[i+index][Y_VEL]*dt;
                }
                else{
                    if (data[i+index][MASS] < data[j][MASS]) {
                        int temp;
                        temp = i;
                        i = j;
                        j = temp; //逻辑错误：i是线程局部索引，j是全局索引，交换后会越界
                    } //动量守恒
                    data[i+index][X_VEL] = (data[i+index][X_VEL] * data[i+index][MASS] + data[j][X_VEL] * data[j][MASS])/(data[i+index][MASS] + data[i+index][MASS]); //分母错误：应为两个天体质量之和
                    data[i+index][Y_VEL] = (data[i+index][Y_VEL] * data[i+index][MASS] + data[j][Y_VEL] * data[j][MASS])/(data[i+index][MASS] + data[i+index][MASS]);
                    //质心位置守恒
                    data[i+index][X_POS] = (data[i+index][X_POS] * data[i+index][MASS] + data[j][X_POS] * data[j][MASS])/(data[i+index][MASS] + data[i+index][MASS]);
                    data[i+index][Y_POS] = (data[i+index][Y_POS] * data[i+index][MASS] + data[j][Y_POS] * data[j][MASS])/(data[i+index][MASS] + data[i+index][MASS]);
                    //质量守恒
                    data[i+index][MASS] += data[j][MASS];
                    //按体积比例增加半径
                    data[i+index][RAD] = powl(powl(data[i+index][RAD], 3) + powl(data[j][RAD], 3), ((long double) 1 / (long double) 3));
                    data[j][DELETED] = 1;
                    data[j][MASS] = 0;
                    data[j][RAD] = 0;
                }
            }
        }

    }
    return NULL;
}

主循环代码 #

int main() {
    pthread_t *thread_array; //线程数组指针
    long *thread_ids;

    short num_obj;
    short sim_time;
    printf("Number of objects to simulate: \n");
    scanf("%hd", &num_obj);
    num_obj = num_obj - num_obj%12;
    printf("Timespan of the simulation: \n");
    scanf("%hd", &sim_time);
    printf("Length of time steps: \n");
    scanf("%f", &dt);
    printf("Relative complexity score: %.2f\n", (((float)sim_time/dt)*((float)(num_obj^2)))/1000); //错误：num_obj^2是按位异或，不是平方

    thread_array = malloc(nthreads*sizeof(pthread_t));
    thread_ids = malloc(nthreads*sizeof(long));

    populate(num_obj);

    int index;
    for (int i = 0; i < nthreads; ++i) { //空循环，无意义

    }
    time_t start = time(NULL);
    print_data();
    for (int i = 0; i < (int)((float)sim_time/dt); ++i) { //模拟主循环
        for (int j = 0; j < nthreads; ++j) {
            index = j*(rows/nthreads);
            thread_ids[j] = j;
            pthread_create(&thread_array[j], NULL, calculate_step, &index); //所有线程共享同一个index，存在竞态条件

        }
        for (int j = 0; j < nthreads; ++j) {
            pthread_join(thread_array[j], NULL);
            //pthread_exit(NULL); //主线程无需调用，会提前退出
        }
    }
    time_t end = time(NULL) - start;

    printf("\n");
    print_data();
    printf("Took %zu seconds to simulate %d frames with %d objects initially, now %d objects.\n", end, (int)((float)sim_time/dt), num_obj, rows);
}

运行错误 #

Number of objects to simulate:
36
Timespan of the simulation:
10
Length of time steps:
0.01
Relative complexity score: 38.00

Process finished with exit code -1073740940 (0xC0000374)

错误码0xC0000374表示堆损坏，调试模式下能运行但常规编译不行，且第一个线程负责的数据块未更新，线程创建后未执行预期操作。

1. 堆损坏的核心原因 #

（1）线程创建时的竞态条件 #

主循环中所有线程共享同一个index变量的地址，线程启动时间不确定，可能多个线程读取到相同或错误的index值，导致访问错误的数组索引，越界写入堆内存，直接损坏堆结构。

修复：为每个线程分配独立的索引存储，用thread_ids传递起始索引：

for (int j = 0; j < nthreads; ++j) {
    thread_ids[j] = j*(rows/nthreads); //存储当前线程的起始索引
    pthread_create(&thread_array[j], NULL, calculate_step, &thread_ids[j]);
}

同时修改线程函数的参数读取：

int index = *(long *)index_val; //匹配thread_ids的long类型

（2）数组越界写入 #

碰撞处理中错误交换局部索引i和全局索引j，导致data[i+index]访问超出线程负责的数组范围，破坏堆内存。

修复：取消索引交换，直接基于原索引处理合并逻辑，确保合并到质量更大的天体上：

if (data[i+index][MASS] < data[j][MASS]) {
    // 将当前天体合并到j对应的天体，标记当前天体为删除
    int target_idx = j;
    int curr_idx = i+index;
    data[target_idx][X_VEL] = (data[curr_idx][X_VEL] * data[curr_idx][MASS] + data[target_idx][X_VEL] * data[target_idx][MASS])/(data[curr_idx][MASS] + data[target_idx][MASS]);
    data[target_idx][Y_VEL] = (data[curr_idx][Y_VEL] * data[curr_idx][MASS] + data[target_idx][Y_VEL] * data[target_idx][MASS])/(data[curr_idx][MASS] + data[target_idx][MASS]);
    data[target_idx][X_POS] = (data[curr_idx][X_POS] * data[curr_idx][MASS] + data[target_idx][X_POS] * data[target_idx][MASS])/(data[curr_idx][MASS] + data[target_idx][MASS]);
    data[target_idx][Y_POS] = (data[curr_idx][Y_POS] * data[curr_idx][MASS] + data[target_idx][Y_POS] * data[target_idx][MASS])/(data[curr_idx][MASS] + data[target_idx][MASS]);
    data[target_idx][MASS] += data[curr_idx][MASS];
    data[target_idx][RAD] = powl(powl(data[target_idx][RAD], 3) + powl(data[curr_idx][RAD], 3), 1.0/3.0);
    data[curr_idx][DELETED] = 1;
    data[curr_idx][MASS] = 0;
    data[curr_idx][RAD] = 0;
    break; //跳过当前天体后续计算
} else {
    // 原合并逻辑，修正分母为两个天体质量之和
    data[i+index][X_VEL] = (data[i+index][X_VEL] * data[i+index][MASS] + data[j][X_VEL] * data[j][MASS])/(data[i+index][MASS] + data[j][MASS]);
    // 其他物理量计算同理修正分母
}

（3）物理计算中的分母错误 #

碰撞合并时重复使用当前天体的质量作为分母，可能导致数值溢出，间接引发堆内存损坏，需修正为两个天体的质量之和。

2. 第一个线程未工作的原因 #

线程创建时的竞态条件导致第一个线程读取到后续线程的index值，比如主线程快速修改index为1*(rows/nthreads)，第一个线程启动后读取到的不是0，而是第二个线程的起始索引，导致它处理错误的数据块，第一个线程的任务被跳过。

3. 其他关键问题 #

（1）多线程数据竞争 #

多个线程同时修改data数组中的同一个天体数据，会导致计算结果混乱，甚至内存访问冲突。

修复：采用双缓冲模式，每个时间步先将所有天体的新状态计算到临时数组，所有线程计算完成后再同步到主数组；或对每个天体的访问加互斥锁（注意锁会影响性能）。

（2）计算逻辑错误 #

• y_dist计算错误：应为data[j][Y_POS] - data[i+index][Y_POS]
• 复杂度计算中的num_obj^2改为num_obj * num_obj
• 线程函数中补充data[i+index][Y_ACCEL] = 0;的初始化

（3）内存泄漏 #

主函数中thread_array和thread_ids分配后未释放，需在程序末尾添加：

free(thread_array);
free(thread_ids);

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