基于BLAS、Fortran与OpenMP的大矩阵加权求和性能优化问询
2026-8-19
大规模速度场加权求和的高效实现优化问题
问题背景
核心操作是将N个速度场各自乘以高斯随机变量后求和为单个速度场，数学表达式为：
$$U(x,y,z) = \sum_{n=1}^{N} \phi_{n}(x,y,z) \xi_{n}$$
其中$\phi$为速度场，$U$为结果。数据规模庞大：每个$\phi$维度为$(100,100,30)$，共720个，因此需要探索高效实现方案提升速度。
测试的两种数据存储方式
- 将所有$\phi$堆叠为$(100,100,30*720)$的三维矩阵
- 将每个$\phi$重塑为向量后拼接成$(10010030,720)$的二维矩阵
第二种方式通过BLAS的dgemv实现矩阵-向量乘积，效率更高。但针对$(300000,720)$的矩阵，尝试用OpenMP手动拆分矩阵并行调用dgemv，即使仅2线程也比串行BLAS操作耗时更长。编译命令：gfortran -O2 -o main.out MWE.F90 -llapack -lblas -fopenmp
后续新增转置矩阵的并行方案，但仍比串行BLAS慢，与预期不符。完整测试代码如下：
PROGRAM new_modes IMPLICIT NONE INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) interface function OMP_get_wtime() INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) real(wp) :: OMP_get_wtime end function OMP_get_wtime function OMP_get_thread_num() integer(kind = 4) :: OMP_get_thread_num end function OMP_get_thread_num end interface INTEGER, PARAMETER :: nn_tlu_nmod = 720 INTEGER, PARAMETER :: jpi = 100! 32 INTEGER, PARAMETER :: jpj = 100! 17 INTEGER, PARAMETER :: jpk = 30! 31 INTEGER, PARAMETER :: np = jpk * jpj * jpi ! REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: A REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:,:) :: B REAL(wp), POINTER, DIMENSION(:,:) :: sub_B REAL(wp), POINTER, DIMENSION(:) :: sub_B_vec REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: U_ref REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: U REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:) :: vec_U REAL(wp), POINTER, DIMENSION(:) :: sub_vec_U REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:) :: tcoeff REAL(wp), POINTER, DIMENSION(:) :: tcoef_pt ! REAL(wp) :: tic, toc INTEGER :: ji, jj, jk, jm INTEGER :: m_idx INTEGER :: cnt ! INTEGER :: ierr ! REAL(wp) :: ddot ! !-------------------------------------------------------------------------------------! ! Declaration for the Parallel parameters !-------------------------------------------------------------------------------------! INTEGER, PARAMETER :: num_threads=1 INTEGER, SAVE :: myID !$OMP THREADPRIVATE(myID) !-------------------------------------------------------------------------------------! ! Declaration for the row scheduling !-------------------------------------------------------------------------------------! INTEGER :: sub_n INTEGER :: istart,inext INTEGER, ALLOCATABLE :: work_sharing(:) WRITE(6,*) 'Number of threads employed in the solver:', num_threads ALLOCATE( A(jpi,jpj,jpk * nn_tlu_nmod), stat=ierr ) ALLOCATE( B(jpi*jpj*jpk, nn_tlu_nmod), stat=ierr ) ALLOCATE( Bt( nn_tlu_nmod, jpi*jpj*jpk), stat=ierr ) ALLOCATE( tcoeff(nn_tlu_nmod), stat=ierr ) ALLOCATE( U(jpi,jpj,jpk), & & U_ref(jpi,jpj,jpk), & & vec_U(jpi*jpj*jpk), stat=ierr ) cnt = 0 DO jm = 1, nn_tlu_nmod ! m_idx = ( jm - 1 ) * jpk ! DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi cnt = cnt + 1 A(ji,jj,m_idx + jk) = cnt * jm END DO END DO END DO B(:,jm) = RESHAPE(A(:,:,m_idx + 1: m_idx + jpk), (/jpi*jpj*jpk/) ) Bt(jm,:) = RESHAPE(A(:,:,m_idx + 1: m_idx + jpk), (/jpi*jpj*jpk/) ) END DO tcoeff = 1._wp ! First option, normal matrix summation (already tested faster than ! with explicit do loop) tic = omp_get_wtime() U = 0._wp ! DO jm = 1, nn_tlu_nmod ! Define zero-th indexed mode m_idx = ( jm - 1 ) * jpk ! U(:,:,:) = U(:,:,:) + A(:,:,m_idx + 1 : m_idx + jpk ) * tcoeff(jm) ! ENDDO ! toc = omp_get_wtime() print '("Time, operation without for loops = ",f13.7," seconds.")', toc-tic U_ref = U ! Second option, exploit a different shape of the matrix in order to use BLAS ! (Usage of pointers speeds-up a lot) tic = omp_get_wtime() U = 0._wp sub_B => B(:, :) sub_n = np sub_vec_U => vec_U(:) tcoef_pt => tcoeff(:) ! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy ) CALL DGEMV( 'n', sub_n, nn_tlu_nmod, 1._wp, sub_B, sub_n, tcoef_pt, 1, 0._wp, sub_vec_U, 1 ) U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) ) ! toc = omp_get_wtime() print '("Time, 2D matrix with BLAS = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref)) ! Third option, Divide the job in a static way between M threads, use pointers to point at the ! portion of the matrix assigned to each processor and use BLAS tic = omp_get_wtime() vec_U = 0._wp !$OMP PARALLEL PRIVATE(istart, sub_n) myID = omp_get_thread_num()+1 istart = work_sharing(myID) sub_n = work_sharing(myID+1)-istart ! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy ) CALL DGEMV( 'n', sub_n, nn_tlu_nmod, 1._wp, B(istart,1), np, tcoeff, 1, 0._wp, vec_U(istart), 1 ) !$OMP END PARALLEL U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) ) toc = omp_get_wtime() print '("Time, Static scheduling BLAS OpenMP = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref)) ! Fourth option, run a parallel Do and use DDOT tic = omp_get_wtime() vec_U = 0._wp !$OMP PARALLEL PRIVATE(ji) !$OMP DO SCHEDULE (STATIC) DO ji=1,np ! res = DDOT( n, x, inc_x, y, inc_y ) vec_U(ji) = DDOT( nn_tlu_nmod, B(ji, :), 1, tcoeff, 1 ) END DO !$OMP END DO !$OMP END PARALLEL U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) ) toc = omp_get_wtime() print '("Time, Parallelized for loop DDOT = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref)) ! Fifth option, Divide the job in a static way between M threads, but using a transposed matrix tic = omp_get_wtime() vec_U = 0._wp !$OMP PARALLEL PRIVATE(istart, sub_n) myID = omp_get_thread_num()+1 istart = work_sharing(myID) sub_n = work_sharing(myID+1) - istart ! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy ) CALL DGEMV( 't', nn_tlu_nmod, sub_n, 1._wp, Bt(1,istart), nn_tlu_nmod, tcoeff, 1, 0._wp, vec_U(istart), 1 ) !$OMP END PARALLEL U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) ) toc = omp_get_wtime() print '("Time, Static scheduling tran OpenMP = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref)) END PROGRAM
问题分析与优化建议
1. 手动拆分dgemv并行效率低下的原因
- BLAS内置并行冲突：绝大多数优化版BLAS（如OpenBLAS、MKL）本身已经实现多线程并行，手动用OpenMP拆分矩阵再调用
dgemv会引发线程嵌套竞争——BLAS内部线程与手动添加的OpenMP线程抢占CPU资源，额外增加调度开销，反而拖慢速度。 - 缓存局部性恶化：拆分矩阵后，每个线程处理的子矩阵无法充分利用CPU缓存，而串行BLAS会针对整个矩阵做缓存优化，内存访问效率更高。
2. 转置矩阵方案慢的原因
- 内存布局不连续：Fortran默认采用列主序存储，原矩阵$B$的列是连续内存块；转置后的$Bt$逻辑上是行优先，但物理存储仍为列主序，调用
DGEMV('t'...)时会访问不连续的内存，缓存命中率大幅下降，性能恶化。 - 同样存在线程竞争问题：若BLAS本身开启多线程，手动OpenMP并行依然会导致调度冲突。
3. 优化方案
方案一：直接使用多线程优化版BLAS
- 替换系统默认BLAS为OpenBLAS或Intel MKL，这类库会自动根据CPU核心数优化
dgemv的并行度，无需手动拆分矩阵。 - 编译示例（OpenBLAS）：
gfortran -O2 -o main.out MWE.F90 -lopenblas -fopenmp - 运行前通过环境变量设置线程数（根据CPU核心数调整）：
export OMP_NUM_THREADS=4
方案二：调整并行策略，避免线程嵌套
- 先禁用BLAS内部多线程（如
export OPENBLAS_NUM_THREADS=1或export MKL_NUM_THREADS=1），再用OpenMP并行处理行维度的DDOT计算。同时优化内存访问连续性：! 转置后Bt为(720, np)，列主序下Bt(:,ji)是连续内存块 !$OMP PARALLEL PRIVATE(ji) !$OMP DO SCHEDULE(STATIC) DO ji=1,np vec_U(ji) = DDOT(nn_tlu_nmod, Bt(:,ji), 1, tcoeff, 1) END DO !$OMP END DO !$OMP END PARALLEL
方案三：利用Fortran数组操作优化
- 直接使用
matmul替代手动循环或拆分：将$\phi$拼接为$(np,720)$的矩阵，执行vec_U = matmul(B, tcoeff)，优化版编译器（如gfortran -O3）会自动调用BLAS或做并行优化。
4. 关键注意事项
- 优先依赖优化后的BLAS/LAPACK库，避免手动实现并行矩阵操作，这类库的优化程度远高于自定义代码。
- 时刻注意Fortran的列主序存储特性，保证内存访问连续，这是提升数值计算性能的核心。
- 测试手动并行方案时，必须关闭BLAS内部多线程，避免线程竞争干扰测试结果。
内容的提问来源于stack exchange，提问作者Kimala
相关产品推荐
相关产品推荐
