openmp并行Fortran后每次结果不一致

? Vita 2020-04-08 10:18:39

进行并行化之后每次结果不一致，对赋值代码进行了atomic保护也不好使。有没有大佬帮忙看看怎么回事。
代码很长。也不知道有没有大佬有耐心看完....

  implicit none

  integer n,m,node,node1,node2,ni,mj,p,q,k

  integer i,j



  real(4) rcm(3),rcn(3),dist

  real r1f(3),r2f(3),r3f(3),nSm(3),rjf(3,SNum_samp),r1s(3),r2s(3),r3s(3),ris(3,SNum_samp)

  real Am,An,lj,li,rmj(3),rni(3),rtmpv(3)

  real Jrm(3,maxRWG),Irn(3,maxRWG),Jlengthm(maxRWG),Ilengthn(maxRWG)



  complex coef1,coef2

  complex Amn_E,Zji,I1ss,I2ss(3),I3ss(3),I4ss

  complex ctmp,ctmp1

    

  real,external:: re_Norm

  integer*4 time1,time2,time

  

  call SYSTEM_CLOCK(time1)

  write(*,*) "Start to compute the SS part of the impedance matrix, on the way..."

  

  allocate(SparseValue(Num_side,Num_side))

  SparseValue=0

      

  coef1=cmplx(0.0,omega*MU0/(16.0*PI))

  coef2=cmplx(0.0,1.0/(4.0*PI*omega*EPS0))

 

!$OMP parallel do default(none) shared(Num_Triangle,triangle_centroid,triangle_area,point_cor,triangle_point,triangle_rwgs,side_index,side_length,SparseValue,threshold_R,coef1,coef2) private(n,m,node,node1,node2,ni,mj,i,j,rcm,rcn,dist,r1f,r2f,r3f,rjf,r1s,r2s,r3s,ris,Am,An,lj,li,rmj,rni,Jrm,Irn,Jlengthm,Ilengthn,Amn_E,ctmp,I1ss,I2ss,I3ss,I4ss)

  do m=1,Num_Triangle    

    rcm=triangle_centroid(:,m)

	   Am=triangle_area(m)

 

    r1f=point_cor(:,triangle_point(1,m))

    r2f=point_cor(:,triangle_point(2,m))

    r3f=point_cor(:,triangle_point(3,m))



	   do k=1,SNum_samp

      rjf(:,k)=epsilS(k)*r1f+etaS(k)*r2f+xiS(k)*r3f

    enddo

 

    do mj=1,triangle_rwgs(m)%N_rwg

      j=triangle_rwgs(m)%rwg(mj)

        

      if(side_index(3,j)==m) then

        !flagm=1.0

        Jlengthm(mj)=side_length(j)

      else  !if(side_index(4,j)==m) then

        !flagm=-1.0

        Jlengthm(mj)=-side_length(j)

      endif

        

   	  node1=side_index(1,j); node2=side_index(2,j)

      do p=1,3

        node=triangle_point(p,m)

        if((node/=node1).and.(node/=node2)) exit

      enddo         

      Jrm(:,mj)=point_cor(:,node)   	          

    enddo 

 

    do n=1,Num_Triangle

      rcn=triangle_centroid(:,n)

      An=triangle_area(n)



      r1s=point_cor(:,triangle_point(1,n))     !source triangle coordinates

      r2s=point_cor(:,triangle_point(2,n))

      r3s=point_cor(:,triangle_point(3,n))



      do k=1,SNum_samp

	       ris(:,k)=epsilS(k)*r1s+etaS(k)*r2s+xiS(k)*r3s

      enddo

                    

      dist=re_Norm(rcm-rcn)          

      

      call Comp_eltsI(I1ss,I2ss,I3ss,I4ss,dist>=threshold_R,An,rjf,ris,r1s,r2s,r3s,SNum_samp,SNum_samp,weightS,weightS)

       

      do ni=1,triangle_rwgs(n)%N_rwg

        i=triangle_rwgs(n)%rwg(ni)

		      if(side_index(3,i)==n) then

	         !flagn=1.0

          Ilengthn(ni)=side_length(i)

        else  !if(side_index(4,i)==n) then

	         !flagn=-1.0

          Ilengthn(ni)=-side_length(i)

	       endif

              

	       node1=side_index(1,i); node2=side_index(2,i)

	       do q=1,3

          node=triangle_point(q,n)

          if((node/=node1).and.(node/=node2)) exit

	       enddo

        Irn(:,ni)=point_cor(:,node)

      enddo  

      

      do mj=1,triangle_rwgs(m)%N_rwg

        j=triangle_rwgs(m)%rwg(mj)

        rmj=Jrm(:,mj)       

        lj=Jlengthm(mj) 

 

        ctmp=coef1*(I4ss-dot_product(rmj,I3ss))-coef2*I1ss

         

        do ni=1,triangle_rwgs(n)%N_rwg

          i=triangle_rwgs(n)%rwg(ni)

                    

          rni=Irn(:,ni)               

          li=Ilengthn(ni)

                    

          Amn_E=lj*li*(coef1*(dot_product(rmj,rni)*I1ss-dot_product(rni,I2ss))+ctmp)



!$OMP atomic

          SparseValue(i,j)=SparseValue(i,j)+Amn_E     

         

	       enddo

      enddo

    enddo



	 enddo

!$OMP end parallel do