Fortran Kmeans代码OpenMP并行优化：多线程聚类寻优与线程报告实现

核心并行逻辑改造 #

你的需求是让每个线程独立执行完整聚类流程，不能在原代码的双重迭代循环中加并行——那是数据拆分式并行，和你的需求不符。正确的做法是把整个聚类流程包裹在OpenMP并行区域内，让每个线程独立初始化、迭代、计算代价。

关键步骤： #

1. 定义线程结果存储类型：

! 存储每个线程的聚类结果、代价、耗时
type :: ThreadResult
    real :: total_cost       ! 代价函数值
    double precision :: run_time  ! 线程运行耗时
    real, allocatable :: centroids(:,:)  ! 该线程最终聚类中心
end type ThreadResult

type(ThreadResult), allocatable :: thread_results(:)
integer :: num_threads

2. 初始化线程结果数组：

! 获取当前OpenMP线程数
!$OMP PARALLEL
    num_threads = omp_get_num_threads()
!$OMP END PARALLEL
allocate(thread_results(num_threads))

3. 并行执行独立聚类：

! 声明线程私有变量：每个线程独立拥有自己的聚类中心、代价数组、迭代计数器等
real, allocatable :: local_centroids(:,:), local_cluster_cost(:)
real :: local_total_cost, prev_cost
integer :: local_iter
double precision :: t_start, t_end

!$OMP PARALLEL PRIVATE(local_centroids, local_cluster_cost, local_total_cost, prev_cost, local_iter, t_start, t_end)
    ! 每个线程独立初始化随机种子（避免所有线程初始聚类中心完全相同）
    call random_seed()
    ! 自定义初始化聚类中心函数（复用你的串行代码逻辑）
    call init_centroids(data, k, local_centroids)

    ! 记录线程开始时间
    t_start = omp_get_wtime()

    ! 线程独立执行完整Kmeans流程（和你的串行代码逻辑完全一致）
    prev_cost = huge(1.0)
    do local_iter = 1, max_iter
        ! 分配样本到聚类中心（复用你的串行代码）
        call assign_samples_to_clusters(data, local_centroids, local_cluster_cost)
        local_total_cost = sum(local_cluster_cost)
        
        ! 判断收敛
        if (abs(local_total_cost - prev_cost) < convergence_threshold) exit
        prev_cost = local_total_cost

        ! 更新聚类中心（复用你的串行代码）
        call update_cluster_centers(data, local_centroids)
    end do

    ! 记录线程结束时间
    t_end = omp_get_wtime()

    ! 将线程结果写入共享数组（必须用临界区避免数据竞争）
    !$OMP CRITICAL
        thread_results(omp_get_thread_num() + 1)%total_cost = local_total_cost
        thread_results(omp_get_thread_num() + 1)%run_time = t_end - t_start
        ! 复制聚类中心到结果结构体
        if (allocated(thread_results(omp_get_thread_num() + 1)%centroids)) deallocate(thread_results(omp_get_thread_num() + 1)%centroids)
        allocate(thread_results(omp_get_thread_num() + 1)%centroids(size(local_centroids,1), size(local_centroids,2)))
        thread_results(omp_get_thread_num() + 1)%centroids = local_centroids
    !$OMP END CRITICAL
!$OMP END PARALLEL

最优结果筛选 #

遍历所有线程的结果，筛选出代价函数最小的；若代价相同，则选运行时间最短的：

integer :: best_thread_idx
best_thread_idx = 1

do i = 2, num_threads
    ! 优先比较代价函数
    if (thread_results(i)%total_cost < thread_results(best_thread_idx)%total_cost) then
        best_thread_idx = i
    ! 代价相同时比较运行时间
    else if (abs(thread_results(i)%total_cost - thread_results(best_thread_idx)%total_cost) < 1e-6) then
        if (thread_results(i)%run_time < thread_results(best_thread_idx)%run_time) then
            best_thread_idx = i
        end if
    end if
end do

! 输出最优结果
print *, "=== 最优聚类结果 ==="
print *, "线程ID: ", best_thread_idx - 1  ! OpenMP线程ID从0开始计数
print *, "代价函数值: ", thread_results(best_thread_idx)%total_cost
print *, "运行时间: ", thread_results(best_thread_idx)%run_time, "秒"

生成多线程运行报告 #

直接遍历thread_results数组，格式化输出每个线程的运行数据：

print *, new_line('a')//"=== 各线程运行报告 ==="
do i = 1, num_threads
    print *, "线程ID: ", i - 1
    print *, "  代价函数值: ", thread_results(i)%total_cost
    print *, "  运行耗时: ", thread_results(i)%run_time, "秒"
    print *
end do

常见错误排查（针对你之前的循环并行错误） #

• 错误逻辑：在聚类的迭代循环中加并行，是让线程分工处理样本分配/中心更新，和“每个线程跑完整流程”的需求完全不符，必然导致数据冲突。
• 常见坑点：
  • 未将局部变量标记为PRIVATE：线程间共享变量会导致数据混乱。
  • 所有线程用相同初始聚类中心：结果完全一致，并行失去意义，必须每个线程独立初始化随机种子。
  • 写入共享结果时未加临界区：多个线程同时写数组会导致数据损坏。
  • 未正确分配/释放局部数组：容易引发内存泄漏或访问越界。

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