Fortran基于OpenMP实现矩阵乘法计算结果错误问题求解

问题定位 #

你的代码存在两处核心错误导致并行版本计算结果异常：

• 局部数组私有问题：子程序multi中的mm1、mm2是局部变量，在OpenMP并行区域内调用子程序时，默认会为每个线程创建独立的私有副本。第一个!$OMP DO拆分循环后，每个线程仅会给自己的副本中对应分片的位置赋值，其余位置为未初始化的垃圾值，后续乘法计算自然出错。
• 打印逻辑位置错误：打印语句放在!$OMP PARALLEL块内部，会导致每个线程都执行一次打印，输出多份重复、甚至可能未完成计算的结果。

修改方案 #

1. 将mm1、mm2调整为并行区域的共享变量，可将其移至主程序声明，在并行块中显式指定为共享属性，作为参数传入子程序即可，避免每个线程生成独立副本。
2. 将打印语句移到并行区域外部，等所有线程计算完成后再输出结果。
3. 可选优化：9x9矩阵的一维数组填充计算量极小，也可直接放在并行区域外串行执行，完全避免共享变量的潜在竞争问题。

修正后的可运行代码 #

program main
    implicit none
    double precision,dimension(0:8,0:8)::MM,MMT,MTEMP
    ! 将mm1、mm2移至主程序声明，作为共享变量使用
    double precision,dimension(0:80)::mm1,mm2

    MM=reshape((/ 1 ,   1   ,   1   ,   1   ,   1   ,   1   ,   1   ,   1   ,   1   ,   &
    -4  ,   -1  ,   -1  ,   -1  ,   -1  ,   2   ,   2   ,   2   ,   2   ,   &
    4   ,   -2  ,   -2  ,   -2  ,   -2  ,   1   ,   1   ,   1   ,   1   ,   &
    0   ,   1   ,   0   ,   -1  ,   0   ,   1   ,   -1  ,   -1  ,   1   ,   &
    0   ,   -2  ,   0   ,   2   ,   0   ,   1   ,   -1  ,   -1  ,   1   ,   &
    0   ,   0   ,   1   ,   0   ,   -1  ,   1   ,   1   ,   -1  ,   -1  ,   &
    0   ,   0   ,   -2  ,   0   ,   2   ,   1   ,   1   ,   -1  ,   -1  ,   &
    0   ,   1   ,   -1  ,   1   ,   -1  ,   0   ,   0   ,   0   ,   0   ,   &
    0   ,   0   ,   0   ,   0   ,   0   ,   1   ,   -1  ,   1   ,   -1  /),shape(MM))

    MMT=1.d0/36.d0 * reshape  ((/ 4 ,   -4  ,   4   ,   0   ,   0   ,   0   ,   0   ,   0   ,   0   ,   &
                    4   ,   -1  ,   -2  ,   6   ,   -6  ,   0   ,   0   ,   9   ,   0   ,   &
                    4   ,   -1  ,   -2  ,   0   ,   0   ,   6   ,   -6  ,   -9  ,   0   ,   &
                    4   ,   -1  ,   -2  ,   -6  ,   6   ,   0   ,   0   ,   9   ,   0   ,   &
                    4   ,   -1  ,   -2  ,   0   ,   0   ,   -6  ,   6   ,   -9  ,   0   ,   &
                    4   ,   2   ,   1   ,   6   ,   3   ,   6   ,   3   ,   0   ,   9   ,   &
                    4   ,   2   ,   1   ,   -6  ,   -3  ,   6   ,   3   ,   0   ,   -9  ,   &
                    4   ,   2   ,   1   ,   -6  ,   -3  ,   -6  ,   -3  ,   0   ,   9   ,   &
                    4   ,   2   ,   1   ,   6   ,   3   ,   -6  ,   -3  ,   0   ,   -9  /),shape(MMT))

    ! 显式指定并行区域变量属性，避免隐式规则出错
    !$OMP PARALLEL DEFAULT(NONE) SHARED(MM,MMT,MTEMP,mm1,mm2)
        call multi(mm,mmt,mtemp,mm1,mm2)
    !$OMP END PARALLEL

    ! 打印移到并行块外，仅执行一次
    PRINT*,MTEMP

end program main

subroutine multi(m1,m2,m3,mm1,mm2)
    implicit none
    double precision,dimension(0:8,0:8)::m1,m2,m3
    double precision,dimension(0:80)::mm1,mm2
    DOUBLE PRECISION::TEMP
    integer::i,j,k
    ! 显式指定循环变量私有
    !$OMP DO PRIVATE(i,j)
    do j=0,8
        do i=0,8
            mm1(j*9+i)=m1(i,j)
            mm2(i*9+j)=m2(i,j)
        enddo
    enddo
    !$OMP ENDDO
    !$OMP DO PRIVATE(TEMP,I,J,K)
    do j=0,8
        do i=0,8
            temp=0
            do k=0,8
                temp=temp+mm1(j*9+k)*mm2(i*9+k)
            enddo
            m3(i,j)=temp
        enddo
    enddo
    !$OMP ENDDO
    return
end subroutine multi

注：当前测试用矩阵尺寸极小，并行调度开销远大于计算收益，仅作为OpenMP语法学习使用。

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