1

I am currently testing OpenMP in a big loop in my FORTRAN code. The code is part of a simulation module which is called from a VB.NET user interface; this interface also does the timing measurements. So I start a simulation, and at the end the software shows me how long it took (I only write this to show that for timing measurements I don't use wtime or cpu_time). Now when I repeatedly start a simulation with my parallelized loop, I always get different simulation times, reaching, in one example, from 1min30sec to almost 3min! The results are always correct.

I tried different schedules for the loop (static, guided, dynamic), I tried to calculate the chunks that are assigned to each thread manually (do i=1,N -> do i=i_start,i_end), I tried to change the number of threads taking part in the calculation of the loop - with no change of the situation. When I remove the OpenMP directives from the code this does not occur, so they must be the reason for this behavior. My machine is a quadcore Intel Xeon(R) CPU X3470 @2.93GHz with Win7 installed. I tried to run the program with both enabled and disabled multithreading (in the bios), however, this also didn't change anything.

Do you have any ideas what could go wrong? A web search for a situation like this showed that similar behavior occured in test environments of other programmers as well, however a solution / reason has never been mentioned. Thanks in advance for your thoughts.

Martin

EDIT

Here's the code:

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(n,k,nk,i,j,l,List,Vx,Vz,cS,AE1,RootCh,Ec1,Ec2,Ec3,FcE,GcE,VxE,VzE,SMuL1,SMuL2) &
!$OMP PRIVATE(W1,W2,W3,Wx,Wz,S,i1,j1,AcE,j2,ic,iB,iBound,i2) &
!$OMP FIRSTPRIVATE(NumSEL) REDUCTION(-:Cum0,Cum1) REDUCTION(+:CumR)
    DO n=1, NumEl
!       Loop on subelements    
        DO k=1, Elements(n)%nCorners-2
            nk = (k-1) * 3
            NumSEL=NumSEL+1
!            
            i=Elements(n)%KX(1)
            j=Elements(n)%KX(k+1)
            l=Elements(n)%KX(k+2)
            List(1)=i
            List(2)=j
            List(3)=l
!            
            IF(Level == NLevel) THEN
                Vx(1)=Nodes(i)%VxO
                Vx(2)=Nodes(j)%VxO
                Vx(3)=Nodes(l)%VxO
                Vz(1)=Nodes(i)%VzO
                Vz(2)=Nodes(j)%VzO
                Vz(3)=Nodes(l)%VzO
            ELSE
                Vx(1)=Nodes(i)%VxN
                Vx(2)=Nodes(j)%VxN
                Vx(3)=Nodes(l)%VxN
                Vz(1)=Nodes(i)%VzN
                Vz(2)=Nodes(j)%VzN
                Vz(3)=Nodes(l)%VzN
            END IF
!            
            cS=cBound(sol,5)
            cS=(MIN(cS,Nodes(i)%Conc(sol))+MIN(cS,Nodes(j)%Conc(sol))+MIN(cS,Nodes(l)%Conc(sol)))/3.0D0       

            AE1=Elements(n)%xMul(k)*Elements(n)%Area(k)*dt*Eps
            RootCh=AE1*cS*(Nodes(i)%Sink+Nodes(j)%Sink+Nodes(l)%Sink)/3.0D0
            Cum0=Cum0-AE1*(Nodes(i)%Gc1+Nodes(j)%Gc1+Nodes(l)%Gc1)/3.0D0
            Cum1=Cum1-AE1*(Nodes(i)%Fc1+Nodes(j)%Fc1+Nodes(l)%Fc1)/3.0D0
            CumR=CumR+RootCh
            Ec1=(Nodes(i)%Dispxx+Nodes(j)%Dispxx+Nodes(l)%Dispxx)/3.0D0
            Ec2=(Nodes(i)%Dispxz+Nodes(j)%Dispxz+Nodes(l)%Dispxz)/3.0D0
            Ec3=(Nodes(i)%Dispzz+Nodes(j)%Dispzz+Nodes(l)%Dispzz)/3.0D0
!
            IF (Level == NLevel) AcE=(Nodes(i)%Ac+Nodes(j)%Ac+Nodes(l)%Ac)/3.0D0
!
            FcE=(Nodes(i)%Fc+Nodes(j)%Fc+Nodes(l)%Fc)/3.0D0
            GcE=(Nodes(i)%Gc+Nodes(j)%Gc+Nodes(l)%Gc)/3.0D0
            VxE=(Vx(1)+Vx(2)+Vx(3))/3.0D0
            VzE=(Vz(1)+Vz(2)+Vz(3))/3.0D0
            SMul1=-Elements(n)%AMul(k)
            SMul2=Elements(n)%Area(k)/20.0D0*Elements(n)%XMul(k)
!
            If (lUpw) THEN
                !W1=WeTab(1,NumSEl)
                !W2=WeTab(2,NumSEl)
                !W3=WeTab(3,NumSEl)
                W1=WeTab(1,(n-1)*(Elements(n)%nCorners-2)+k)
                W2=WeTab(2,(n-1)*(Elements(n)%nCorners-2)+k)
                W3=WeTab(3,(n-1)*(Elements(n)%nCorners-2)+k)
                Wx(1)=2.0D0*Vx(1)*(W2-W3)+Vx(2)*(W2-2.0D0*W3)+Vx(3)*(2.0D0*W2-W3)
                Wx(2)=Vx(1)*(2.0D0*W3-W1)+2.0D0*Vx(2)*(W3-W1)+Vx(3)*(W3-2.0D0*W1)
                Wx(3)=Vx(1)*(W1-2.0D0*W2)+Vx(2)*(2.0D0*W1-W2)+2.0D0*Vx(3)*(W1-W2)
                Wz(1)=2.0D0*Vz(1)*(W2-W3)+Vz(2)*(W2-2.0D0*W3)+Vz(3)*(2.0D0*W2-W3)
                Wz(2)=Vz(1)*(2.0D0*W3-W1)+2.0D0*Vz(2)*(W3-W1)+Vz(3)*(W3-2.0D0*W1)
                Wz(3)=Vz(1)*(W1-2.0D0*W2)+Vz(2)*(2.0D0*W1-W2)+2.0D0*Vz(3)*(W1-W2)
            END IF
!
            DO j1=1, 3
                i1=List(j1)
!$OMP           ATOMIC
                Nodes(i1)%F=Nodes(i1)%F+Elements(n)%GMul(k)*(GcE+Nodes(i1)%Gc/3.0D0)
                IF (Level==NLevel) then
!$OMP               ATOMIC
                    Nodes(i1)%DS=Nodes(i1)%DS+Elements(n)%GMul(k)*(Ace+Nodes(i1)%Ac/3.0D0)
                end if
                iBound=0
                IF (Nodes(i)%Kode/=0) THEN
                    BP_Loop : DO id=1, NumBP
                        IF((KXB(id)==i1) .AND. (KodCB(id) > 0)) THEN
                            iBound=1
                            EXIT BP_Loop
                        END IF
                    END DO BP_Loop
                END IF 
!
                DO j2=1, 3
                    i2=List(j2)
                    S(j1,j2)=SMul1*(Ec1*Elements(n)%dz(nk+j1)*Elements(n)%dz(nk+j2)+ &
                                    Ec3*Elements(n)%dx(nk+j1)*Elements(n)%dx(nk+j2)+ &
                                   Ec2*(Elements(n)%dz(nk+j1)*Elements(n)%dx(nk+j2)+ &
                                        Elements(n)%dx(nk+j1)*Elements(n)%dz(nk+j2)))

                    S(j1,j2)=S(j1,j2)-(Elements(n)%dz(nk+j2)/8.0D0*(VxE+Vx(j1)/3.0D0)+ & 
                                       Elements(n)%dx(nk+j2)/8.0D0*(VzE+Vz(j1)/3.0D0))*Elements(n)%xMul(k)

                    IF(lUpw) S(j1,j2)=S(j1,j2)-Elements(n)%xMul(k)* &
                                              (Elements(n)%dz(nk+j2)/40.0D0*Wx(j1)+ &
                                               Elements(n)%dx(nk+j2)/40.0D0*Wz(j1))

                    ic=1
                    IF (i1==i2) ic=2
                    S(j1,j2)=S(j1,j2)+SMul2*ic*(FcE+(Nodes(i1)%Fc+Nodes(i2)%Fc)/3.0D0)
                    IF (iBound==1) then
                        if(j2.eq.1) then
!$OMP                      ATOMIC
                                Nodes(i1)%Qc(sol)=Nodes(i1)%Qc(sol)-Eps*S(j1,j2)*Nodes(i2)%Conc(sol)-Eps*Elements(n)%GMul(k)*(GcE+Nodes(i1)%Gc/3.0D0)
                        else
!$OMP                       ATOMIC
                            Nodes(i1)%Qc(sol)=Nodes(i1)%Qc(sol)-Eps*S(j1,j2)*Nodes(i2)%Conc(sol)
                        end if
                    end if
                    IF (Level/=NLevel) THEN
!$OMP                   ATOMIC
                        B(i1)=B(i1)-alf*S(j1,j2)*Nodes(i2)%Conc(sol)
                    ELSE
                        IF (lOrt) THEN
                            CALL rFIND(i1,i2,kk,NumNP,MBandD,IAD,IADN)
                            iB=kk
                    ELSE
                        iB=MBand+i2-i1
                    END IF
!$OMP                   ATOMIC
                        A(iB,i1)=A(iB,i1)+epsi*S(j1,j2)
                    END IF
                END DO
            END DO
        END DO
    END DO
!$OMP END PARALLEL DO    
5
  • You need to provide some code! If your loop has a high frequency and there might be problems with scheduling your threads in critical sections.
    – FFox
    Jul 4, 2011 at 11:22
  • I edited my original post, it includes the code now.
    – Martin
    Jul 7, 2011 at 8:43
  • I have a question about your code. Variable kk is shared and this seems like it could be a problem (if lOrt is "true", then multiple threads could call rFind and set kk - leading to a race condition).
    – ejd
    Jul 7, 2011 at 14:30
  • Right, I see that - I changed that now, but I already tried a couple of modifications of this code, including making the variable kk private, and copying the few lines of rFIND-code directly into the loop - it didn't change anything. Seems like I falsely removed kk from private-list again.
    – Martin
    Jul 8, 2011 at 5:09
  • I think your problem is related to the bunch of "atomic" statements. Since they do not allow a parallel processing of this part, they block all other threads and effectively synchronise them. Depending on how much relative delay the OMP threads accumulate at the startup this is a performance problem or not (just as you measure it).
    – FFox
    Jul 9, 2011 at 11:32

1 Answer 1

0

If you want to check the performance in the program i would suggest you did timings in the program with the OpenMP timing functions. See OpenMP Ref. sheet. So you need to do something like:

USE omp_lib

t1 = omp_get_wtime()
! Big loop
t_final = omp_get_wtime() - t1

I some time find these to reflect the actual parallization timings better. Do you use those?

As FFox says it can simply be due to the ATOMIC statements which is delaying in different manors on each run. Remember that the threads are created at run time, so the layout of the threads may not be the same for each run.

With such a loop i would try to see if you could gain speed by splitting it up. Of course this is not needed if the speedup is around 2 for 2 threads. Just an idea.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.