openmp矩阵乘法
The openmp matrix multiplication
我尝试编写一个基于 Openmp 的矩阵乘法代码。矩阵mm与矩阵mmt的乘积为对角矩阵且等于1。我尝试正常计算和 Openmp。正常结果是正确的,但是 Openmp 结果是错误的。我觉得应该是跟Openmp利用率有关。
program main
implicit none
double precision,dimension(0:8,0:8)::MM,MMT,MTEMP
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
call multi(mm,mmt,mtemp)
PRINT*,MTEMP
!$OMP END PARALLEL
endprogram main
subroutine multi(m1,m2,m3)
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
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
endsubroutine
我在下面给出了正常版本,如果你计算这个,结果是对角 1 矩阵。
program main
implicit none
double precision,dimension(0:8,0:8)::MM,MMT,MTEMP
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))
call multi(mm,mmt,mtemp)
PRINT*,MTEMP
endprogram main
subroutine multi(m1,m2,m3)
double precision,dimension(0:8,0:8)::m1,m2,m3
double precision,dimension(0:80)::mm1,mm2
DOUBLE PRECISION::TEMP
integer::i,j,k
do j=0,8
do i=0,8
mm1(j*9+i)=m1(i,j)
mm2(i*9+j)=m2(i,j)
enddo
enddo
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
return
endsubroutine
要解决您的问题,一种解决方案是将并行区域放在 multi
子例程内。此代码给出与串行代码相同的结果:
subroutine multi(m1,m2,m3)
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 PARALLEL
!$OMP DO
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
!$OMP END PARALLEL
return
注意工作量很小,所以可能不会比串口代码快。另请注意,如果以后不需要 mm1
和 mm2
数组,则不必计算它们:
subroutine multi(m1,m2,m3)
double precision,dimension(0:8,0:8)::m1,m2,m3
DOUBLE PRECISION::TEMP
integer::i,j,k
!$OMP PARALLEL
!$OMP DO PRIVATE(TEMP,I,J,K)
do j=0,8
do i=0,8
temp=0
do k=0,8
temp=temp+m1(k,j)*m2(i,k)
enddo
m3(i,j)=temp
enddo
enddo
!$OMP ENDDO
!$OMP END PARALLEL
return
我尝试编写一个基于 Openmp 的矩阵乘法代码。矩阵mm与矩阵mmt的乘积为对角矩阵且等于1。我尝试正常计算和 Openmp。正常结果是正确的,但是 Openmp 结果是错误的。我觉得应该是跟Openmp利用率有关。
program main
implicit none
double precision,dimension(0:8,0:8)::MM,MMT,MTEMP
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
call multi(mm,mmt,mtemp)
PRINT*,MTEMP
!$OMP END PARALLEL
endprogram main
subroutine multi(m1,m2,m3)
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
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
endsubroutine
我在下面给出了正常版本,如果你计算这个,结果是对角 1 矩阵。
program main
implicit none
double precision,dimension(0:8,0:8)::MM,MMT,MTEMP
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))
call multi(mm,mmt,mtemp)
PRINT*,MTEMP
endprogram main
subroutine multi(m1,m2,m3)
double precision,dimension(0:8,0:8)::m1,m2,m3
double precision,dimension(0:80)::mm1,mm2
DOUBLE PRECISION::TEMP
integer::i,j,k
do j=0,8
do i=0,8
mm1(j*9+i)=m1(i,j)
mm2(i*9+j)=m2(i,j)
enddo
enddo
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
return
endsubroutine
要解决您的问题,一种解决方案是将并行区域放在 multi
子例程内。此代码给出与串行代码相同的结果:
subroutine multi(m1,m2,m3)
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 PARALLEL
!$OMP DO
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
!$OMP END PARALLEL
return
注意工作量很小,所以可能不会比串口代码快。另请注意,如果以后不需要 mm1
和 mm2
数组,则不必计算它们:
subroutine multi(m1,m2,m3)
double precision,dimension(0:8,0:8)::m1,m2,m3
DOUBLE PRECISION::TEMP
integer::i,j,k
!$OMP PARALLEL
!$OMP DO PRIVATE(TEMP,I,J,K)
do j=0,8
do i=0,8
temp=0
do k=0,8
temp=temp+m1(k,j)*m2(i,k)
enddo
m3(i,j)=temp
enddo
enddo
!$OMP ENDDO
!$OMP END PARALLEL
return