24小时热门版块排行榜    

CyRhmU.jpeg
查看: 388  |  回复: 1
当前主题已经存档。
【有奖交流】积极回复本帖子,参与交流,就有机会分得作者 futurehero 的 100 个金币

futurehero

新虫 (正式写手)

[交流] 【求助】诚求助并行计算中的一个问题,谢谢大家!

这个程序大致如下,我想用并行计算解决“do循环中包含积分”的计算问题。如果将并行的线程设为1,则和串行计算的结果一直,但是线程大于1时,计算的结果就不正确了。小弟刚刚接触并行计算两天,做了很多种方案都发现线程大于1时计算的结果都不对,请诸位大侠帮助,我想得到争取的并行计算结果。如果将程序中的call omp_set_num_threads(4 )线程改为call omp_set_num_threads(1 )就和串行结果是一致的。(可能您只需要看一下并行区就可以发现问题所在了)
程序如下:

program main  
implicit none
            
integer:: Num, i, p, j, k ,var, r ,Y
real(kind=8):: D, epsl21, epsl22, T, eta, Gama, GamaL, GamaR, epsl0, VL, VR

parameter(D=1.0_8, Num=10001, Gama=0.01_8, GamaL=0.005_8, GamaR=0.005_8, epsl0=-0.1_8, VL=0.0_8, VR= 0.0_8)

real(kind=8):: pi, Dc, Hc, fL(Num), fR(Num), f0(Num), epsl(Num), epsl11, epsl12, dIdV(Num)
complex(kind=8):: sigma0, FG11(Num), FG12(Num), FG21(Num), FG22(Num)

complex(kind=8):: G011(Num), G012(Num), G021(Num), G022(Num), G11(Num), G12(Num), G21(Num), G22(Num)
complex(kind=8):: DDOS(Num), DOS, D11, D12, D21, D22
complex(kind=8):: n11, n12, n21, n22

complex(kind=8):: delta1211(Num), delta2111(Num), delta2211(Num)
complex(kind=8):: delta1112(Num), delta2112(Num), delta2212(Num)

complex(kind=8):: delta1121(Num), delta1221(Num), delta2221(Num)
complex(kind=8):: delta1122(Num), delta1222(Num), delta2122(Num)

complex(kind=8)::   FZB11(Num),   FZB12(Num),   FZB21(Num),   FZB22(Num)

complex(kind=8)::   B1211(Num),   B2111(Num),   B2211(Num)  
complex(kind=8)::   FB1211(Num),  FB2111(Num),  FB2211(Num)

complex(kind=8)::   B1112(Num),   B2112(Num),   B2212(Num)  
complex(kind=8)::  FB1112(Num),  FB2112(Num),  FB2212(Num)  

complex(kind=8)::  B1121(Num),   B1221(Num),   B2221(Num)
complex(kind=8)::  FB1121(Num),  FB1221(Num),  FB2221(Num)

complex(kind=8)::   B1222(Num),   B2122(Num),   B1122(Num)
complex(kind=8)::  FB1222(Num),  FB2122(Num),  FB1122(Num)

real(kind=8)::    m, Field, SOC, Tk   

Dc=2*D / (Num-1)  
eta =2.0*Dc                       
pi=4.0_8*datan(1.0_8)  

m=0.00_8  
T=0.0_8   
SOC=0.00000_8   
Field=0.00000_8  
epsl11=epsl0 - SOC/2.0_8 - m*Field - Field
epsl12=epsl0 + SOC/2.0_8 + m*Field - Field  
epsl21=epsl0 + SOC/2.0_8 - m*Field + Field
epsl22=epsl0 - SOC/2.0_8 + m*Field + Field
do i=1, Num, 1
   epsl(i)=-D + (i-1)  *Dc                                                                           
end do
do i=1, Num, 1        
   if(  epsl(i)           <  VL       )  fL(i)=1.0_8         
   if(  epsl(i)           >  VL       )  fL(i)=0.0_8               
   if(  dabs(epsl(i)-VL)  <  0.01_8*Dc ) fL(i)=0.5_8   
   if(  epsl(i)           <  VR       )  fR(i)=1.0_8                     
   if(  epsl(i)           >  VR       )  fR(i)=0.0_8         
   if(  dabs(epsl(i)-VR)  <  0.01_8*Dc ) fR(i)=0.5_8            
end do                                                                                             
f0=(    fL*GamaL  +  fR*GamaR    ) /Gama         
delta1211=dcmplx(  epsl12 - epsl11,      eta  ) - epsl                 
delta2111=dcmplx(  epsl21 - epsl11,      eta  ) - epsl              
delta2211=dcmplx(  epsl22 - epsl11,      eta  ) - epsl                 
      
delta1112=dcmplx(  epsl11 - epsl12,      eta  ) - epsl         
delta2112=dcmplx(  epsl21 - epsl12,      eta  ) - epsl         
delta2212=dcmplx(  epsl22 - epsl12,      eta  ) - epsl  
delta1121=dcmplx(  epsl11 - epsl21,      eta  ) - epsl         
delta1221=dcmplx(  epsl12 - epsl21,      eta  ) - epsl      
delta2221=dcmplx(  epsl22 - epsl21,      eta  ) - epsl            
delta1122=dcmplx(  epsl11 - epsl22,      eta  ) - epsl                  
delta1222=dcmplx(  epsl12 - epsl22,      eta  ) - epsl         
delta2122=dcmplx(  epsl21 - epsl22,      eta  ) - epsl      
!================================
sigma0=dcmplx(0.0_8, -Gama)                                                                        
G011  =epsl-epsl11-sigma0                                                              
G012  =epsl-epsl12-sigma0                                                              
G021  =epsl-epsl21-sigma0                                                              
G022  =epsl-epsl22-sigma0                                                              
!********************************************
call omp_set_num_threads(4 )
!$OMP parallel                                                                        
!$OMP do private (Y)
do Y=1, Num, 1                                                                              
   FB1211=f0/( epsl(Y) + delta1211  )                                                                    
   FB2111=f0/( epsl(Y) + delta2111  )                                                                    
   FB2211=f0/( epsl(Y) + delta2211  )                                                                    
   FB1112=f0/( epsl(Y) + delta1112  )                                                                    
   FB2112=f0/( epsl(Y) + delta2112  )                                                                    
   FB2212=f0/( epsl(Y) + delta2212  )                                                                    
   FB1121=f0/( epsl(Y) + delta1121  )                                                                    
   FB1221=f0/( epsl(Y) + delta1221  )                                                                    
   FB2221=f0/( epsl(Y) + delta2221  )                                                                    
   FB1122=f0/( epsl(Y) + delta1122  )                                                                    
   FB1222=f0/( epsl(Y) + delta1222  )                                                                    
   FB2122=f0/( epsl(Y) + delta2122  )                                                                    
call INTEGRAL(Num, Dc,  FB1211, B1211(Y))
call INTEGRAL(Num, Dc,  FB2111, B2111(Y))   
call INTEGRAL(Num, Dc,  FB2211, B2211(Y))
call INTEGRAL(Num, Dc,  FB1112, B1112(Y))
call INTEGRAL(Num, Dc,  FB2112, B2112(Y))
call  INTEGRAL(Num, Dc,  FB2212, B2212(Y))
call  INTEGRAL(Num, Dc,  FB1121, B1121(Y))
call  INTEGRAL(Num, Dc,  FB1221, B1221(Y))
call  INTEGRAL(Num, Dc,  FB2221, B2221(Y))
call  INTEGRAL(Num, Dc,  FB1122, B1122(Y))      
call  INTEGRAL(Num, Dc,  FB1222, B1222(Y))   
call  INTEGRAL(Num, Dc,  FB2122, B2122(Y))
end do                                                                             
!$OMP end do
!$OMP end parallel  
!******************************************
B1211=B1211*Gama/pi
B2111=B2111*Gama/pi
B2211=B2211*Gama/pi   
                     
B1112=B1112*Gama/pi  
B2112=B2112*Gama/pi   
B2212=B2212*Gama/pi   
                     
B1121=B1121*Gama/pi   
B1221=B1221*Gama/pi   
B2221=B2221*Gama/pi   
                        
B1122=B1122*Gama/pi     
B1222=B1222*Gama/pi     
B2122=B2122*Gama/pi      
G11= (   1.0_8 - 3.0_8*0.245_8    )/(   G011 - B1211-B2111-B2211   )   
G12= (   1.0_8 - 3.0_8*0.245_8    )/(   G012 - B1112-B2112-B2212   )
G21= (   1.0_8 - 3.0_8*0.245_8    )/(   G021 - B1121-B1221-B2221   )
G22= (   1.0_8 - 3.0_8*0.245_8    )/(   G022 - B1122-B1222-B2122   )

FG11=-f0*G11/pi              
FG12=-f0*G12/pi   
FG21=-f0*G21/pi         
FG22=-f0*G22/pi                 
call INTEGRAL(Num, Dc, FG11, n11)  
call INTEGRAL(Num, Dc, FG12, n12)
call INTEGRAL(Num, Dc, FG21, n21)
call INTEGRAL(Num, Dc, FG22, n22)

write(*,*) "n11=", dimag(n11), "n12=", dimag(n12)
write(*,*) "n21=", dimag(n21), "n22=", dimag(n22)

stop
end program main

subroutine INTEGRAL(Num, Dc, Y, answer)                                                                  
implicit none                                                         
integer::  Num                                                               
real(kind=8):: Dc                                                         
complex(kind=8):: Y(Num), answer
integer:: i                                                                           
complex:: s1                                                                     
                                                                                               
s1=(Y(1)+ Y(Num))/2                                                                              
do i=2, Num-1, 1                                                                       
   s1=s1+ Y(i)                                                                                   
end do                                                                        
answer=s1*Dc                                                                       
                                                                    
return                                                                                       
end subroutine INTEGRAL  

stop
end program main

[ Last edited by futurehero on 2010-3-4 at 23:09 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

futurehero

新虫 (正式写手)

已经解决了!错误原因在于不同的线程会同时将中间变量写到同一个内存单元FB1211,造成结果错误.只需用将再为每个线程定义不同的变量就解决问题了.
2楼2010-03-20 16:08:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 futurehero 的主题更新
普通表情 高级回复(可上传附件)
信息提示
请填处理意见